FC2ブログ

スポンサーサイト

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

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

一次元シュレディンガー方程式の数値計算の例として、
調和振動子の固有状態を計算してみようと思います。

解くべきシュレディンガー方程式は、\[
-\frac{\hbar^2}{2m} u'' - \frac{1}{2} m\omega^2 x^2 u = Eu
\tag{1}
\]
ここで、エネルギー E と 位置座標 x に関して、\[
E \rightarrow \frac{\hbar\omega}{2} E
\tag{2}
\]\[
x \rightarrow \left( \frac{\hbar}{m\omega} \right)^{1/2} x
\tag{3}
\]という無次元化を行うと、式は簡単になり、\[
u'' = (x^2 - E) u
\tag{4}
\]となる。この時、ポテンシャルも\[
V = x^2
\tag{5}
\]と無次元化される。

解析解はよく知られている通り、非負の量子数 n = 0, 1, 2, ・・・ に対して、\[
E = 2n + 1
\tag{6}
\]となる時だけ、無限遠で発散しない解となり、規格化因子を無視すると、\[
u(x) = H_n(x) \exp(-x^2/2)
\tag{7}
\]となる( $H_n$ はエルミートの多項式)。

まず、n = 0 すなわち E = 1 (基底状態)の近傍で、E を掃引しながら
4次ルンゲクッタ法で数値計算して、確認してみる。
偶関数が期待されるので、初期条件として、$u(0) = 1$、$u'(0) = 0$ とした。

1DTISEHarmonicOsc01.png

E = 0.5 から 0.1 間隔で E = 1.2 まで振ってみた。
確かに、E = 1.0 のところで発散しない解が得られているようだ。
n = 0 の解析解 $u = \exp (-x^2/2)$ (緑の点線)とも一致している。

今度は、n = 1 すなわち E = 3 (第一励起状態)の近傍で
奇関数の初期条件 $u(0) = 0$、$u'(0) = 1$ として計算してみる。
1DTISEHarmonicOsc02.png
やはり、 E = 3.0 のところで発散しない解が得られる。
n = 1 の解析解 $u = x \exp(-x^2/2)$ (緑の点線)とも一致している。

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

二階微分方程式の数値求解 (2)

とりあえず、クラスへの実装の話は置いておいて、
無骨にプログラミングして、計算してみました。

刻みを h = 0.1 に設定して、計算した結果。
偶関数の初期条件の場合。
RungeKuttaCos01.png
奇関数の初期条件の場合。
RungeKuttaSin01.png
赤が計算結果で、青が解析解。
比較してみると、きちんと計算できてそうです。

せっかくなので、刻みをどれだけ粗くできるかチェックしてみることに。

h = 0.2
RungeKuttaCos02.png

h = 0.5
RungeKuttaCos03.png

h = 1
RungeKuttaCos04.png

さすがに、ここまで粗くすると、少しずれが見えてきましたが、
結構、粗くても大丈夫なんですね。
振動周期の1/10程度なら大丈夫ということが分かりました。

というわけで、次回、シュレディンガー方程式の計算を試してみたいと思います。
数値計算>微分方程式 | コメント(0) | 2014/12/17 19:38

二階微分方程式の数値求解 (1)

シュレディンガー方程式を解くために、
まずは、二階の微分方程式を数値的に解く練習。

最も簡単なもので、\[
y'' = -y
\tag{1}
\]を解くことにする。

解析解は、奇関数の初期条件 $y(0) = 0$, $y'(0) = 1$ で\[
y = \sin x
\tag{2.1}
\]偶関数の初期条件 $y(0) = 1$, $y'(0) = 0$ で\[
y = \cos x
\tag{2.2}
\]
$y' = z$ とおくと、(1) は以下の連立微分方程式に帰着。\[
y' = z
\tag{3.1}
\]\[
z' = -y
\tag{3.2}
\]
計算には、4次ルンゲ・クッタ法を使うことにする。\[
y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \\
k_1 = f(x_n, y_n) \\
k_2 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h \right) \\
k_3 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h \right) \\
k_4 = f \left( x_n + h, y_n + k_3 h \right)
\tag{4}
\]
これを2関数の連立の場合\[
y' = f(x, y, z) \\
z' = g(x, y, z)
\tag{5}
\]に適用すると、\[
y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4) \\
z_{n+1} = z_n + \frac{h}{6} (l_1 + 2l_2 + 2l_3 + l_4)
\tag{6}
\]となり、各係数は、\[
k_1 = f(x_n, y_n, z_n) \\
l_1 = g(x_n, y_n, z_n) \\
k_2 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h, z_n + \frac{1}{2} l_1 h \right) \\
l_2 = g \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_1 h, z_n + \frac{1}{2} l_1 h \right) \\
k_3 = f \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h, z_n + \frac{1}{2} l_2 h \right) \\
l_3 = g \left( x_n + \frac{1}{2}h, y_n + \frac{1}{2} k_2 h, z_n + \frac{1}{2} l_2 h \right) \\
k_4 = f \left( x_n + h, y_n + k_3 h, z_n + l_3 h \right) \\
l_4 = g \left( x_n + h, y_n + k_3 h, z_n + l_3 h \right)
\tag{7}
\]
さて、このアルゴリズムでプログラミングしていくのですが、
以前にプログラミングの記事で書いていた
自作ウェーブ型クラスの DWave に
微分方程式求解機能をインプリメントしたら便利だろう
と思うのです。

