FC2ブログ

調和振動子の固有状態の数値計算 (4)

自動的に直線補間して、固有値を一気に求めるプログラムができた。
エネルギーを掃引しながら、無限遠での零点を見つけては、
直線補間して、零点の正確な値を出力していく。

零点の求め方は、まず、符号が反転したところを見つける。
順次エネルギーを掃引していき、$E_k$ と$E_{k+1}$ に対して、
$u_k u_{k+1} < 0$ となったら、符号が反転していると判定できる。

その時に、直線補間で零点の正確な値を求める。
その値を E とすると、\[
E = \frac{u_{k+1} E_k - u_k E_{k+1}}{u_{k+1} - u_k}
\]
1DTISEHarmonicOscEValuesLinInterpConcept01.png

以上のアルゴリズムで、E = 0.5~100 の範囲で求めた固有値一覧がこちら。

1.06019474008538
5.00001293981931
9.00000003919728
13.0000001157887
17.00000025828
21.0000004862375
25.0000008197467
29.0000012788781
33.0000017062101
37.0000026542532
41.0000036105981
45.0000047727627
49.0000054600211
53.0000077946515
57.0000085103253
61.000011879947
65.0000143713303
69.0000171885402
73.0000203514719
77.0000205266203
81.0000277945322
85.0000273957349
89.0000368597403
93.0000420505301
97.0000477065965


エネルギーの刻みは前と同様、0.01に設定し、
無限遠も前と同様、 x = 2√E に設定した。

僕の使っているコンパイラでは、double 型は仮数部が 52 bit らしいので、
$\log_{10} 2^{52} \simeq 15.7$ より
一応、15 桁ぐらいが有効だろうということで、
出力ストリームで precision(15) と指定して出力している。
(なぜか、15桁になっていないところがあるのは謎)

E > 81 で発散していると思っていたが、ちゃんと E = 100 まで計算できている。
前は、振動の振幅が大きすぎて、グラフソフトがうまく表示できていなかっただけらしい。

求めた固有値をグラフにすると、
1DTISEHarmonicOscEvaluesLinInterp01.png
表を見ても明らかであるが、E = 4n + 1 の理論直線上に並んでいるのが分かる。
(偶関数のみで計算しているので、4n + 3 の系列は現れない)

次に、理論値からの相対誤差 $|\Delta E|/E$ を調べてみた。
1DTISEHarmonicOscEvaluesLinInterp02.png
始めの2点の誤差が大きいのは、無限遠の値 2√E が小さすぎるためだろう。
それ以外は、6桁以上の精度が出ているので、とりあえず、十分。
スポンサーサイト



数値計算>量子力学 | コメント(0) | 2015/01/21 19:20
コメント

管理者のみに表示