nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

二階微分演算子の変数変換 for Numerov method

プログラミング 数学

Numerov法は一階微分方程式を含まない二階微分方程式を解くのに便利な方法である。
この方法を適用するには、等間隔の変数刻み（メッシュ or グリッド）が必要である。
しかし、変化が穏やかなところでは粗いメッシュで十分だが、一方で変化が急なところでは細かなメッシュを用いなければ誤差が大きくなってしまうため、等間隔メッシュの使用は無駄が多い。
そこで、「細かいところは細かく、そうでないところは粗い、理想的なメッシュ」と対応する「等間隔なメッシュ」に変数変換し、変換後の方程式をNumerov法で解くことを考える。

変換前と後のメッシュを r, \rhoとして、 \rho = m[ r ] という関係で結ばれているとする。
\rhoは等間隔なので、最初の点を \rho_0として、
\displaystyle \rho_n = \rho_0 + n h \quad ( n = 0, 1, \cdots, N-1 )
と与えられる。 h, Nはそれぞれ \rhoのメッシュと点数である。
ここから、
\displaystyle r = m^{-1}[ \rho ] \quad ( m ( m^{-1}( \rho ) ) = \rho )
と逆変換したいところだが、関数系が複雑だと、解析的に書けない。
そこで、 y(r) = m[ r ] - \rho_n ととして、 y(r) = 0を与える rを数値的に求めて、それを r_nとして保存しておけば良い。
y(r) = 0を解くための数値計算法として、Newton法が簡便であるが、ここでは立ち入らない。
ニュートン法とは何か？？ニュートン法で解く方程式の近似解 - Qiita

次に、二階の微分演算子を変数変換する。
数値計算をするのは \rhoに対してだから、 \rhoでの二階微分から出発して、 rの二階微分とどう繋がるのかを見る。
df/dx = f'とすると、
\displaystyle \frac{d^2f}{d\rho^2} = \frac{d}{d\rho} \left( \frac{dr}{d\rho} f' \right) = \frac{d^2 r }{d \rho^2} f' + \frac{dr }{d \rho} \frac{df'}{d\rho} = \frac{d^2 r }{d \rho^2} f' + \left( \frac{dr }{d \rho} \right)^2 f''

重要なポイントとして、一次微分 f'が残っているとNumerov法が使えない。
そのため、 f(r) = a(r) g(r)として、上手いこと f'の項を消すような a(r)を考える
（ fを欲しい関数とすると上手く行かないから、本当に欲しい関数を gとして、そこに aを掛けた仮の関数を fと置き直すという流れ）。
準備として、
\displaystyle \frac{ dr }{ d \rho} = (\rho')^{-1} \\ \displaystyle \frac{ d^2r }{ d \rho^2} = \frac{ d }{ d \rho} (\rho')^{-1} = ( {\rho'} )^{-1} \left( (\rho')^{-1} \right)' = ( {\rho'} )^{-1} \left( - \frac{ \rho'' }{ (\rho')^2 } \right) = - \frac{ \rho'' }{ (\rho')^3 } \\ f' = a' g + a g' \\ f'' = a'' g + 2 a' g' + a g''

したがって、
\displaystyle \frac{d^2 f}{ d \rho^2} = \frac{d^2 r}{ d \rho^2 } ( a' g + a g' ) + \left( \frac{d r}{ d \rho } \right)^2 ( a'' g + 2 a' g' + a g'' ) \\ \displaystyle \quad = \left( \frac{d r}{ d \rho } \right)^2 a g'' + \left( \frac{d^2 r}{ d \rho^2} a + 2 \left( \frac{d r}{ d \rho } \right)^2 a' \right) g' + \left( \frac{d^2 r}{ d \rho^2} a' + \left( \frac{d r}{ d \rho } \right)^2 a'' \right) g

よって、 f'を残さないためには、
\displaystyle \frac{d^2 r}{ d \rho^2} a + 2 \left( \frac{d r}{ d \rho } \right)^2 a' = 0 \\ \displaystyle - \frac{ \rho'' }{ (\rho')^3 } a + \frac{2}{ (\rho')^2 } a' = 0 \\ \displaystyle \frac{ a' }{ a } = \frac{1}{2} \frac{ \rho'' }{ \rho' }
両辺を積分すれば、
\displaystyle \ln( a ) = \frac{1}{2} \ln( \rho' ) \\ \therefore \displaystyle a = \sqrt{ \rho' }
\rho=m[r]という形で \rhoが分かっているから、その微分は容易である。

今、欲しい関数は gに置き直していて、 gの二階微分方程式が g''(r) = b(r) g(r) + G(r)で与えられているとすれば、
\displaystyle \frac{ d^2 f }{d \rho^2} = \frac{ 1 }{ (\rho')^2 } a g'' + \left( - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a} \right) a g \\ \displaystyle = \frac{ 1 }{ (\rho')^2 } a ( b g + G ) + \left( - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a} \right) f \\ \displaystyle = ( c^{(1)}_g + c ) f + c^{(2)}_g G \\ \displaystyle c = - \frac{ \rho'' }{ (\rho')^3 } \frac{a'}{a} + \frac{1}{(\rho')^2} \frac{a''}{a} \\ \displaystyle c^{(1)}_g = \frac{ b }{ (\rho')^2 } \\ \displaystyle c^{(2)}_g = \frac{ a }{ (\rho')^2 }
c^{(1)}_g, c^{(2)}_gは元の微分方程式にあった係数が修正されたもので、一方 cは変数変換によって新たに付け足された項となる。
これをNumrov法で解けば良い。

具体例として、linear-log mesh \rho = \alpha r + \beta \ln r を使ってみる。
係数 cを手計算するのは面倒くさいので、pythonのsympyにやらせることにする。

from sympy import *

r, alpha, beta = symbols( "r α β" )

rho = alpha * r + beta * log( r )
rhop = Derivative( rho, r ).doit()
rhopp = Derivative( rhop, r ).doit()

a = sqrt( rhop )
ap = Derivative( a, r ).doit()
app = Derivative( ap, r ).doit()

c = simplify( - rhopp / (rhop**3) * ap / a + 1/(rhop**2) * app / a )

結果
\displaystyle c = \frac{ \beta ( \alpha r + \beta / 4 )}{ ( \alpha r + \beta )^4 }

sympyをもっと有効活用して行きたい。