1

我正在尝试修改旧代码以解决新问题，但我遇到了一些问题。我有这段代码可以使用松弛方法解决区域中任何点的潜力。

from matplotlib.pylab import *
nx=100
V=zeros((nx,nx))
V[nx-1,:]=1
err=1.0
while (err > 1e-3):
    Vold=V.copy()
    V[1:-1,1:-1]=0.25*(V[:-2,1:-1]+V[2:,1:-1]+V[1:-1,:-2]+V[1:-1,2:])
    u = (V - Vold)
    err = sqrt(sum(u*u))
imshow(V, extent=[0,1,0,1])
title('Electrostatic Potential')
colorbar()
show()

我试图求解长度为 2b 的接地立方体的体积内的电势，该立方体的电荷 Q 均匀分布在较大立方体中心的另一个立方体（边 = 2a）的体积上。其中 Q=1.5*10^-7 C, b=3m, a=1m N=40 个分区，每个分区的长度为 3/20 m。由于 Q 均匀分布在中心立方体上，因此当 a=1 m 时，较小立方体内各处的 p（电荷密度）是 Q/Vol= Q/8 C/m^3 的常数值。让 p 是一个 3d 数组给出了 p= Q/8 在 x,y,z=(N/2-(1/3)N,N/2+(1/3)N) 和 p=0 之间的其他任何地方。我到目前为止的代码是

from matplotlib.pylab import *
N=40
b=3.0
a=1.0
e=1.256*10**-6
h=(2.0*b)/N
q=1.5*10**-7
Vol=8*a**3
V=zeros((N,N,N)) 
p=zeros((N,N,N))
V[:,:,:]=0
p[(N/2-(a*N)/b):(N/2+(a*N)/b),(N/2-(a*N)/b):(N/2+(a*N)/b),(N/2-(a*N)/b):(N/2+(a*N/b)]=q/Vol
err=1.0
while (err > 1e-3):
    Vold=V.copy()
    V[1:-1,1:-1,1:-1]=(V[:-2,1:-1,1:-1]+V[2:,1:-1,1:-1]+V[1:-1,:-2,1:-1]+V[1:-1,2:,1:-1]+V[1:-1,1:-1,:-2]+V[1:-1,1:-1,2:])/6+p[1:-1,1:-1,1:-1]*h**2/(6*e)
    u = (V - Vold)
    err = sqrt(sum(u*u))

print V[N/2,N/2,N/2]

但我没有得到 V(center) ~= 1166 Volts 的预期值。而是得到 V(center) ~= 0.02 Volts。程序是否正确，我的预期值是否正确，或者计算中是否有错误？

4

0 回答 0