今度は別の問題を考えてみましょう。 やはり変数分離できる場合を考えます。

この問題では、円柱座標(r, θ, z)での Laplace方程式

の解を、以下の境界条件で解く事になります。

変数分離解 φ(r, θ, z) = R(r) Θ(θ) Z(z) の形で探します。 ただし、θ方向に関して電位分布は一定なので、Θ(θ)=1 (Constant) と置くことができるので、Laplace方程式は次の2つの常微分方程式

になります。ここでα、βは α²+β²=0を満たす任意の定数です。 1つめの常微分方程式をMathematica で DSolveで解いてみると、

となります。ここで、BesselJ や BesselY は Bessel関数と呼ばれる 特殊関数です。 これについて少し解説しておきます。

Bessel 関数は Bessel方程式と呼ばれる2階の常微分方程式

の解として定義され、その二つの独立解として J_n(r)、Y_n(r) と書かれます。一般解は R(r) = C₁ J_n(r) + C₂ Y_n(r)です。 J_n をn 次の第一種 Bessel 関数、 Y_n をn 次の第二種 Bessel 関数(または n 次の Neumann 関数) と呼びます。 Y_n(r) は時にN_n(r)とも書かれます。

具体的な関数形は省略しますが、上の微分方程式が確定特異点を 持つことに 注意してその周りで級数解法で解けば、2つの独立解 J_n、J_-n の形を求めることができます。ただし nが整数の場合は J_n、J_-n の間に関係式が生じてしまうため独立ではなくなってしまいます。 そこで J_nとJ_-n の線型結合より Neumann関数 Y_n を定義します。また以下の常微分方程式

を Modified Bessel 方程式と呼び、その独立解 I_n(r)、K_n(r) をそれぞれ n 次の第一種、第二種 Modified Bessel 関数と呼びます。 このとき Modified Bessel 方程式の一般解は R(r)=C₁I_n(r)+C₂K_n(r) となります。 また、I_n、K_nは上の J_n や Y_nとは独立ではなく Modified Bessel 方程式の解は R(r)=C'₁J_n(I r)+C'₂K_n(I r) とも書けます。

これらの J_n、Y_n、I_n、K_n など（他にもその派生関数があります。）を総称して Bessel 関数と呼び、Mathematica ではそれぞれ BesselJ[n, x]、 BesselY[n, x]、 BesselI[n, x]、 BesselK[n, x] と書きます。どんな関数形をしているか、最初の数項（n=0,1,2,3） をプロットしてみると、それぞれ

となります。それぞれの x → 0 の振舞いなど、基本的な特性を よく覚えておいてください。

さて、円柱座標(r, θ, z)での Laplace方程式は Bessel関数を使えば

となります。α、βについての重ね合わせは形式的に 和の記号を用いましたが一般的には積分になることもあります。

あとは先ほどと同様、任意係数A(α),B(α),C(β),D(β) を境界条件から決定してやるだけです。

まず円筒の内部では電位は至るところ正則だと考えられます。 そこで、円筒内部のz軸に相当する r=0,-L <= z <= L の範囲の電位を見てみます。すると、 J₀(0) は上の グラフのように有限な値ですが、Y₀(0) は発散して います。したがって未定係数B(α)=0 である必要があ ります。

次に第二、三の境界条件 φ(ｒ, -L)=φ(ｒ, L)=φ₀, (0<r<a) を考えると、C(β)-D(β)=0 が分かります。

さらに第一の境界条件φ(a,z)=0, (-L<=z<=L)より

とならねばなりません。これは J₀(- I α a) = 0 とならねばならない事を言っています。 実は Bessel 関数 J₀(z), (z 複素数) のゼロ点は実軸上 にしかなく、飛び飛びに孤立しています。Mathematica で確かめて みます。ContourPlot を使って BesselJ[0, x + I y] の実部と虚部がともにゼロになるところを見てみます。

と打つと、以下の図が得られます。

青が実部がゼロとなる所で、赤が虚部がゼロとなるところです。 赤と青が重なっているところは確かに実軸上に飛び飛びにしかない ようです。実軸上でのBesselJ[0,x] の振舞いは上で見た ように振動しています。この x 軸を横切っている点の値 (J₀(x)の根)を それぞれ、 ... p_-3,p_-2,p_-1, p₁,p₂,p₃,..., p_n, と値の小さい方から名前をつけます(n 番目の根 = p_n)。

そうすると、第一番目の境界条件から得られた J₀(- I α a) = 0 から次の式が得られます。

(n=...-3,-2,-1,1,2,3,...、ゼロを除く)であるこ とが分かります。 またα²+β²=0 でしたので β_n= + p_n/a、または β_n= - p_n/aとなります。 またしても α、βは連続的な値を取れないようです。 これらをφ(x,y)に代入し、指数関数を双曲余弦関数に置 き換え、n についての和で負の整数の部分は正の整数 の和と同じである事などを用いて比例係数などを吸収するよう に未定係数を再定義すると、 (A(α_n)C(β_n)=a_n 等の置き換えを行った。)

