2017-09

目次

≪ WKBJ 3 ALL シュレディンガー方程式の数値計算 その3 ≫

スポンサーサイト

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

シュレディンガー方程式の数値計算 その2

schrodinger[pot[ポテンシャルの式、座標変数],{エネルギーの初期値、刻み幅},{ルンゲクッタのスタート点、波動関数、波動関数微係数}]

と入力します。例えば

schrodinger[pot[1/2 x^2、x],{-1, 0.001},{0, 1, 0}]

と入力すれば以下の出力が得られます。

調和振動子問題V(x)=1/2x^2のエネルギー固有値を求める計算です。基底状態は偶関数だと知っているので、初期条件は原点z=0での波動関数ψ(0)=1と微係数ψ'(0)=0を与えています。エネルギー固有値の可能性としてE=-1から始めて+0.001の刻み幅でシュレディンガー方程式を満足するかチェックしています。基底状態はe=1/2が答えで偶関数の第一励起状態はe=3/2が答えになるはずです。結果を見ると励起状態は誤差が大きいのですが、高度な事はするつもりがないので、これで満足しておきます。


schrodingerkeisankekka(1).gif

20070910082752.gif

20070910082737.gif

20070910082837.gif

20070910082818.gif

スペル・チェックをかけてませんので、早速おかしな単語を・・・・。見逃してください。


そのまま丸写しでも動作するはずです。適当にいじってもらって結構です。使用結果のフィードバックがあると嬉しいですが、必要条件では有りません。


===========マスマティカ・プログラム==========

schrodinger[pot[func_, x_], {esta_, estep_}, {z0point_, psi0val_, dpsi0val_}] :=
Module[{NumSol},
(* 初期値その他 *)
vpot[z_] := func /. x -> z;
e0 = esta;
de = estep;
zrange = 5;
mpoint = 3;
imax = 10000;
Print[" *** User supplied initial conditions ψ[", z0point, "]=", psi0val, " and ψ'[", z0point, "]=", dpsi0val];
Print[" *** Potential is v[z]=", vpot[z]];

(* ルンゲ・クッタ*)
PsiSol[z_, en_] := psi[z] /. NumSol[en, {z0point, psi0val, dpsi0val}, {zrange}];
NumSol[en_, {z0_, psi0_, dpsi0_}, {zwidth_}] := First[NDSolve[{(-2^(-1))*Derivative[2][psi][z] + (vpot[z] - en)*psi[z] == 0,
psi[z0] == psi0, Derivative[1][psi][z0] == dpsi0}, psi, {z, 0, zwidth}]];

(* z=∞での波動関数の制御:基底状態 *)
e = e0;
For[i = 0, i < imax && Sign[PsiSol[mpoint, e0]] === Sign[PsiSol[mpoint, e + de]], i++, e = e + de];
Print[" Result>> Ground State Energy = ", N[e], " ; starting e0= ", N[e0], " ; de step = ", N[de]]; firstenergy = N[e]; e0 = e + de;
(* z=∞での波動関数の制御:第一励起状態 *)
For[i = 0, i < imax && Sign[PsiSol[mpoint, e0]] === Sign[PsiSol[mpoint, e + de]], i++, e = e + de];
Print[" Result>> First Excited State Energy = ", N[e], " ;starting e0= ", N[e0], " ; de step = ", N[de]]; secondenergy = N[e];

(* グラフ出力 *)
Plot[{vpot[z], PsiSol[z, firstenergy], PsiSol[z, secondenergy]}, {z, 0, 2}, PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, PlotRange -> {{0, 2}, {-1, 2}}, Frame -> True,
FrameLabel -> {"z", "", StringJoin["\!\(ψ\_\(0, 1\)\) and V(z) for e= ", ToString[{firstenergy, secondenergy}]], ""}]; ]


コメント

コメントの投稿

管理者にだけ表示を許可する

トラックバック

http://letsphysics.blog17.fc2.com/tb.php/318-e7469440

«  | HOME |  »

CATEGORIES

RECENT ENTRIES

RECENT COMMENTS

RECENT TRACKBACKS

APPENDIX

アトム 

アトム 

趣味   近所散策と物理

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