スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 | --/--/-- --:--

一次元クーロン場の固有状態の数値計算 (1)

解析的に解けない例として(たぶん解けないと思う)、計算してみたかったのが
一次元クーロン場の固有状態。

一次元クーロンポテンシャルはこんな感じのもの。\[
V(x) = -\frac{e^2}{4\pi\varepsilon_0 x}
\tag{1}
\]
こんなものは実際には存在しませんが、
3次元クーロンポテンシャルを疑似的に扱うのに使えるようです。

ただし、原点に特異性があると扱いにくいので、\[
V(x) = -\frac{e^2}{4\pi\varepsilon_0 \sqrt{x^2 + a^2}}
\tag{2}
\]と修正する( soft-core potential と呼ぶそうです)

原点の特異性が取り除かれて、
x >> a の遠方では (1) のクーロンポテンシャルに漸近する。
こうすることによって、原点を通過する電子の動きが計算できるとのこと。
確かに、便利そうですね!

さっそく、シュレディンガー方程式を立てる。\[
-\frac{\hbar^2}{2m} u'' - \frac{e^2}{4\pi\varepsilon_0 \sqrt{x^2+a^2}} u = Eu
\tag{3}
\]原子単位系を用いて、長さの単位をボーア半径\[
x_0 = \frac{4\pi\varepsilon_0 \hbar^2}{me^2}
\tag{4}
\]エネルギーの単位をハートリーエネルギー\[
E_0 = \frac{e^2}{4\pi\varepsilon_0 x_0}
= \frac{me^4}{(4\pi\varepsilon_0)^2 \hbar^2}
\tag{5}
\]としてスケーリングすると、\[
-\frac{1}{2} u'' - \frac{1}{\sqrt{x^2+a^2}} u = Eu
\tag{6}
\]と簡単化できる。u'' について解いて、\[
u'' = -2\left( E + \frac{1}{\sqrt{x^2+a^2}} \right) u
\tag{7}
\]を4次ルンゲクッタで計算する。


参考文献
[1] J. H. Eberly et al. "Numerical Experiments in Strong and Super-Strong Fields"
in M. Gavrila Ed. "Atoms in Intense Laser Fields"

スポンサーサイト
数値計算>量子力学 | コメント(0) | 2015/01/27 13:04

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

ちょっと本題からずれるのですが、
ついでに、もう少し高い励起状態までグラフにしてみました。
確率密度がどういう形になってるんだろう?と興味があったので。

数値計算ではなく、以下の解析解で計算しています。
規格化も入れています。\[
u_n(x) = (\pi^{1/2} 2^n n!)^{-1/2} H_n(x) e^{-x^2/2}
\tag{1}
\]

こちらが波動関数 $u(x)$ そのもの。
1DTISEHarmonicOscAnalyticEigenFunctions01.png

こちらが波動関数の2乗の確率密度 $|u(x)|^2$ をプロットしたもの。
1DTISEHarmonicOscAnalyticEigenFunctions02.png

なんとなく、イメージわきました(笑)

参考文献
[1] シッフ 「量子力学」(上)
数値計算>量子力学 | コメント(0) | 2015/01/27 12:38

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

固有値を求めるだけでなく、固有関数の方もちゃんと計算できてるか?
ということを確認。

初めの6個の固有状態を計算した結果。
青が偶関数で、赤が奇関数。規格化はしていません。
折り返した方 (x < 0) も書くと分かりやすいんでしょうけど、面倒なので・・・

1DTISEHarmonicOscEigenFunction01.png

固有エネルギーの分だけ下駄をはかせるお決まりの方法で表示してみました。

白い部分が古典的に許容される領域で、
影の部分は古典的には存在しえないが、
量子論的には波動関数が少し浸み出しているというのがこの表示だとよくわかります。

解析解\[
u(x) = a_n H_n(x) e^{-x^2/2}
\]の方も同様に書いてみて、比較してみました。
1DTISEHarmonicOscEigenFunction02.png

数値計算の方は規格化していないので、こちらの方の係数 $a_n$ を
数値計算の結果に合うように、適当に調整しています。

まあ、見た目ほぼ一致してるようですね。
どの程度の精度で一致しているかまでは評価してないけど、
まあこんなところでいいってことで・・・

ところで、実は調和振動子の計算をしたかったわけではありません。

解析解が分かるもので、正しく計算できることをテストしておいて、
次に、解析的には解けないようなポテンシャルで計算してみたいのです。

というわけで、次回からいよいよ解けないポテンシャルで計算してみようかと思います。
数値計算>量子力学 | コメント(0) | 2015/01/22 20:15

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

前に発散してしまうからという理由で、x = 2√E を無限遠に設定していたが、
前回の結果で、意外にも多少発散していても、零点は見つけられるということが分かったので、
無限遠を x = 30 に固定して、固有値の計算をやり直してみた。

これがその結果。

1.0000000001061
5.00000000687847
9.00000003899501
13.0000001166616
17.0000002600491
21.0000004893016
25.0000008245405
29.0000012858667
33.0000018933623
37.0000026670919
41.000003627103
45.0000047934266
49.0000061860783
53.0000078250579
57.0000097303504
61.0000119219255
65.0000144197385
69.0000172437302
73.0000204138268
77.0000239499403
81.0000278719684
85.0000321997946
89.0000369532885
93.0000421523053
97.0000478166863


固定したので、最初の数個もちゃんと精度が出ている。
ちなみに、x = 40 にしてみたら、さすがにやりすぎで、
低エネルギーの方では固有値が求められなかった。
(エネルギー掃引の刻みをもっと細かくすれば行けるかもしれないが・・・)

前回のように、固有値をグラフにしてみると、
1DTISEHarmonicOscEvaluesLinInterpInf30-01.png

相対誤差は、
1DTISEHarmonicOscEvaluesLinInterpInf30-02.png

今回は、すべての点で6桁以上の精度が得られている。

数値計算>量子力学 | コメント(0) | 2015/01/22 12:39

調和振動子の固有状態の数値計算 (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
 | HOME | Next »

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。