Potentials and Fields
Problem 5.7
背景：
In regions of space that do not contain any electric charges, the electric potential obeys Laplace's equation:
1.gif
We can get approximation as：
2.gif
For two dimensional problem：
3.gif
Jacobi method：
4.gif
Gauss-Seidel method：
5.gif
SOR method：
6.gif
Potentials and Fields
In regions of space that do not contain any electric charges, the electric potential obeys Laplace's equation:
1.gif
We can get approximation as：
2.gif
For two dimensional problem：
3.gif
Jacobi method：
4.gif
Gauss-Seidel method：
5.gif
SOR method：
6.gif
正文：
import numpy as np
from pylab import *
from math import *
import mpl_toolkits.mplot3d
global error
error=1e-5
def Jacobi(L):
V0=[[0 for i in range(L)]for j in range(L)]#i represents x, j represents y
a=int(2*(L-1)/5)
b=int(3*(L-1)/5)
for i in [a]:
for j in range(a,b+1):
V0[j][i]=1.0# j,i
for i in [b]:
for j in range(a,b+1):
V0[j][i]=-1.0
VV=[]
VV.append(V0)
s=0
dx=0.1
#iteration
while 1:
VV.append(V0)
for i in range(1,L-1):
for j in range(1,L-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 [a]:
for j in range(a,b+1):
VV[s+1][j][i]=1.0
for i in [b]:
for j in range(a,b+1):
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,L-1):
for j in range(1,L-1):
dV=dV+abs(dVV[i][j])
s=s+1
print dV
if dV<error*(L-1)**2 and s>1:
#if dV<error and s>10:
break
return s
def GS(L):
V0=[[0 for i in range(L)]for j in range(L)]#i represents x, j represents y
a=int(2*(L-1)/5)
b=int(3*(L-1)/5)
for i in [a]:
for j in range(a,b+1):
V0[j][i]=1.0# j,i
for i in [b]:
for j in range(a,b+1):
V0[j][i]=-1.0
VV=[]
VV.append(V0)
s=0
dx=0.1
#iteration
while 1:
VV.append(V0)
for i in range(1,L-1):
for j in range(1,L-1):
VV[s+1][i][j]=(VV[s][i+1][j]+VV[s+1][i-1][j]+VV[s][i][j+1]+VV[s+1][i][j-1])/4.0
for i in [a]:
for j in range(a,b+1):
VV[s+1][j][i]=1.0
for i in [b]:
for j in range(a,b+1):
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,L-1):
for j in range(1,L-1):
dV=dV+abs(dVV[i][j])
s=s+1
if dV<error*(L-1)**2 and s>1:
#if dV<error and s>1:
break
return s
def SOR(L):
a=int(2*(L-1)/5)
b=int(3*(L-1)/5)
V0=[[0 for i in range(L)]for j in range(L)]#i represents x, j represents y
for i in [a]:
for j in range(a,b+1):
V0[j][i]=1.0# j,i
for i in [b]:
for j in range(a,b+1):
V0[j][i]=-1.0
VV=[]
VV.append(V0)
alpha=2.0/(1+pi/L)
s=0
dx=0.1
#iteration
while 1:
VV.append(V0)
for i in range(1,L-1):
for j in range(1,L-1):
VV[s+1][j][i]=(VV[s][j+1][i]+VV[s+1][j-1][i]+VV[s][j][i+1]+VV[s+1][j][i-1])/4.0
if i==a and j>a-1 and j<b+1:
VV[s+1][j][i]=1.0
if i==b and j>a-1 and j<b+1:
VV[s+1][j][i]=-1.0
VV[s+1][j][i]=alpha*(VV[s+1][j][i]-VV[s][j][i])+VV[s][j][i]
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,L-1):
for j in range(1,L-1):
dV=dV+abs(dVV[i][j])
print dV,L
s=s+1
if dV<error*(L-1)**2 and s>1:
break
return s
L=[]
NJ=[]
NGS=[]
NSOR=[]
f=open('problem5.7.txt','w')
print >> f,'L','J','GS','SOR'
for i in range(6,61,5):
J=Jacobi(i)
G=GS(i)
S=SOR(i)
L.append(i)
NJ.append(J)
NGS.append(G)
NSOR.append(S)
print >> f,i,J,G,S
f.close()
plot(L,NJ)
plot(L,NGS)
plot(L,NSOR)
scatter(L,NJ)
scatter(L,NGS)
scatter(L,NSOR)
legend(('Jacobi method','GS method','SOR method'),'upper left')
title('3 different methods',fontsize=15)
xlabel('L')
ylim(0,1000)
ylabel('N')
savefig('different methods.png')
show()
7.gif
从图中可以看出，SOR方法所需迭代次数最少，且随着计算区域尺度的增加仅仅是线性增大的，另两种方式则是二次方增加的。
