ラベル Mathematics の投稿を表示しています。 すべての投稿を表示
ラベル Mathematics の投稿を表示しています。 すべての投稿を表示

2016年8月30日火曜日

微分方程式を行列方程式に


まず正規直交関数\( \lbrace \phi_{k} \rbrace^{\infty }_{k=1} \)を考える。直交多項式は\( k=0 \)から始まるのでほんのすこし面倒だが、ここでは正規直交関数の列を便宜的に\( \phi_{k} \; , k=1,2,3 \ldots \)とする。 例えば、正規直交型のルジャンドル多項式は \begin{align} \phi_{k}(z) = \sqrt{ \frac{2k -1}{2} } P_{k-1} (z) \end{align} などとして取り扱えば良い。このような正規直交型の関数を使って関数を展開する。例えば、基底の数をnとすると \begin{align} F(z) = \sum^{n}_{k=1} c_{k} \phi_{k} (z) \end{align} といった感じだ。

ここでは正規直交化された一般ラゲール関数\( \phi^{(\alpha)}_{n}(z) = \sqrt{ \frac{n!}{(n + \alpha)!} } z^{ \frac{\alpha }{2} } e^{ - \frac{z}{2} } L^{(\alpha ) }_{n } (z) \) を例にとって、その1階、2階微分を計算すると \begin{align} \frac{\partial \phi^{(\alpha ) }_{n } (z) }{\partial z} & = \sqrt{ \frac{n!}{(n + \alpha)!} } \left[ \left( \frac{\alpha }{2z} - \frac{1 }{2 } \right) e^{ \frac{\alpha }{2} \log{z} - \frac{z}{2} } L^{(\alpha ) }_{n} ( z ) + e^{ \frac{\alpha }{2} \log{z} - \frac{z}{2} } (-1 ) L^{(\alpha +1 ) }_{n-1} (z) \right] \\ \frac{\partial^{2} \phi^{(\alpha ) }_{n } (z) }{\partial z^{2} } & = \sqrt{ \frac{n!}{(n + \alpha)!} } \left[ \left\lbrace - \frac{\alpha }{2z^{2}} + \left( \frac{\alpha }{2z} - \frac{1}{2} \right)^{2} \right\rbrace e^{ \frac{\alpha}{2} \log{z} - \frac{z}{2} } L^{(\alpha ) }_{n} (z) + \left( - \frac{\alpha }{z } + 1 \right) e^{ \frac{\alpha}{2} \log{z} - \frac{z}{2} } L^{(\alpha +1) }_{n-1} (z) + e^{ \frac{\alpha}{2} \log{z} - \frac{z}{2} } L^{ ( \alpha + 2 ) }_{ n - 2 } (z) \right] \end{align} ここで2階の微分方程式のための演算子を\( \hat{L} \)とでも書くと、その斉次微分方程式は \begin{align} \hat{L} F = E F \end{align} というようになる。 Fを展開して \begin{align} \hat{L} F = A_{2} \frac{\partial^{2} F }{\partial z^{2} } + A_{1} \frac{\partial F }{\partial z } + A_{0} F = \sum_{k} \hat{L} c_{k} \phi_{k} = E \sum_{k} c_{k} \phi_{k} \end{align} といったようになる。この演算子の行列要素を、 \begin{align} L_{jk} = \int \phi^{*}_{j} \hat{L } \phi_{k} dz \end{align} といった感じで定義する。 よって行列方程式のj行目は \begin{align} \sum^{n}_{k=1} L_{j k } c_{k} = E c_{j} \; , \quad (j=1,2, \ldots , n) \end{align} となって連立方程式を解くことになる。基底の数 n と連立方程式の数は同じである。

2016年8月27日土曜日

接続法は実用的ではなさそう

シュレディンガー方程式を数値的に解いて固有値を得ようとするときに、前方・後方からRK4などを使って解いてとある接続点で接続させるという手法は実用的ではないという印象です。 二階の微分方程式のうち、他方の余分な解数値解が入ってきて、解が接続点に依存しがちなことが問題だなあという感じです。