となります。

最後に残る条件は φ(r,-L)=φ(r, L) = φ₀, (0 <= r < a) です。これは、 以下の境界条件の式を満たす様にa_nを決めることです。

この式からa_n は以下のようにして求めることができます (Bessel 関数の直交性を用いる)。

まず、両辺に r J₀(p_kr/a) をかけて r で 0 から a まで積分します。

---(*)

となります。右辺の積分は k と n が等しい場合と異なる場合で評価します。

まず、p と k が異なる時を見てみましょう。 p_n、p_k をそれぞれpn、pkと 表して Mathematica に計算させると、

となります。pk、pnの定義より (それぞれ Bessel J₀ 関数の相異なるゼロ点)、 J₀(p_k)=0、 J₀(p_n)=0 であり、pk とpn が異なることからこれは ゼロです。

次に k=n の時は、あらわに pn=pk とおいて Mathematica で計算すると、

となり、p_nの定義より

となります。 結局、式(*) の右辺の n についての和は n = k の所だけが残ります。

また左辺の積分も同様に Mathematica で実行すれば、

となります。従って a_n は(上の k を n と改めておき直して)、

となることが分かります。 ここではa_n をあらわに p_n の関数としてan[pn] と定義しました。

これで全ての未定係数は求まりました。

どんな解が求まったのか、またお絵描きしましょう。

それにはまず、p_n を具体的に求めなければなりません。それには先週やった FindRootが使えます。初期値としては原始的ですが、 グラフから読み取った大まかな数字を手で与えてやることにします。 グラフ

より、{x, 0, 50} の範囲にあるｘ軸との交点の座標をざっと読み取り、 それを Newton 法の初期値とします。すなわちグラフから読み取った ｘの大まかな値をリスト形式で与えておいて、それを使ってFindRoot します。

guess1[[i]] とするとリスト guess1 の i 番目の要素を参照できる ことを利用しています。 これは前出の Part［guess1, i] の短縮形に相当します。 しかし FindRoot を16回も行うのは大変なので、関数 Table[ ] を使用して楽をしています。 Table[hogehoge, {i, 1, 16} ] は、1 から 16まで i を1ずつ変えながら hogehoge を評価し、 その値のリストを出力する関数 です。これは前出の Do[ ] に似ていますが、結果をリスト形式 で出力してくれるところが相違点です。

このリストから例えば 1 番目の x の値を引き出すには、list1[[1]] でよさそうですが、やってみると

となって余計なものがたくさんついています。 この出力はよく見ると要素が1個のリストです。つまり List1 はリストのリスト、すなわち 16ｘ1 の行列なのです。 したがって、1行1列目を出力させる、すなわち

とすれば、少し余計なものを減らすことができます。 ここからさらに→以降だけを取り出すには、 前出の関数 First[ ] に類似した（ある意味で逆の）、 Last[ ]という関数 を使います。まとめると p_nは

になります。関数の定義でよく使う:= が登場しましたが、 := は正確には SetDelayed と言って、右辺が評価されるまで =（通常のSet）を保留するためのもの です。今の段階では n が特定の数に決まっていないので:= を使用しています。

あるいは同じことが

としても可能です。list1 の n行1列目の要素は ｘ→hogehoge ですが、 →の前が list1[[n,1,1]] で、→の後が list1[[n,1,2]] です。 シンプルなので、慣れるとこちらのほうがわかりやすくなると思います。

これで p_n の具体的な値がわかったので、プロットに移りましょう。 pn として 16 個しか用意できなかったので、無限項の和は16項目までの和 で近似します。

こんな図が描けたでしょうか？

今回も両端にギザギザが現れてしまいました。 これはやはり無限項の和を有限項の和で近似したことより発生しています。 これをよりスムーズにするには pn としてもっとたくさんの値を用意しな ければなりません。{x, 0, 50} のかわりに {x, 50, 100} にして上と同じことを繰り返します。 グラフの交点を読み取り、それをもとに Newton法で解を求めると

となります。全体のリストは

とすればよさそうですが、結果をよく見ると 2つのリストが並んで いるだけで、1つのリストとして融合してはいません。 こんなときは リストをフラットにする(融合する)関数 Flatten[ ] を用います。1段階だけフラットにすればよいので、

で、望ましい出力が得られます。これを使って

とすれば、

となります。

実は、Bessel関数 J_n(x) のゼロ点を与えてくれる コマンド BesselZeros があります。BesselJZeros[n, num] とすると J_n(x) のゼロ点を原点に近い順番にnum個 計算してリストにして返してくれます。あるいは、N[BesselJZero[n, num]] （最後に s がないことに注意）とすると、J_n(x) のゼロ点を原点に近い方から num 番目のゼロ点の近似値を返してくれます。 これらを使えば上で行ったFindRoot の部分をもっと簡単に出来ます。皆さんやってみて下さい。