2019年09月16日

一次元調和振動子のシュレディンガー方程式の数値解に対して知見が広がった てへ

櫻井捷海 パーソナルコンピュータを用いた量子力学入門 裳華房 1989年 を参考にしてプログラムをいじっていたのである。著者の時代はパソコンはNECの国民機とそれの互換機がありDOSが主流だった。
本のBASICのコードをCで書いたり、gFortranで書いたり、Pythonで書いたりしたが結果はいずれも右側で発散が生じた。
一次元調和振動子のシュレディンガー方程式は2階微分方程式であり単純化された次のようなものだ。
y’’ = (t^2 – E) y
これを4次R-K法で数値計算で解く。この方程式は固有値・固有関数を求める問題でもあり数値解の中から固有値・固有関数を満たすものを探す。固有値固有関数を満たす場合には境界条件が満たされる、すまわち遠方での波動関数は0になる。そこには粒子はないからだ。
そこで左端の初期値を設定して右側に計算していく。固有値・固有関数を満たす場合には右側で0に収束するのである・・・・はずである…・収束してぇ!! 
 
そう、収束するのだが、さらに計算を続けると発散するのである。そのグラフは過去に何度か載せた。
私のプログラムは本のBASICを移植しただけである。gFortran, Cではそうだ。4次R-K法は単純な計算式だし。PythonではScipy というモジュールのodeint()を使った。これは4次R-K法の改良版だと思う。
そして、wxMaximaでも上の微分方程式の数値解を出して見た。wxMaximaには4次R-K法の関数が用意されている。そこでも結果は右側で発散した。

この数値計算で発散するのはどうしてか?   うーーむ・・・・4次R-K法の限界か?  うーーむ。

本では右側で収束すると直ぐに計算を打ち切っているのでグラフに出てこない。また、サイト記事で「物理のかぎしっぽ」にC言語で似たような計算しているのがあったが、これも右側で発散している。
微分方程式の数値解で右側に発散が見られる場合に、その発散は解ではなくアルゴリズムによるものなのか、それとも解であるのか?  そこの区別をどうするか。ノータリンのワタシにはわからないである(笑)

というときに別の解法を見つけた次第である。
https://qiita.com/sci_Haru/items/edb475a6d2eb7e901905#%E5%95%8F%E9%A1%8C%E8%A8%AD%E5%AE%9A

一次元シュレディンガー方程式の数値解を求めている。手法は、ヌメロフ法と書いてある。初耳である。不勉強が身に染みるである(´・ω・`)
  結果のグラフの例である。Pythonである。numpy + matplotlibを使っている。Scipyは使っていないである。

https://qiita-user-contents.imgix.net/https%3A%2F%2Fqiita-image-store.s3.amazonaws.com%2F0%2F192457%2Fae74add3-642f-8c0d-418d-dce797067954.png?ixlib=rb-1.2.2&auto=compress%2Cformat&fit=max&s=2d16b8172969b743d71ad38307a05255
このグラフを見る限り、横を20に増やしても発散しない。。。。気がする。

というわけで、ヌメロフ法とは?   上のサイトにリンクがあったのでクリックしたである。
https://qiita.com/sci_Haru/items/338a5bb68ce17189917b

この方程式の解法として,ヌメロフ法と呼ばれる非常に簡単かつ効率的な陽解法のアルゴリズがあ る。この方法は4次のルンゲ・クッタ法[1]よりも1次だけ精度が高い[2]。

なんて書いてあるである。すばらしい!!   1次だけ精度が高いというのであれば、4次R-K法を改良したバージョンもあるはずだなあ。それだと発散しないのかなあ。どうだろなあ。
ただ、2階微分方程式で1階尾分を含まない方程式で使えるというので・・・・なによ、それ。

では、別問題。Scipyにはodeint()の他にも積分器があるって話だが、ようするに各種のアルゴリズムを揃えているのか?  もし各種のアルゴリズムがあるのなら試してみたいである。
自分で書くの・・・・やだーーー。。。。。コピペだけどさ(笑)

以上、数値計算の入門あたりでウロウロしたりオロオロしたりしているワタクシの紹介でした。

posted by toinohni at 11:08| 東京 ☔| Comment(0) | ソフト系雑学 | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
コチラをクリックしてください