solve the Poisson equation using the Red & Black-SOR method

 from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes 3D
import numpy as np
import itertools
import csv
 # condition setting
delta_L = 1.0
delta = 1.0
n_inter = 0
nx = 100
ny = 100
xmin = 0
xmax = 100 * delta
ymin = 0
ymax = 100* delta

# even odd determination
default_even(Ryouiki)
  if Ryouiki %2 == 0:
    red
  else:
    black


convergence_criterion = 1.0e-10
 
dx=(xmax-xmin)/(nx-1)
dy=(ymax-ymin)/(ny-1)
 
# initial state
phi=np.zeros(ny,nx))
electric=np.zeros(ny,nx))
x = np.linspace (xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)

# for SOR method
aa_recta=0.5*(np.cos(np.pi/nx)+np.cos(np.pi/ny))#
omega_SOR_recta=2/(1+np.sqrt(1-aa_recta**2)# Optimal acceleration parameters for rectangular areas
print(" omega_SOR_rect=", omega_SOR_recta)
 
# charge
eps0 = 1
charge=np.zeros(ny,nx))
Q1 = 1
Q2 = Q1/delta**2

charge [50,50] = Q2

#Red

conv_check=[]
while delta>convergence_criterion:
    phi_in=phi.copy()
    if n_inter%50 == 0:
      print("interaction No=", n_inter, "delta=", delta)
    conv_check.append([n_inter::1, delta])
    for i in range (nx-1):
        for jin range (ny-1):
          if i == xmin or i == xmax or j == ymax or j == ymin:
            phi[i,j] = 0
          else:
            phi[i,j] = phi[i,j] + omega_SOR_recta* (omega_SOR_recta*1/2*)

    delta = np.max (abs(phi-phi_in))
    
    n_inter+=1


print("The number of totalization=",n_inter)
print("data_points=",nx*ny)
    
body=(phi)
header=('Electrical potential of 2D array')
 # Open File in Write Mode
with open('Deni1.csv', 'w') as f:
 
  writer=csv.writer(f)#Create writer object
  writer.writerows(header)# Write the header
  writer.writerows(body)#Write content

body2=(x)
header2=('Value1 of Axis')
with open('jikunoatai1.csv', 'w') as f1:
  writer=csv.writer(f1)
  writer.writerow (header2)
  writer.writerow (body2)


defplot2D(x,y,phi):
    config=plt.figure(figsize=(11,7), dpi=100)
    ax=fig.gca(projection='3d')
    X,Y = np.meshgrid(x,y)
    surf=ax.plot_surface(X,Y,phi[:],rstride=1,cstride=1,cmap=cm.viridis,linewidth=0,antialized=False)
    ax.view_init(30,225)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('phi')
 
plot2D(x,y,phi)
plt.show()


# Generate Object
config=plt.figure(figsize=(11,7), dpi=100)
figure
 
xmin = 0
xmax = 100
ymin = 0
ymax = 100
 
x = np.linspace (xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)
 
X,Y = np.meshgrid(x,y)
plt.contourf(X,Y,phi,alpha=0.5,cmap=cm.viridis) 
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')

#potential calculation
print("Keisan V")
for i in range (1, 6):
    print("x=", i*delta_L, "V=", phi[50+i,50] - phi[50+int(1/delta_L),50])

# electric field calculation
L = 99
Ey = np.zeros ([L, L])
for i in range (L):
    for jin range (L):
        Ey[i,j]=-(phi[i+1,j]-phi[i-1,j])/2/delta_L
print("Keisan E")
for i in range (1, 6):
    print("x=", i*delta_L, "E=", Ey[50+i,50])