調和振動子のエネルギー固有状態の式についてオーバーフローを起こさない計算法

daruma3940.hatenablog.com

この前の続きです

この前は $(2^n n!)^{-1/2}$ をオーバーフローを起こさずに計算する方法についてと
それが量子的な調和振動子のエネルギー固有状態に出てくるという話をしました。
↓エネルギー固有状態の式
$\phi_n (x)= AHn( \xi )exp(- \xi^2 /2 )$
ここで
$\xi = \frac{x}{x_0}$
$A= \frac{C}{ \sqrt{n! 2^n} }$
としてみました。
Cとかx0が何かはここに書いてる
調和振動子 - Wikipedia

Hnはエルミート多項式なのですが、これがやべーやつ
$H_{n+1}(\xi )=2 \xi H_{n}(\xi ) +2 n H_{n-1}(\xi )$

なので $(2^n n!)^{-1/2}$ が計算できるようになっても結局オーバーフローします
さっきのΦについての式をHnについて書き直してみましょう
$Hn( \xi )= \sqrt{ \frac{n! 2^n}{C} } exp( \xi^2 /2) \phi_n(x)$
これをさっきの
$H_{n+1}(\xi )=2 \xi H_{n}(\xi ) -2 n H_{n-1}(\xi )$
に代入してまとめると
$\phi_{n+1} =\sqrt{ \frac{2}{n+1}} \xi \phi_{n} - \sqrt{ \frac{n}{n+1}} \phi_{n-1}$

なのでn=0 とn=1について計算すれば計算出来ちゃいます

def zenka_oscn(n, x):
    if (n == 1 or n == 0):
        return occeogenfanclist[n](x)
    else:
        return ((2 / n) ** 0.5) * x/x0 * zenka_oscn(n - 1, x) - (((n - 1) / n) ** 0.5) * zenka_oscn(n - 2, x)

これを単純に計算しようとすると重いのでフィボナッチ数列の場合のようにメモ化（キャッシュ化？）してみましょう
pythonのメモ化はクッソ簡単です
頭に@lru_cache(maxsize=hoge)をつけるだけです。やっぱりpython便利

from functools import lru_cache

@lru_cache(maxsize=None)
def zenka_oscn(n, x):
    if (n == 1 or n == 0):
        return occeogenfanclist[n](x)
    else:
        return ((2 / n) ** 0.5) * x/x0 * zenka_oscn(n - 1, x) - (((n - 1) / n) ** 0.5) * zenka_oscn(n - 2, x)