index.html

電界中の誘電体の電束線

一定電界中に誘電率εで半径aの球状誘電体を入れたときの電気力線描くことを考えよう。電界は次の式のようになっている。
　　　　　=(1+ $(\[Epsilon] - \[Epsilon]_0)/(\[Epsilon] + 2\[Epsilon]_0)$ ) cosθ    r>a
　　　　　= $(3\[Epsilon]_0)/(\[Epsilon] + 2\[Epsilon]_0)$ cosθ    r<a
           =-(1- $(\[Epsilon] - \[Epsilon]_0)/(\[Epsilon] + 2\[Epsilon]_0)$ )sinθ      r>a
           = $(-3\[Epsilon]_0)/(\[Epsilon] + 2\[Epsilon]_0)$ sinθ    r<a
電束密度は
           D=e0 E r>a
           D=e E r<a
           である。

プログラム

定義

定数の定義

E0 = 1 ; (* 磁界の強さ *)
a = 0.87 ; (* 球の半径 *)
dr = 0.02 ;
ep = 3 ; (* 誘電体の誘電率 *)
ep0 = 1 ; (* 真空の誘電率 ep/ep0 = 3とする*)

電界の定義。ここではIf文を使って定義することにする。

Clear[Dr, Dth] ; Dr[r_, th_] := If[r>a, ep0 * (1 + (ep - ep0)/(ep + 2ep0) 2a^3/r^3) E0 Cos[th], ep * 3ep0/(ep + 2ep0) E0 Cos[th]] ;

Dth[r_, th_] := If[r>a, -ep0 * (1 - (ep - ep0)/(ep + 2ep0) 2a^3/r^3) E0 Sin[th], -3 * ep * ep0/(ep + 2ep0) E0 Sin[th]] ;

球座標でφ=0のとき、直交座標に変換する。

r[x_, y_] := Sqrt[x^2 + y^2] ; th[x_, y_] := Which[x>0, ArcTan[y/x], x<0, Pi + ArcTan[ ... 0] ; Dx[x_, y_] := Dr[r[x, y], th[x, y]] Cos[th[x, y]] - Dth[r[x, y], th[x, y]] Sin[th[x, y]] ;

作画のためのプログラム

(x,y)にいるときに、つぎに電気力線がどちらの方向に進むか計算する。それは
,
で与えられる。

$(* x軸方向にどのくらい移動するか *) calx[x_, y_] := Module[{length, ret = 0}, len ... ; ret = y + Dy[x, y]/r * dr ; Return[ret] ; ] ;$

計算が終了したかどうかを判定する関数。この場合はx<0の領域のみ計算するので、下のようなプログラムになる。たとえば、r=を計算して、r<aの時に終了するとか変更することができる。

$limit[x_, y_] := Module[{r}, Return[x>0] ;] ;$

実際に作画するところ。liというリストにLine[{{x1,y1},{x2,y2}, ---,{xn,yn}}]という線を描くグラフィック要素を付け加えていく。はじめは(-2.66,-2)という点から始まって次の点を計算し、limit[x,y]で無い限り点を付け加える。左右対称を利用している。
ysの範囲を変更して、プログラムを実行し、li, tsl, tsrを?で表示させてみるとよい。

cal[] := Module[{x, y},
li = {} ;
For[ys = 2, ys >= -2, ys -= 0.25 ... evel[0.5], Disk[{0, 0}, a]}, li}, AspectRatio -> 3/4, PlotRange -> {{-2.66, 2.66}, {-2, 2}}]] ;
]