アトムの物理ノート

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

前回までにシュレディンガー方程式の固有値問題を解くという事について説明してきた。　私の知る限り、問題に応じて多少バリエーションはあるが、数値計算のアルゴリズムとしてはルンゲ･クッタやヌメロフ法で尽きている。　（ハミルトニアンを強引に対角化するという方法も有りえる。）　ところが、波動関数の時間発展を知りたいと思ったら、多少面倒な事に遭遇する。　

何が問題か？　差分方程式を与えられた初期条件で解く（時間発展させてゆく）だけではないか。数学的には全くその通りなのだが、life is more complex　と言っておこう。　論よりも証拠、先ず、シュレディンガー方程式を単純に離散化して解いてみた動画を見てもらおう。　問題はシュレディンガー方程式

i d/dt ψ(t,x) = [-∇² + V(x)] ψ(t,x)　;

V(x) ＝　 V₀ ( a < x < b )
＝　０　（その他の領域）

を離散化公式

∇² ψ(t,x)=dx^-2 [ ψ(t,x+dx) + ψ(t,x-dx) - 2ψ(t,x)]
dψ(t,x)/dt = dt^-1 [ ψ(t+dt,x) - ψ(t,x)]

を用いて解くことである。



（黄色はポテンシャル障壁がある領域）。　時刻ｔ＝０では適当に作った波束を用意した。　プロットは波動関数の２乗、つまり確率密度です。さて、この結果、おかしい。　どうおかしいって、有り得ない事が起こっている、確率が保存していないのである。　波動関数の裾がポテンシャル障壁に触れるまでは何とか形になっているのだが、暫くしてポテンシャル障壁にぶつかったとたんおかしな事になってくる。　実はこの離散化のアルゴリズムは最悪の欠点を持っているのである。　よってこの方法をあまり追求しても得るものがない（補足１）。　では何故私がこのアルゴリズムでプログラミングしたかというと、「うまくいかない事を実感するため」という目的があったからだ(笑）。　その問題点とは、シュレディンガー方程式を離散化するアルゴリズムがユニタリーにならないという事である。　よって量子力学で最も重要な確率の保存が成立しない。少し説明しよう。

シュレディンガー方程式は　idψ/dt　=　Hψ　と書ける。　この式の離散化バージョンは

ψ(t+dt) = ( 1 -i H dt) ψ(t)

である。　（1 -i H dt)≡U(dt) は波動関数に掛かって時間発展させる演算子だが

U(dt)^† U(dt) = ( 1 + i H dt) ( 1 -i H dt) = 1+ H² dt²

でユニタリー発展になっていない事が確認できる。　正確な時間発展演算子 U_exact=exp(-i H dt)　と比べると、シュレディンガー方程式を離散化して得られた U(dt) は dt,の一次までしか正しくない。　ここで多少注意しておかねばならないのはｄｔを十分小さくとれば高次の項は問題にならずに望む精度で時間発展が実現できると思うかもしれないが、現実の数値計算は常に誤差との戦いであり、理想的な極限操作は実現できないことである　（補足２）。よって何らかの改善を施さなければならないのだが、それは次回としよう。

（補足１）正直に明かすと、これがアルゴリズム自体の問題点なのか、それとも２次的な効果が利いているのか確信を持つにはいたらなかった。時間の刻み幅を相当に細かくしてみたがこの状況は改善しなかったため、最終的にこれはアルゴリズム自体が持つ問題点だろうと思ったわけだが、他の要因（数値誤差が増幅？　単なるバグ？）を完全には排除できなかっった。

（補足２）　理想論としてU_exactの演算子展開に現れるHⁿが離散化した系で小さく抑えられるという保障があれば、ｄｔを小さくとる極限で高次項は無視できるはずなのであるが、これが実際にはなかなか大変なのである。この辺は実際にコードを作ってみないと分からない部分であろう。