Codes

Code 1

import numpy as np
from pylab import *
from math import *
import mpl_toolkits.mplot3d
V0=[[0 for i in range(21)]for j in range(21)]#i represents x, j represents y
for i in [6]:
    for j in range(6,15):
        V0[j][i]=1.0# j,i
for i in [15]:
    for j in range(6,15):
        V0[j][i]=-1.0
#check
#for j in range(len(V0)):
#    for i in range(len(V0[1])):
#        print V0[i][j], 
#    print
VV=[]
VV.append(V0)
s=0
dx=0.1
#iteration
#for k in range(100):
while 1:
    VV.append(V0)
    for i in range(1,len(V0)-1):
        for j in range(1,len(V0[1])-1):
            VV[s+1][i][j]=(VV[s][i+1][j]+VV[s][i-1][j]+VV[s][i][j+1]+VV[s][i][j-1])/4.0
    for i in [6]:
        for j in range(6,15):
            VV[s+1][j][i]=1.0
    for i in [15]:
        for j in range(6,15):
            VV[s+1][j][i]=-1.0
    VV[s]=np.array(VV[s])
    VV[s+1]=np.array(VV[s+1])
    dVV=VV[s+1]-VV[s]
    dV=0
    for i in range(1,len(V0)-1):
        for j in range(1,len(V0[1])-1):
            dV=dV+abs(dVV[i][j])
    #print dV 
    s=s+1
    if dV<0.0001 and s>10:
        break
print s
V=array(VV[-1])
Ex=array(V0)
Ey=array(V0)
for j in range(1,len(V0)-1):
    for i in range(1,len(V0[1])-1):
        Ex[j][i]=-(V[j][i+1]-V[j][i-1])/(2*dx)
        Ey[j][i]=-(V[j+1][i]-V[j-1][i])/(2*dx)
#for i in range(len(V)):
#    for j in range(len(V[1])):
#        print V[i][j], 
#    print
figure(figsize=[8,8])
x=np.arange(-1.0,1.01,dx)#-1.0,-0.9,...,1.0
y=np.arange(-1.0,1.01,dx)
X,Y=np.meshgrid(x,y)
CS = contour(X,Y,V,15)
clabel(CS, inline=1, fontsize=10)
xlim(-1,1)
xlabel('x')
ylim(-1,1)
ylabel('y')
title('contours of electric potential')
savefig('electric potential J.png')
fig = figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, V,rstride=1, cstride=1,linewidth=0)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('V')
title('electric potential')
savefig('electric potential J 3D.png')
figure(figsize=[8,8])
Q=quiver(X,Y,Ex,Ey,scale=100)
xlim(-1,1)
xlabel('x')
ylim(-1,1)
ylabel('y')
title('electric field')
savefig('electric field J.png')
show()

