クーロンポテンシャルを場合分けによって書くことができたが,実際に計算する場合には,いちいち分けるのは面倒だ。といっても高々六類型なのでどうってことはないといえばない。
軸対称電荷分布(2)でやったように,ヘヴィサイドの階段関数を使えば,r<, r> が表現できるので,クーロンポテンシャルのルジャンドル成分 ˜Vn(r) を一つの関数を使って表わせる。そこで,R(t)=1√1+t2, a=1.0, c=0.707107の場合について具体的に計算してみた。
rs[s_, r_] := s* HeavisideTheta[r - s] + r*HeavisideTheta[s - r];
rl[s_, r_] := r* HeavisideTheta[r - s] + s*HeavisideTheta[s - r];
R[t_] := 1/Sqrt[1 + t^2];
v[n_, r_] := 3*NIntegrate[ 1 + 0 * LegendreP[n, t]
*rs[s, r]^n * s^2/rl[s, r]^(n + 1), {t, 0, 1}, {s, 0, R[t]}]
v0 = v[0, 0]
2.64412
Plot[{v0, v[0, r]}, {r, 0, 3}]
Plot[{v0, v[2, r]}, {r, 0, 3}, PlotRange -> {2.55, 2.65}]
Plot[{v0, v[4, r]}, {r, 0, 3}, PlotRange -> {2.64, 2.665}]
Plot[{v0, v[6, r]}, {r, 0, 3}, PlotRange -> {2.638, 2.648}]
ここで,被積分関数の中に 1 を加えたのは,積分値が 0 になって計算が不安定になるのを防ぐためである。その部分の定数 v0 を取り除くために,ルジャンドル関数の前に 0 をかけて積分することで v0 を求め,実際の計算ではこれを取り除く。
その結果は次のようになった。˜V4(r), ˜V6(r)が,r<c でゼロになるのは,∫ds からでてくる R(t)n−2 がルジャンドル関数のn次未満の線形結合になるためである。
図:クーロンポテンシャルのルジャンドル成分 ˜Vn(r) (n=0,2,4,6)
0 件のコメント:
コメントを投稿