電界中の誘電体の電束線

一定電界(Subscript[E, 0])

      中に誘電率ϵで半径aの球状誘電体を入れたときの電気力線描くことを考えよう。電界は次の式のようになっている。
　　　　　(Subscript[E, r]) = (1+((ϵ - (Subscript[ϵ, 0]))/(ϵ + 
                2  (Subscript[ϵ, 0]))) ((2  (a^3))/(r^3))) (Subscript[E, 0])cosθ    r>a
　　　　　(Subscript[E, r])=\(TraditionalForm\`\(3  ϵ\_0\)\/\(ϵ + 2  ϵ\_0\)\) (Subscript[E, 0])cosθ    r<a
	       (Subscript[E, θ])=-(1-((ϵ - (Subscript[ϵ, 0]))/(ϵ + 
                2  (Subscript[ϵ, 0]))) ((2  (a^3))/(r^3)))(Subscript[E, 0])sinθ      r>a
	       (Subscript[E, θ])=\(TraditionalForm\`\(\(-3\) ϵ\_0\)\/\(ϵ + 2  ϵ\_0\)\) (Subscript[E, 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;

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

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

\(TraditionalForm\`球座標でφ = 0\)のとき、直交座標に変換する。

r[x_, y_] := Sqrt[x^2 + y^2]; 
    th[x_, y_] := 
      Which[x > 0, ArcTan[y/x], x < 0, Pi + ArcTan[y/x], x == 0, 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]];
Dy[x_, y_] := 
        Dr[r[x, y], th[x, y]] Sin[th[x, y]] + 
          Dth[r[x, y], th[x, y]] Cos[th[x, y]];

作画のためのプログラム

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

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

    \(caly[x_, y_] := 
        Module[{r, rdx, rdy, ret = 0}, 
	
          r = Sqrt[
              Dx[x, y]^2 + 
                Dy[x, y]^2]; 
\ \ \ \ \ \ \ \ \ ret = 
            y + Dy[x, y]/r*dr; 
	
          Return[ret];
];\)\)\)\)

計算が終了したかどうかを判定する関数。この場合はx<0の領域のみ計算するので、下のようなプログラムになる。たとえば、r=\(TraditionalForm\`\(\@\(x\^2 + 
              y\^2\)\) を計算して、\)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, 
	
            x = \(-2.66\); 
	y = ys; 
	
            tsl = {{x, 
                  y}}; \  (*\ 左からの計算*) 
	
            tsr = {{\(-x\), 
                  y}};  (*\ 右からの計算*) 
	
            While[Not[limit[x, y]], 
	
              x = calx[x, y]; 
	
              y = caly[x, y]; 
	
              tsl = Append[tsl, {x, y}]; 
	
              tsr = Append[
                  tsr, {\(-x\), 
                    y}];
\ \ \ \ \ \ ]; 
	li = Append[li, {Line[tsl], 
                  Line[tsr]}];
\ ]
 (*\ 計算が終わったので作画する\ *) \ 
\ \ \ Show[
            Graphics[{{GrayLevel[0.5], Disk[{0, 0}, a]}, li}, 
              AspectRatio -> 3/4, 
              PlotRange -> {{\(-2.66\), 2.66}, {\(-2\), 
                    2}}]];
]\)\)\)

実行

epとep0の値を変化させて計算させると面白い。ep/ep0<1の時にははじくようになる。

ep = 3; ep0 = 1; cal[]

ep = 1; ep0 = 3; cal[]

ep = 0; ep0 = 3; cal[]

PlotVectorFieldを使うと簡単に表示できるが、美しさは今一歩。

Needs["Graphics`PlotField`"] ;

ep = 3; ep0 = 1;
PlotVectorField[{Dx[x, y], Dy[x, y]}, {x, -2, 2}, {y, -2, 2}, 
      ScaleFunction -> (0.5 &)]