第十二次作业

---

一、背景
1.Laplace'equation

$$\triangledown ^2V=0$$
It can be written in a rectangular coordinate system as$\frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2}+\frac{\partial^2 V}{\partial z^2}=0$.
2.Poisson'equation
$$\triangledown ^2V=f$$
It is a Non-homogeneous partial differential equations.
二、代码
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)]for k in range(21)]#i represents x, j represents y
VV=[]
VV.append(V0)
s=0
dx=0.1
p0=[[[0 for i in range(21)]for j in range(21)]for k in range(21)]
for i in [10]:
    for j in [10]:
        for k in [10]:
            p0[k][j][i]=1/dx**3
p=[]
p.append(p0)
while True:
    VV.append(V0)
    p.append(p0)
    for m in range(1,len(V0)-1):
        for j in range(1,len(V0[1])-1):
            for i in range(1,len(V0[1][1])-1):
                VV[s+1][m][j][i]=(VV[s][m][j][i+1]+VV[s][m][j][i-1]+VV[s][m][j+1][i]+VV[s][m][j-1][i]+VV[s][m+1][j][i]+VV[s][m-1][j][i])/6.0+p[s][m][j][i]*dx**2/6.0
    VV[s]=np.array(VV[s])
    VV[s+1]=np.array(VV[s+1])
    dVV=VV[s+1]-VV[s]
    dV=0
    for k in range(len(V0)):
        for j in range(1,len(V0[1])-1):
            for i in range(1,len(V0[1][1])-1):
                dV=dV+abs(dVV[k][j][i])
    s=s+1
    if dV<0.0001 and s>1:
        break
print s
V=np.array(VV[-1][10])
Ex=np.array(V0[10])
Ey=np.array(V0[10])
for j in range(1,len(V0[10])-1):
    for i in range(1,len(V0[10][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)
figure(figsize=[8,8])
x=np.arange(-1.0,1.01,dx)
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(-0.5,0.5)
xlabel('x')
ylim(-0.5,0.5)
ylabel('y')
title('Equipotential lines around a point charge in 3D')
fig = figure(figsize=[8,8])
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')
figure(figsize=[8,8])
Q=quiver(X,Y,Ex,Ey,scale=100)
xlim(-1,1)
xlabel('x')
ylim(-1,1)
ylabel('y')
title('Field lines around a point charge in 3D')
show()

2

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)]for k in range(21)]#i represents x, j represents y
VV=[]
VV.append(V0)
s=0
dx=0.1
p0=[[[0 for i in range(21)]for j in range(21)]for k in range(21)]
for i in [10]:
    for j in [10]:
        for k in [10]:
            p0[k][j][i]=1/dx**3
p=[]
p.append(p0)
while True:
    VV.append(V0)
    p.append(p0)
    for m in range(1,len(V0)-1):
        for j in range(1,len(V0[1])-1):
            for i in range(1,len(V0[1][1])-1):
                VV[s+1][m][j][i]=(VV[s][m][j][i+1]+VV[s][m][j][i-1]+VV[s][m][j+1][i]+VV[s][m][j-1][i]+VV[s][m+1][j][i]+VV[s][m-1][j][i])/6.0+p[s][m][j][i]*dx**2/6.0
    VV[s]=np.array(VV[s])
    VV[s+1]=np.array(VV[s+1])
    dVV=VV[s+1]-VV[s]
    dV=0
    for k in range(len(V0)):
        for j in range(1,len(V0[1])-1):
            for i in range(1,len(V0[1][1])-1):
                dV=dV+abs(dVV[k][j][i])
    s=s+1
    if dV<0.0001 and s>1:
        break
print s
V=[]
x=np.linspace(0.1,1.0,10)
X=np.linspace(0.01,1.00,99)
Y=[]
for i in range(len(X)):
    Y.append(1/(4*pi*X[i]))
for i in range(11,21):
    V.append(VV[-1][10][10][i])
figure(figsize=[8,8])
plot(X,Y)
scatter(x,V)
title('Potential near a point charge')
legend(('V versus x(y=0,z=0)'),'upper right')
xlabel('x')
ylabel('V')
xlim(0,1.1)
ylim(0,1.1)
show()

三、致谢
感谢李鹏对我本次作业的帮助！