時間非依存シュレディンガー方程式を行列形式に書き換えることもでき、この方が実用的かもしれない。
ある状態関数\( | \psi \rangle \)がヒルベルト空間に属しているならば、その空間の正規直交基底で展開できるので、\( | \psi \rangle = \sum_{k} c_{k} | k \rangle \)となる。時間非依存のシュレディンガー方程式は、 \begin{align} H | \psi \rangle = E | \psi\rangle \; , \quad \rightarrow H \sum_{k} c_{k} | k \rangle = E \sum_{k} c_{k } | k \rangle \end{align} となる。 ここで \( \hat{1} = \sum_{m} | m \rangle \langle m | \)であるから \begin{align} \hat{1} H \sum_{k} c_{k} | k \rangle = E \sum_{k} c_{k } | k \rangle \; , \quad \rightarrow \sum_{m, k} \langle m | H | k \rangle c_{k} | m \rangle = E \sum_{k} c_{k } | k \rangle \end{align} となる。ここで \( \langle m | H | k \rangle = H_{mk} \)とすると \begin{align} \sum_{m, k} H_{mk} c_{k} | m \rangle = E \sum_{k} c_{k } | k \rangle \end{align} となり、左から \( | j \rangle \)を作用させると、\( \langle j | m \rangle = \delta_{jm} \)であるので、 \begin{align} \sum_{k} H_{jk} c_{k} = E c_{j} \; , \quad \rightarrow {\bf H} {\bf c} = E {\bf c} \end{align} を得る。ここでは\( {\bf H} \)は行列、\( {\bf c} \)は列ベクトルである。
これは行列の固有値問題になります。あとはこれを解くだけで固有値が求まります。

関数を直交関数で展開する時に、指数関数の近似があまり良くないことはしばしばなのではと思います。\( e^{ -\frac{z}{2} } \left( z^{2} - 2 z +2 \right) \)を正規直交型の一般ラゲール関数 \begin{align} \sqrt{ \frac{n!}{ (\alpha + n) ! } } z^{ \frac{\alpha }{2} } e^{ - \frac{z }{2} } L^{(\alpha )}_{n} (z) \end{align} で展開してやるとグラフのようになります。ここでは\( \alpha=1 \)で計算しています。左が自然スケールで右が対数スケールです。右のグラフでは関数の絶対値をとっています。グラフにみられる#BFsというのは基底関数の数です。



基底関数の数が上がれば近時の精度は上がりますが、zの値が大きい領域ではむしろ基底関数の数が少ないほうが精度がよいというような感じになっています。これは当然の話しで、一般ラゲール多項式\( L^{(\alpha ) }_{n} (z) \)のもっとも優位な項は \( z^{n} \)ですから、正規直交化させたラゲール関数のもっとも優位な項は \begin{align} \sqrt{ \frac{ 1 }{ n! (\alpha + n) ! } } (-1)^{n} z^{\frac{\alpha }{2} } e^{- \frac{z}{2} } z^{n} \end{align} となります。基底関数の数が上がるほどnの値が大きくなるので、\( e^{- \frac{z}{2} } \)との差は大きくなります。元の関数は\( e^{ -\frac{z}{2} } \left( z^{2} - 2 z +2 \right) \)だったので、zがかなり大きい領域では基底の数が上がるほど、近似関数と元の関数との差が大きくなります。

2016年8月24日水曜日

ウィットテーカー関数