Code 2

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
class jacobi:
    """
    Jacobi method
    """
    def __init__(self):
        pass
    def initialization(self,initial_file):
        itxt= open(initial_file)
        self.lattice_in = []
        self.lattice_out = []
        for lines in itxt.readlines():
            lines=lines.replace("\n","").split(",")
            #print lines[0].split(" ")
            self.lattice_in.append(lines[0].split(" "))
            self.lattice_out.append(lines[0].split(" "))
        itxt.close()
        #print self.lattice_in
        #print self.initial_lattice
        for i in range(len(self.lattice_in)):
            for j in range(len(self.lattice_in[i])):
                self.lattice_in[i][j] = float(self.lattice_in[i][j])
                self.lattice_out[i][j] = float(self.lattice_out[i][j])
        return 0
    def evolution(self):
        delta = 10
        self.N = 0
        delta_record = []
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
            delta_record.append(delta)
            for i in range(len(self.lattice_in)):
                for j in range(len(self.lattice_in[i])):
                    self.lattice_in[i][j] = self.lattice_out[i][j]
            self.N+= 1
        plt.plot(delta_record,label = 'jacobi method')
        print self.N
        print len(self.lattice_in)
        return 0
    def electric_field(self):
        self.ex = np.zeros([len(self.lattice_in),len(self.lattice_in)])
        self.ey = np.zeros([len(self.lattice_in),len(self.lattice_in)])
        deltax = 0.3333
        deltay = 0.3333
        for i in range(1,len(self.lattice_in) - 1):
            for j in range(1,len(self.lattice_in[0]) - 1):
                self.ey[i][j] = (-(self.lattice_in[i + 1][j] - self.lattice_in[i - 1][j])/(2*deltax))
                self.ex[i][j] = (-(self.lattice_in[i][j + 1] - self.lattice_in[i][j - 1])/(2*deltay))
        return 0
    def plot_field(self):
        self.electric_field()
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        #plt.figure(figsize = (8,8))
        plt.quiver(X, Y, self.ex, self.ey, pivot = 'mid', units='width')
        plt.title('Electric field between two metal plates')
        #plt.show()
        return 0
    def plot_contour(self): 
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        #plt.figure(figsize = (8,8))
        CS = plt.contour(X, Y, self.lattice_in, 13)
        plt.clabel(CS, inline=1, fontsize=10)
        plt.title('Equipotential lines')
        #plt.show()
        return 0
    def plot_surface(self):
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        surf = ax.plot_surface(X, Y, self.lattice_in, rstride=1, cstride=1,cmap = cm.coolwarm_r,
                       linewidth=0, antialiased=False)
        ax.set_zlim(-1.01, 1.01)
        ax.zaxis.set_major_locator(LinearLocator(10))
        ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        fig.colorbar(surf, shrink=0.5, aspect=5)
class SOR(jacobi):
    def evolution(self):
        delta = 10
        self.N = 0
        delta_record = []
        alpha = 2/(1+np.pi/len(self.lattice_in))
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
                        self.lattice_in[i][j] = alpha * (self.lattice_out[i][j] - self.lattice_in[i][j]) + self.lattice_in[i][j]
            #print self.lattice_out
            delta_record.append(delta)
            self.N+= 1
        plt.plot(delta_record, label ='SOR algorithm')
        print self.N
        print len(self.lattice_in)
        #print self.lattice_in
        return 0    
class SOR_2(jacobi):
    def evolution(self,alpha):
        delta = 10
        self.N = 0
        delta_record = []
        #alpha = 2/(1+np.pi/len(self.lattice_in))
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
                        self.lattice_in[i][j] = alpha * (self.lattice_out[i][j] - self.lattice_in[i][j]) + self.lattice_in[i][j]
            #print self.lattice_out
            delta_record.append(delta)
            self.N+= 1
        #plt.plot(delta_record, label ='SOR algorithm')
        print self.N
        #print len(self.lattice_in)
        #print self.lattice_in
        return 0  
#A = jacobi()
A = SOR()
A.initialization('initial_capacitor_26.txt')
A.evolution()
#plt.figure(figsize = (8,8))
plt.figure(figsize = (16,8))
plt.subplot(121)
A.plot_contour()
plt.subplot(122)
A.plot_field()
#A.plot_surface()
plt.savefig('chapter5_5.7_26.png', dpi = 256)
plt.show()
### ----------------------
'''A = jacobi()
B = SOR()
A.initialization('initial_capacitor_26.txt')
B.initialization('initial_capacitor_26.txt')
plt.figure(figsize = (8,8))
plt.xlabel('The number of iteration, N',fontsize = 14)
plt.ylabel(r'$\Delta V$', fontsize = 14)
A.evolution()
B.evolution()
plt.legend()
plt.savefig('chapter5_5.7_alpha.png', dpi = 256)
plt.show()'''
### ----------------------
'''A = SOR_2()
alpha = 1.00
for i in range(20):
    A.initialization('initial_capacitor_26.txt')
    A.evolution(alpha)
    alpha+=0.05'''