Exercise 12

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方法所需迭代次数最少，且随着计算区域尺度的增加仅仅是线性增大的，另两种方式则是二次方增加的。