たとえば、

double f(double x, double y);

DWave d1;
・・・ //d1 の初期設定

d1. RungeKutta4(1, f) //初期値 1, 関数 f

みたいな感じで呼び出せると便利!

と思ったのですが、考えてみたら、連立になっているので、
2値のベクトルとして扱わなければならないんですよね。

というわけで、

1.DWaeve をベクトルが扱えるように改良するか、
2.一般の連立ではなく、一関数の高階微分方程式だけに特化するか
3.DWave に実装するのはあきらめて、別でやるか。

1 は、すごく大変なので、すぐにやる気がしない^^;
2 もいいですが、とりあえず、DWave は置いておいて、3でやろうかと思っています。
数値計算>微分方程式 | コメント(0) | 2014/12/12 12:39

一次元シュレディンガー方程式の数値計算 (2)

というわけで、まずは、前記事1番の微分方程式を解く方法で挑戦してみます。

解くべき方程式は、\[
-\frac{\hbar^2}{2m} u''(x) + V(x) u(x) = Eu(x)
\tag{1}
\]
変形して、\[
u''(x) = -\frac{2m}{\hbar^2}\{ E-V(x) \} u(x)
\tag{2}
\]
ここで、\[
\alpha(x) \equiv \frac{2m}{\hbar^2}\{ E-V(x) \}
\tag{3}
\]と置いてしまうと、解くべき方程式は、\[
u''(x) = -\alpha(x) u(x)
\tag{4}
\]α(x) は長さの2乗の逆数(つまり、波数の2乗)の次元。

これを4次ルンゲ・クッタ法で解くことにしましょうか・・・

4次のルンゲクッタといえば、証明に難航しておりましたが、
まあ、公式は与えられているので、とりあえず、使うことにします。

2階微分まで含んでるから、処理が必要ですね。
$u' = v$ とおいて、(4) 式と合わせて、
\[
u' = v \\
v' = \alpha u
\tag{5}
\]という一階の連立微分方程式に帰着させる。

さて、初期条件をどうすればよいか。

とりあえず、原点に対して対称なポテンシャルのみを考えることにして、
パリティの決まった解のみを探せばいいんだろうか。

縮退がなければ、それでよさそう。
縮退があっても、2つの解を線形に組み合わせて、パリティの決まった解に持って行けたはず。
3つ以上とか無限に縮退している場合は、どうなるんだろう?

まあ、そんなことは後で考えよう・・・
とりあえず、縮退のない系を解くことにする。

縦方向のスケールは、後から規格化すれば任意でよいので、
$x=0$ における初期条件を
奇関数に対しては、\[
u_o(0) = 0, \ u_o'(0) = v_o(0) = 1
\tag{6.1}
\]偶関数に対しては、\[
u_e(0) = 1, \ u_e'(0) = v_e(0) = 0
\tag{6.2}
\]として、解けばよいのだろうか。

あまり自信がないですが、この方針で。
数値計算>量子力学 | コメント(0) | 2014/12/11 12:53

一次元シュレディンガー方程式の数値計算 (1)

以前、やりたいと言っていた一次元シュレディンガー方程式の数値計算

時間依存シュレディンガー方程式で波束の時間発展を見るプログラムは
クランク・ニコルソン法を用いて、以前に書いたことがあります。

今回は、時間に依存しないシュレディンガー方程式を解いて、
定常状態のエネルギー固有値と固有関数を求めたい。

2種類ぐらい方法がありそうです。

1.固有値をパラメータとして変化させながら、
  微分方程式をルンゲ・クッタなどで数値的に解いて、
  無限遠境界で発散しない解を探す。
  
  それと等価ですが、
  原点位置から解いたものと無限遠から解いたものが
  滑らかに接続するような解を探すのでもよい。

2.空間を離散化して、ハミルトニアンを行列表示する。
  この時、空間2階微分の項は三重対角(←たぶん)になるので、
  このハミルトニアンを対角化すれば、固有値と固有ベクトルが求まる。

2の方法にも興味はありますが、まずは1の方法を試したいと思っています。

文献 [1] を見ると、この2通りの方法の詳細が載っていますが、
まだちゃんと読んでません。


参考文献
[1] 倉澤治樹 千葉大講義資料「数値計算」
http://physics.s.chiba-u.ac.jp/~kurasawa/compute.pdf

数値計算>量子力学 | コメント(0) | 2014/12/11 12:09
« Prev | HOME | Next »

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