Plotting potentials of 2d parallel plate capacitor in 0V box. Plots coming out wrong?

Question

I'm trying to plot the potential, equipotential, and electric field lines of a parallel plate capacitor in a box of 0V 2 dimensionally. Only the voltage on the plates and the distance between the plates are given. When I initially made the plot I assumed that the charge density was zero since the plates are equal and opposite. But the plot came out suspicious (potentials are skewed to the left and vectors point towards the positive V plate. Aren't the vectors suppose to point towards the negative?). Problem Linked here

I'm using the Gauss-Seidel method with overrelaxation to calculate the potentials.

L=100
a=1
e0=8.85e-12
p=0
w=0.9
iterations=0
tolerance=1e-6
dvmax=100
dec=4
vpot=np.zeros((L,L))
vpot[20:80,20:21]=1
vpot[20:80,80:81]=-1
xy = np.linspace(0,L-1,L)
xydec=xy[::dec]
while dvmax>tolerance:
      iterations+=1
      dvmax=0
      for i in range(1,L-1):
          for j in range(1,L-1):
              vnew=(vpot[i-1,j]+vpot[i+1,j]+vpot[i,j-1]+vpot[i,j+1]+(p*(a**2)/e0))/4
              dv=(1+w)*(vnew-vpot[i,j])
              if abs(dv)>dvmax:
                 dvmax=abs(dv)
              vpot[i,j]+=dv
print(f'Number of iterations: {iterations}')
plt.imshow(vpot, origin='lower')
plt.plasma()
plt.colorbar()
plt.title(f'Plot of Positive Potentials and Centered Charge Density with Tolerance {tolerance}')
plt.figure()
cv=plt.contour(vpot,15,colors='k', linestyles='solid')
plt.clabel(cv, fontsize=8, inline=True)
plt.title('Contour of Potentials')
plt.figure()
Ex=np.zeros((L,L))
Ey=np.zeros((L,L))
for j in range(1,L-1):
          for k in range(1,L-1):
              if k==0:
                 Ex[j,k] = -(vpot[j,k+1] - vpot[j,k])/a
              elif k==L-1:
                 Ex[j,k] = -(vpot[j,k] - vpot[j,k-1])/a
              else:
                 Ex[j,k] = -(vpot[j,k+1] - vpot[j,k-1])/(2*a)
              if j==0:
                 Ey[j,k] = -(vpot[j+1,k] - vpot[j,k])/a
              elif j==L-1:
                 Ey[j,k] = -(vpot[j,k] - vpot[j-1,k])/a
              else:
                 Ey[j,k] = -(vpot[j+1,k] - vpot[j-1,k])/(2*a)
Exdec = Ex[::dec,::dec]
Eydec = Ey[::dec,::dec]
plt.quiver(xydec,xydec,Exdec,Eydec)
plt.title('Electric Field Vectors')
plt.show()

I'm honestly stuck on what is wrong other than my assumption that the charge density rho, or p in my code, is zero. But I don't know how to calculate the charge considering I'm only given the voltage on each plate and their separation distance.

Plot of potentials and equipotential

Plot of field vectors

Answer 1

You are updating ALL internal potentials in the Gauss-Seidel iteration. You are supposed to keep those that you set as 1 and -1 constant. Probably the easiest way of dealing with this is to mark those points for non-update:

vpot=np.zeros((L,L));   vfixed=np.zeros((L,L))
vpot[20:80,20:21]=1;    vfixed[20:80,20:21] = 1
vpot[20:80,80:81]=-1;   vfixed[20:80,80:81] = 1

Then, in your loop you can just skip any points for which vfixed[] is non-zero:

  for i in range(1,L-1):
      for j in range(1,L-1):
          if vfixed[i,j] != 0: continue

Then your output looks reasonable: