2019年09月15日

一次元調和振動子のシュレディンガー方程式の数値解で四苦八苦して苦しー!!

y'' = (v(t) - E)*y  これ。無次元化した後の微分方程式。2階微分方程式である。初期値を与え、4次R-K法で解く。数値計算する。v( t) = t^2
櫻井捷海 パーソルコンピュータを用いた量子力学入門 裳華房 1989年 を参照のこと

その本ではBASICでプログラムを書いておった。それを自分でc言語で、gFortranで、wxMaximaで、Pythonで・。・・・やってみた。Pythonでの結果の例を次に示す。
image

右側で発散しているのは無視する・・・としていたのだが気になった。wxMaxima, c, gFortranでも同じように右側では発散する。計算は左から始める。

で、これは4次R-K法でも誤差の蓄積によって発散に至るのだろう・・・・と思っていたが実はそうではないのではないか、と思い始めている。いや、誤差の蓄積なんだろうけどさ。
この微分方程式は解析的に解がでているので、それだと発散しないのだろ・・・確認したわけではないが。
y'' = (v(t) - E)*y  この微分方程式を y’’ = y という形に無理やり単純化すれば、解は簡単で、y = e^x  だ。それ指数関数なのでxが増えると発散するよ。
なのでxがある領域を超えると、この微分方程式の解は発散する・・・・と思った。解析的な解が発散しないのであれば数値解では何が原因でそうなるのか解明したい。
固有値・固有関数としてはE = 5は解析的にも解であって数値解でもグラフは左と右で収束していて。発散の前の0担っているところで計算を止めるといいので(笑)

この微分方程式は初期値を x= 0で0に近い値に設定して数値解を出せば発散するんでね?  と思った。ようするに4次R-K方だろうが修正したアルゴリズムだろうが発散するのだ。と思う。
x = 0で初期条件を0近くにして計算したらこうなった。
あるてころで、解がy = e^xの形に近づいているのでそうなるのではないか? 正負はともかくさ。

image

ここらの話は教科書てあれば右側で収束した領域で計算を止めているし、サイト記事もそういうのがヒットした・・・きがするぞ。

だが、左側から計算を始めて右側も収束した後も計算を続けると発散するのだ。wxMaxima, C, gFortranでも同じ結果。
発散まで計算した例は「物理のかぎしっぽ」に載っていた。Cで計算したとか。
というわけで、疑問は発散するのが本来の解ではないのか?   いや、解析解が発散はしないのだから数値計算上の問題だろ、その発散は・・・・ とか。
で、右側で微分方程式がy'' = (v(t) - E)*y  から y’’ = y という形に近似できる領域では発散して当然だと、さっきビール飲んでいて気づいたのワタクシ。

発散のように見えるのが本来の解であるという可能性はどうなのか? 

ボテンシャルが対称な場合には波動関数も対称になる・・・として右側で一度収束したら計算は止めるって妥当か?

この微分方程式の数値解で発散しない計算結果って、どっかにある?  


というわけで今の所、ワケワカメなのであります。賢人・先人の解説を探しても右側をずーと計算した例が「物理のかぎしっぽ」しか見つからなかった次第であります。
科学技術計算の達人に知り合いはいないしよ。
もっとも、ここらが理解できたとしても、それがどしたのって話だけどね。オレ、引きこもりでアル中 予備軍だし(笑)

posted by toinohni at 15:39| 東京 ☁| Comment(0) | 物理科学雑学 | このブログの読者になる | 更新情報をチェックする

分かりが遅いボクのために例題はたくさん・・という本やサイトを探す Python

クラスがまるで理解できない。分かりが遅いボクのためには例題がたくさん必要だ。たくさん見ているうちに、ある時点で、あ、そういうことか!!  と気づくとその一瞬が理解の瞬間だ。
サイトにしても本にしてもクラスを説明するために短いコードを出す。初心者向けとしてはいいのだろう。だが、何かがちんぷんかんププんがブン。
要するに、なんでこんなことするの、別の方法でもできるじゃん、関数定義してやれば簡単じゃん、・・・・・とかね。
原理的なところの説明のために短いコードを用意する著者の努力はわかるが、そんな約立たない例題ではなくてさ。多少、長くてもいいから実用になりそうな例を持って来いやーー てへ。

クラスの説明をカラスがカーと鳴くとか、カラスは黒いとか、そんな説明のものを見ると、オレとは次元が違うわな、こいつらは・・・・と思う次第である。

というわけだが、本を探す。図書館で。中古はAmazonで送料込みで500円未満を探す。

考えればわかることだがオレがプログラムを書くとしてもクラスを使うほどのものは書かない。せいぜい、100行前後。長くても500行ぐらいだ。クラスは使う側であり作る側ではないのだ、オレは・・・・・Python + Scipy + matplotlib で科学技術計算の入門みたいなところでアーーソブってのが趣味だし(笑)
クラスを理解して操れるようになろうって妄想であり、高望みであった ぎゃっふーーん!!

posted by toinohni at 09:08| 東京 ☀| Comment(0) | ソフト系雑学 | このブログの読者になる | 更新情報をチェックする