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()

より大きな立方体の中心にある別の立方体 (辺 = 2a) の体積に均等に分布する電荷 Q を持つ、長さ 2b の接地立方体の体積内のポテンシャルを解こうとしています。ここで、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 を 3 次元配列とすると、x,y,z=(N/2-(1/3)N,N/2+(1/3)N) の間で p= Q/8 となり、それ以外では 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.I の期待値を取得していません。代わりに V(center) ~= 0.02 Volts を取得しています。プログラムは正しく、期待値は正しくないか、計算に誤りがありますか??

4

0 に答える 0