ヌメロフ法

$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が大きくなるほど収束は遅くなることがわかる。
全体的に解析解と数値解はかなり近い値となっている。計算速度もかなり速かったと思う。