合流型超幾何関数(confluent hypergeometric function)は合流型超幾何微分方程式の解である。これは2階の微分方程式であって、 \begin{align} z \frac{ d^{2} \omega }{ d z^{2 } } + ( b - z ) \frac{d \omega }{d z } - a \omega = 0 \; \end{align} というような型である。
この方程式を変形させるために \begin{align} \omega = z^{ \frac{-b}{2 } } e^{ \frac{z}{2 } } \chi \; \\ a = \frac{1}{2} - \kappa + \mu \; , \quad b = 1 + 2 \mu \; \end{align} としてやる。これを合流型超幾何微分方程式に代入してやると、 \begin{align} \chi'' + \left( \frac{- \mu^{2 } + \frac{1}{4} }{z^{2 } } + \frac{\kappa }{z } - \frac{1}{4} \right) \chi = 0 \; \end{align} となる。これをウィットテーカーの標準形あるいはウィットテーカーの方程式と呼ぶ。また、ここではプライムの記号は \( z \) での微分を意味する。 \( \chi' = \frac{\partial \chi}{\partial z} \).
この方程式の解はウィットテーカー関数と呼ばれ、以下のような式で表される。 \begin{align} M_{\kappa, \mu} & = z^{\mu + \frac{1}{2} } e^{ - \frac{z}{2} } M(\mu + \frac{1}{2} - \kappa, 2 \mu + 1, z ) = z^{ \frac{b}{2} } e^{ - \frac{z}{2} } M(a,b; z ) \; , \\ W_{\kappa, \mu} & = z^{\mu + \frac{1}{2} } e^{ - \frac{z}{2} } U(\mu + \frac{1}{2} - \kappa, 2 \mu + 1, z ) = z^{ \frac{b}{2} } e^{ - \frac{z}{2} } U(a,b; z ) \; , \end{align} ここでは \( {\rm Arg}(z)\in (- \pi, \pi ] \)、でありパラメータは \( \kappa = \frac{b}{2} - a \; \; , \mu = \frac{ b - 1 }{2} \)と表される。
例として\( \kappa=1.0, \; \mu=3 \) と \( \kappa=2.0, \; \mu=4 \) の場合のウィットテーカー方程式を解き、その中で\( M_{\kappa, \mu} (z) \)関数をグラフに示す。ここではヌメロフ法を使う。


解析値と数値解は近い値になっていることがわかる。 また、ウィットテーカー関数の\( M_{\kappa, \mu} (z) \)関数の一階導関数は \begin{align} \frac{\partial }{\partial z} M_{\kappa , \mu} (z) = \left( -\frac{1}{2} + \frac{\mu + \frac{1}{2} }{z} \right) M_{\kappa, \mu} (z) + e^{ - \frac{z}{2} } z^{\mu + \frac{1}{2 } } \frac{\mu + \frac{1}{2} - \kappa }{1 + 2 \mu} M ( \mu + \frac{1}{2} - \kappa + 1, 2 \mu + 1 + 1; z ) \; , \end{align} のようになる。ここではクンマーの合流型超幾何関数の微分の性質 \begin{align} \frac{\partial }{\partial z} M(a, b ; z ) = \frac{a}{b} M(a+1, b+1; z) \; \end{align} を使っている。

ここで、アダムス・バッシュフォース・モールトン法を使って\( \kappa=2.0 , \; \mu = 4.5 \) と \( \kappa=3.0, \; \mu=5.0 \)の場合のウィットテーカー関数を解いてやり、 \( M_{\kappa , \mu} \) 関数の数値解を求める。


するとグラフのようになり、ウィットテーカー関数の解析解と数値解が近い値になっていることがわかる。

2016年8月21日日曜日

ヌメロフ法


\( y''(z) = P(z) = g(z) y(z) \)というようなタイプの微分方程式の数値解を得ようとする時には、ヌメロフ法が便利だと思う。


ヌメロフ法の簡単な説明としては、まずテイラー展開で \begin{align} y_{n+1} + y_{n-1} - 2 y_{n} = \frac{2}{2!} y''(z) h^{2} + \frac{2}{4!} y''''(z) h^{4} + \mathcal{O}(h^{6}) \; , \end{align} とする。 \( y_{n+1} = y(z+h), \; y_{n-1} = y(z-h) \)
そして\( y''''(z) \)については、 \( y''(z) = P(z) \)であることを使って、 \begin{align} P_{n+1} + P_{n-1} - 2 P_{n} = \frac{2}{2!} P''(z) h^{2} + \mathcal{O}(h^{4} ) \; , \\ h^{2} y''''(z) = h^{2} P''(z) = P_{n+1} + P_{n-1} - 2P_{n} - \mathcal{O}(h^{4 } ) \; , \end{align} となることがわかる。それを利用して、 \begin{align} y_{n+1} + y_{n-1} - 2 y_{n} & = P_{n} h^{2} + \frac{ 2 h^{2} }{4! } \left( P_{n+1} + P_{n-1} - 2P_{n} - \mathcal{O}(h^{4 } ) \right) + \mathcal{O}(h^{6 } ) \; , \\ & = \frac{h^{2} }{12 } P_{n+1} + \frac{5h^{2}}{6} P_{n} + \frac{h^{2}}{12 } P_{n-1} + \mathcal{O}(h^{6} ) \; \end{align} となることがわかる。よって、 \begin{align} \left( 1 - \frac{h^{2} }{12} g_{n+1} \right) y_{n+1} + \left( 1 - \frac{h^{2} }{12} g_{n-1} \right) y_{n-1} = \left( 2 - \frac{5 h^{2 } }{6 } g_{n} \right) y_{n} + \mathcal{O}(h^{6} ) \; \end{align} を得る。これがヌメロフ法である。

ヌメロフ法を使ってこの手の微分方程式を解いてみる。 \begin{align} \frac{d^{2} \chi_{n} }{d z^{2}} + \left( \frac{ - \frac{\alpha^{2} }{4} + \frac{1}{4} }{z^{2} } + \frac{ n + \frac{\alpha}{2} + \frac{1 }{2} }{z} - \frac{1}{4} \right) \chi_{n} = 0 \; \end{align} といった方程式はどうだろう。この解析解は \begin{align} \chi_{n } (z) = z^{ \frac{\alpha + 1}{2} } e^{ \frac{-z}{2} } L^{(\alpha ) }_{n } (z) \end{align} である。
例として、 \( \alpha=1.0 \; , \; N=2 \)と \( \alpha=2.5 \; , N=3 \), \( \alpha=4.5, \; N=5 \), \( \alpha=5, \; N=7 \) の場合では、


\( \alpha=1.0, \; N=10 \) と \( \alpha=1.0, \; N=20 \)の場合は



のようになり、Nが大きくなるほど収束は遅くなることがわかる。
全体的に解析解と数値解はかなり近い値となっている。計算速度もかなり速かったと思う。

2016年8月19日金曜日

ロドリゲス公式よりも便利な公式

有名な多項式を計算する場合にロドリゲスの公式をつかって計算する人がいるはず。ロドリゲス公式よりもはるかに速く計算できる方法というと漸化式だろうと思う。

漸化式を使わずに、一般ラゲール多項式を計算する公式はある。

\begin{align} L^{ (\alpha) }_{n } (z) = \sum^{n}_{k=0} \frac{ (-1 )^{k} }{ k! (n-k)! } z^{k } \frac{ (\alpha + n)! }{ (\alpha + k)! } \quad \end{align} 便宜上この公式をformula GLとしよう。この公式はロドリゲス公式からライプニッツの公式を使って導出できる。
この公式をつかって \begin{align} L^{(\alpha ) }_{0} (z) =1 \; , \quad L^{(\alpha )}_{1} (z) = \alpha + 1 - z \; , \quad L^{(\alpha ) }_{2} (z) = \frac{(\alpha + 2) (\alpha + 1 ) }{2 } - (\alpha + 2) z + \frac{ z^{2 } }{2} \; , \\ L^{(\alpha ) }_{3} (z) = \frac{1}{6} (\alpha + 3 ) (\alpha + 2 ) (\alpha + 1) + \frac{(-1)}{2} (\alpha +3)(\alpha + 2) z + \frac{( \alpha + 3) }{2} z^{2} + \frac{ (-1) }{6} z^{3} \end{align} を計算してみた。 また、4次5次の多項式も \begin{align} L^{(\alpha ) }_{4} (z) & = \frac{1}{24} (\alpha + 4) (\alpha + 3) (\alpha + 2) (\alpha + 1) + \frac{(-1)}{6} (\alpha + 4)(\alpha + 3)(\alpha +2) z + \frac{1}{4} (\alpha + 4)(\alpha + 3)z^{2} + \frac{(-1)}{6} (\alpha +4) z^{3} + \frac{1}{24} z^{4 } \\ L^{(\alpha ) }_{5} (z) & = \frac{1}{120}(\alpha + 5)(\alpha + 4) (\alpha +3) (\alpha +2) (\alpha + 1) + \frac{(-1)}{24} (\alpha +5)(\alpha +4) (\alpha +3) (\alpha +2) z + \frac{1}{12} (\alpha +5)(\alpha +4) (\alpha +3) z^{2} \nonumber \\ & \quad + \frac{(-1) }{12} (\alpha +5) (\alpha +4) z^{3} + \frac{1}{24} (\alpha +5) z^{4} + \frac{(-1) }{120} z^{5 } \end{align} となる。
この公式の利点は漸化式にくらべて各項の係数を把握しやすいというところだと思う。もちろん低次の項に限っての話でしょうが。 こららの他にクンマー級数(クンマーの超幾何関数)をもちいて \begin{align} L^{(\alpha ) }_{n} (z) = \frac{ (\alpha + n) !}{ n! \alpha ! } M(-n, \alpha+1; z ) \end{align} とかくこともできる。この公式をformula Kとすると、もちろんformula GLで計算される結果と一致する。 たとえば\( \alpha=1.5, \; N=4 \) と \( \alpha=5.0, \; N=7 \)の場合を計算しグラフに示すと
といったようになる。

またformula GLを使うことで \begin{align} \frac{\partial L^{(\alpha) }_{n } (z) }{\partial z} = - L^{(\alpha + 1) }_{n-1} (z) \end{align} という公式も証明できる。

単純に、\( L^{(\alpha ) }_{n} (z ) \)の微分を考えてやると、 \begin{align} \frac{\partial L^{(\alpha) }_{n } (z) }{\partial z} & = \sum^{n}_{k=1} \frac{ (-1)^{k} }{(k-1)! (n-k)!} z^{k-1} \frac{(\alpha + n)!}{(\alpha + k)! } \\ & = - \sum^{n}_{k=1} \frac{ (-1)^{k-1} }{(k-1)! ( (n -1) - (k-1) )!} z^{k-1} \frac{(\alpha + 1 + (n-1) )!}{(\alpha + 1 + (k-1) )! } \end{align} となる。ここで\( \nu = k-1 \) としてやると \begin{align} \frac{\partial L^{(\alpha) }_{n } (z) }{\partial z} = - \sum^{n-1}_{\nu = 0} \frac{ (-1)^{\nu} }{ \nu ! (n-1 - \nu) ! } z^{\nu } \frac{(\alpha + 1 + (n-1) )!}{(\alpha + 1 + \nu )! } = - L^{(\alpha + 1) }_{n-1} (z) \; \end{align} となることがわかる。