鏡像法その3

下の図の様に電荷qを置くとき鏡像法を用いて電位、電界を可視化する。電位は
                   ϕ=(q/(4(Subscript[πϵ, 0])(Sqrt[(((x-a))^2)+(((y-a))^2)+(z^2)])))+(q/(4(Subscript[πϵ, 0])(Sqrt[(((x+a))^2)+(((y+a))^2)+(z^2)])))-(q/(4(Subscript[πϵ, 0])(Sqrt[(((x+a))^2)+(((y-a))^2)+(z^2)])))-(q/(4(Subscript[πϵ, 0])(Sqrt[(((x-a))^2)+(((y+a))^2)+(z^2)])))           (1)	
で与えられる。

定義

電位phiを定義しておく。

Clear[phi];
phi[x_,y_,z_]:=q/(4 Pi ep0 Sqrt[(x-a)^2+(y-a)^2+z^2]) +q/(4 Pi ep0 Sqrt[(x+a)^2+(y+a)^2+z^2]) -q/(4Pi ep0  Sqrt[(x+a)^2+(y-a)^2+z^2])-
q/(4Pi ep0  Sqrt[(x-a)^2+(y+a)^2+z^2])
ep0=1;q=1;a=1;

計算がちゃんとできるか確かめておこう。

N[phi[0.5,0.5,0]]

電位の表示

ほとんどのMathematicaの関数は直交座標系を想定しているので、この問題のように直交座標系になっているものは簡単に表示できる。PlotPointsを省いて表示させるとPlotPointsの役割がわかるだろう。

graph1=ContourPlot[phi[x,y,0],{x,0,3},{y,0,3},ContourShading->False,PlotPoints->30]

Plot3D[phi[x,y,0],{x,0,3},{y,0,3}]

電位は電荷がある(a,0)のところで高くなっており金属表面では一定値になっていることが分かる。

電界の計算

E=-grad ϕなので、ExとEyを計算する。このときには遅延割り当て(:=)ではなくて即時割り当て(=)を用いる。

Clear[Ex,Ey,L];
Ex[x_,y_,z_]=-D[phi[x,y,z],x]
Ey[x_,y_,z_]=-D[phi[x,y,z],y]
L[x_,y_,z_]:=Sqrt[Ex[x,y,z]^2+Ey[x,y,z]^2];

電界の表示

ベクトルの表示にはパッケージを読み込む必要がある。詳しくはヘルプブラウザーを使って調べること。

Needs["VectorFieldPlots`"]

ベクトルの長さLを使って長さを規格化する。Lを使わないとどうなるか試してみるとよい。

graph2=VectorFieldPlot[{((Ex[x,y,0])/(L[x,y,0])),((Ey[x,y,0])/(L[x,y,0]))},{x,0,3},{y,0,3}]

電位の等高線図と重ねて表示する.
電界ベクトルの向きは必ず等電位面と垂直になっていることが分かる。
また導体表面(x=0,y=0)のところでは電界ベクトルは垂直になっている。

Show[graph1,graph2]

次のようにScaleFunctionを利用すれば、長さを規格化しなくても、一定値の0.5をとらせることができる。

VectorFieldPlot[{Ex[x,y,0],Ey[x,y,0]},{x,0,3},{y,0,3},ScaleFunction->(0.5`&)]

ちなみに、４つの電荷がどのような電界や電位を作っているかということを見ると次の様になる。これらの表示の第一象限のみを考えて解として用いている。

graph1=ContourPlot[phi[x,y,0],{x,-3,3},{y,-3,3},ContourShading->False,PlotPoints->30]

Plot3D[phi[x,y,0],{x,-3,3},{y,-3,3}]

graph2=VectorFieldPlot[{((Ex[x,y,0])/(L[x,y,0])),((Ey[x,y,0])/(L[x,y,0]))},{x,-3,3},{y,-3,3}]

Show[graph1,graph2]