Microsoft PowerPoint コンピュータ物理2_第13回.pptx

Transcription

1 コンピュータ物理学 第 3 回 06.. 第 回 0/ 金 ガイダンス 第 回 0/ 9 金 数値表現と誤差 第 3 回 0/6 金 第 4 回 0/3 金 数値微分 積分 第 5 回 0/30 木 第 6 回 /3 金 第 7 回 /0 金 常微分方程式 第 8 回 /7 金 第回 / 4 金 休講 第 9 回 / 金 常微分方程式 第 0 回 /8 金 偏微分方程式 第 回 /5 金 第 回 / 8 金 第 3 回 / 金 量子力学 第 4 回 / 9 金 最終レポート 第 5 回 / 5 金 モンテカルロ法

2 課題 4 解答例 Nma 0 N ma 00 36points 980points y 逐次 469 回 : δma= e 06 N ma points 逐次 847 回 : δma=9.9956e 06 収束の具合 y N ma 00 y 逐次 8045 回 : δma=9.998e 06 Nma 0 N ma 00

3 課題 4 解答例 課題 4 3 解答例 U,y Mes: iter=360 errorma=9.9839e 06 y y U,y y= 中心 +cm y=.5cm 中心面 cm

4 課題 4 3 補足 無限に広い電極なら同一平面上では同電位だが +5V 有限サイズの電極では. +5V 中心面 A での電位は 0 になるはずだが 実際は計算誤差が生じる 収束条件を厳しくすれば誤差は減らせる deltama=e 08 deltama=e 07 5V 3 3 5V y=+cm deltama=e 06 電極間内部 cm iter=739 deltama= e 07 電極内部でも電極端付近では一様性が悪い iter=8645 deltama=9.9844e 08 iter=9898 deltama= e 09 電極中心付近では一様性が良い cm

5 課題 4 3 プログラム補足 電極線上の電位 +5000V, 5000V は値が書き換えられてはならない 電極 電極上の点では電位計算を行わない wile errorma > e 6{ errorma = 0; forj=; j<=n_ma ; j++{ fori=; i<=n_ma ; i++{ ifi>=0 && i<=80&&j==0 j==80 continue; delta = U[i][j] 0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ]; U[i][j] = 0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ]; iffabsdelta>errorma errorma = fabsdelta; printf"iter=%d errorma=%e n", iter, errorma; iter++; 電極上の電位値を設定値にリセットする foriter=0;iter<=0000;iter++{ delta_ma=0.; fori=;i<=nma ;i++{ forj=;j<=nma ;j++{ delta=u[i][j] 0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ]; U[i][j]=0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ]; 値のリセット ifi>=0 && i<=80 { U[i][0]= 5000; U[i][80]=5000; ifj!=0&&j!=80 j==0 j==80&&i<0 i>80{ iffabsdelta> delta_ma delta_ma=fabsdelta; 収束判断は電極以外の場所のみで printf"%d %e n", iter, delta_ma; ifdelta_ma <.0e 6 brea;

6 例題 の program sample /* Partial Differential equation */ /* Sample problem; Poisson */ #include <stdio.> #include <mat.> main{ int i, j, iter; int N_ma=00; double, e0=8.85e, q=.0e 9; double errorma, corr; double U[N_ma+][N_ma+], ro[n_ma+][n_ma+]; FILE *outfile; car filename[0]; =.0/N_ma; /* bin widt range=0 */ /* set output file */ printf"output file?"; scanf"%s", &filename; outfile = fopenfilename, "w"; /* calculation */ for iter=0; iter<0000; iter++{ errorma = 0; forj=; j<=n_ma ; j++{ fori=; i<=n_ma ; i++{ corr = U[i][j] 0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ] ro[i][j]**/e0; U[i][j] = 0.5*U[i+][j]+U[i ][j]+u[i][j+]+u[i][j ]+ro[i][j]**/e0; iffabscorr>errorma errorma = fabscorr; printf"iter=%d errorma=%e n", iter, errorma; /* convergence cec */ iferrorma < e 5 brea; fori=0; i<=n_ma; i++{ forj=0; j<=n_ma; j++{ fprintfoutfile, "%f %f %f n", i*, j*, U[i][j]; /* Initialize and boundary condition */ forj=0; j<=n_ma; j++{ fori=; i<=n_ma; i++{ U[i][j]=0; ro[i][j]=0.0; ro[n_ma/][n_ma/]=q;

7 課題 4 4 /* y[i][j] : i= 座標, j= 時刻 =0,, */ /* 初期変位 y[i][0], y[i][] を代入 */ fori=0; i<8; i++{ y[i][0] = 0.005*i; /* i=80 */ fori=8; i<0; i++{ y[i][0] = *i 80; /* i=80 */ fori=0; i<00; i++{ y[i][] = y[i][0] + 0.5*63.0/00*63.0/00* y[i+][0]+y[i ][0] *y[i][0]; if%3==0{ fori=0; i<0; i++{ fprintfoutfile, "%d %d %f n",, i, y[i][]; if==75 fprintfoutfile, "%d %f n", i, y[i][]; /* 最初のつの時刻は初期値として書き出せる */ fori=0; i<00; i++ fprintfoutfile, "0 %d %f n", i, y[i][0]; fori=0; i<00; i++ fprintfoutfile, " %d %f n", i, y[i][]; /* 次々と時刻を進めて変位を計算していく */ for=; <ma; ++{ fori=; i<00; i++{ y[i][] =.0*y[i][] y[i][0]+63.0/00*63.0/00* y[i+][]+y[i ][].0*y[i][]; /* 最新の y[i][] が計算されたので y[i][] を y[i][0] に y[i][] を y[i][] に代入し直す 次のループで y[i][] を再び計算する */ fori=; i<00; i++{ y[i][0] = y[i][]; y[i][] = y[i][];

8 量子力学における数値計算 ; 次元井戸型ポテンシャル問題 原子レベルまたはそれ以下のサイズでの世界を扱う理論 特殊な系ではマクロサイズもあり得る 物理量について確率的に記述されるある粒子が ~ +d の間に発見される確率は粒子の状態 波動関数 ψ を用いて と表せる P d 粒子が定常状態であるとき 波動関数は時間に依存しない たとえば 中心にある原子核からのクーロンポテンシャル内に束縛された電子 原子核の一体場ポテンシャルに束縛された陽子 中性子 光電場にトラップされたイオンこのとき閉じ込められた粒子の波動関数は Scrodinger 方程式に従う : m d V E d ここで m は粒子の質量 V はポテンシャル E は粒子のエネルギー ある与えられた V のもとで 許される ψ と E を求めるのが解くべき課題である 粒子のエネルギーが量子化されているかを確認する ここでは粒子が井戸型ポテンシャルに存在するとする 陽子や中性子が原子核内に束縛されている場合等 d m d d d : E V0 a V0 a m V E a 0 a a V V 0 a

9 ポテンシャルの井戸の端, =±a において 粒子を発見する確率が連続で 粒子の流れも連続 である 波動関数とその微分が連続 である必要がある 接続条件 従って Ce a B cos a a と仮定して Ce a B cosa Ce a Bsina Ce a が必要となる tana 0 a, a とおくと tan 0 この両辺の比をとれば ここで 一方 ξ と η は を解くことに帰着する mv0a a V 0 m 83 MeV 940MeV 97MeVfm MeV - fm - ξ, η 平面上で この つのグラフを描いて交点を求めれば良い おおまかに求めると :,.5, 3.80, 3.60, 解が つ存在する 反対称な波動関数は考慮していない tan

10 数値計算で解く ; ヌメロフ法微分方程式の数値解は 4 次のルンゲクッタ法が推奨できるが Srodinger 方程式のように 次導関数項を含まない方程式の場合 方程式そのものを利用して極めて高精度に計算することができる ヌメロフ法と呼ばれ 実際に幅広く使われている Numerov 法の他に Fo Goodwin 法と呼ばれることもある 0 d d a a E m a a V E m, 0 まず Scrodinger 方程式を変形 : 4! 3!! 波動関数をテーラー展開する この式について つの微小区間 + について積分を行う際にポテンシャルを一定とみなす 4! 3!! O の周りでの展開も考える両者を足すと奇数次の項がキャンセルされる : : : 数値的に Scrodinger 方程式を解く ; ヌメロフ法 : :

11 これから 次の導関数は 4 4 O 階微分を中心差分で表現誤差 :~ 4 0 に上式を代入すると O しかし 4 階微分項はどう表せるのか? 0 4 d d そこで Srodinger 方程式に 階微分を作用させる : 4 d d もとの方程式 :

12 0 4 O 0 6 O 6 6 O O 5 6 O この 4 階微分の表現を先ほどの Scrodinger 方程式に代入すると + での ψ を計算するのに と - での値を使えば良い 前 ステップでの ψ の値 プログラムでは 離散的な空間座標 = i を使って = 刻み幅, i=0,,, 計算を進めていけば良い ルンゲクッタに比べても精度良く計算できる しかも式の数も少なく ψ+/ は計算不要 適用できる方程式の形が限定されているとはいえ 量子力学等において重要な計算方法である

13 5 この式は漸化式と見なせるので Runge=Kutta 法と同様に 変数領域の一方の端点で与えられた初期値から次々と方程式を積分していって解を求めることができる 実際の量子力学の固有値問題では 無限遠で波動関数が 0 に収束するような漸近条件が課されることが多い このような場合は まず適当なエネルギー固有値を与え 一方の端点から波動関数を計算していき 他方の端点での値が条件を満たすように エネルギー固有値を調整して再び波動関数を計算し直す というプロセスを繰り返す方法が良く使われる 0 において Numerov 法において波動関数の満たすべき関係式

14 井戸型ポテンシャル問題のヌメロフ法での解き方 解析的な手法では : ポテンシャル内とポテンシャル外との関数形を仮定して 接続条件を課してパラメータを決定する 数値的な解法では : 解析的手法と同様に上記の条件を利用する エネルギー E を推定してある値に仮定する 大きな正の で ψ= ep β であるとして出発し 中心に向かって だけステップを進めてヌメロフ法で波動関数を計算する 3 と並行して同じ微分方程式を 今度は大きな負の において ψ= epβ として出発し 中心に向かって だけステップを進める 4, 3 を交互に繰り返して進めていけば適当な位置 接続条件を判断する位置 において 左から来た解と右から来た解に整合性があるかどうかを判断する 4 整合性が悪ければ 分法を適用して整合性を満たす E の値を追い詰めていく からやり直す 5 左からの解と右からの解の値と傾きが許容範囲内で一致する場合に そのときの E の値が解 固有値 となり ただし今の場合 値が等しいという条件だけで十分 main 関数 : 分法 後述 により推定 E 値を追跡していく 計算が終わったら E と波動関数を出力する diff 関数 : main 関数から受け取った E を用いて 波動関数を左右両側からヌメロフ法により計算していき 整合性のチェックを行う 計算に出てくる は 関数から持ってくるのが便利 関数 : diff 関数内の計算から呼び出される の値に応じて 種類の値を返す plot 関数 : 固有値 E が決定されたらその E を使って全領域の波動関数を計算しなおしてファイル出力 main 関数の最後で呼び出す

15 例題 : 実際にプログラムを作成する 波動関数 のひとつ データを作成せよ 固有エネルギーのひとつを求めよ V V min = 000 ~ ma = 000 の区間 000 ステップ を考える 従って波動関数は 000 個の数値の並びとして表わされる i=0,,, と左から出発して波動関数 ψ left を計算する ヌメロフ法では 左から ψ の値を one, two, tree とおき one, two から tree を求める one, two には初期値が必要 例えば one = 0.0; two = とおいて良い 同様に右から i=0,, と波動関数 ψ rigt を計算する こちらも one, two, tree を今度は右から定義する ポテンシャルの深さは V 0 = 0.00 とする 整合性の良さを確認するための関数を用意する f E E E left rigt 右側 プラス側 のポテンシャル壁の位置で整合性のチェックを行うことにする 左側から = 0,,, 500. 右側から = 0,,, 500. と計算していき 重なる位置での両者の違い f を見る ここでは f が e 8 以下かどうかで判別する まずは E = と仮定して 波動関数を計算してみる 両側からの波動関数が接続しないことを確認する V 0 = 0.00 としたので エネルギー固有値の推定値として E min = 0.00 と E ma = を 分法の両端の初期値とし 整合性チェック時に左右の波動関数の差が e 8 以下になるまで 分法を繰り返す i= 0 i= 500 i= 500 m 簡単のため d 0 d E V E 0 a a, a a i= 0

16 分法によるゼロ点の探索 ある関数 fξ がゼロになる点 ξ = ξ 0 を効率的に探したい f ここを効率的に見つけたい 例えば f tan 0となるを探す ある区間 ξ - < ξ < ξ + において関数 fξ が符号を変えた場合 f f 0 f 区間 この条件を満たしながら ξ - と ξ + を徐々に接近させる c ξ と ξ+ の差がある許容限界 ε 以下になったとき 解に到達したとする 実際に行う作業は 各区間の中間点を考える c f c f 0 なら c そうでないなら c 区間 区間 3 c 区間 としてあらたな区間 [ξ, ξ+] を設定する このように区間次々と半分にし その都度どちらの側に解が含まれるかを調べ ξ - - ξ + < ε となるまで繰り返す

17 波動関数の計算 整合点 ポテンシャルから十分遠方では粒子に力が働かず 束縛状態においては の増加に伴い波動関数は ep β に従って減少する 大きな正の では Ψ = ep β であるとして中心に向 かって 左向きに だけステップを進める 大きな負の では Ψ = epβ として出発して中心に 向けて 右向き だけステップを進める 上記のつのステップを繰り返していけば いづれは 適当な位置で左から来た解と右から来た解に整合性 があるかを確かめられる 整合性が悪ければ 分法によりやり直し 整合性が十分 i= 0 i= 500 良ければそれを解とする 波動関数およびβ つまりE 今 波動関数の大きさは問題にしてないので 満たす逆べき整合条件はただ一つで良い 波動関数の値自体の連続性 i= 500 i= 0 [tree] [two] [one] i 0 i i 5 tree left one two rigt i i i 整合性のある解を探してエネルギー固有値を求めるだけならばわざわざ波動関数の値を配列にに記憶する必要はない しかし整合点での Ψleft と Ψrigt はそれぞれ独立に記憶させる必要はある

18 量子力学固有値問題計算プログラムの大枠 ; E を固定バージョン #define error e 8 #define V #define Emin #define Ema main{ E= ; まず E を固定してテストプログラム その E を使って波動関数を計算 ファイル出力 ; 関数 plot /* 波動関数計算中に ^ を代入するのに分かりやすくここで場合分けて計算する */ double int i, double E{ if i<500 return なんたら ; /* ポテンシャル外側 */ if i>=500 return なんたら ; /* ポテンシャル内側 */ V0 の符号に注意 /* 左からと右からの ψ を計算し 接続の良さを計算 */ double diffdouble E{ /* 左からの計算 */ fori=; i<=500; i++{ 波動関数の計算 ; tree = *****; 関数 の呼び出し i= 500 での値が決まる ; /* 右からの計算 */ fori=; i<=500; i++{ 波動関数の計算 ; tree = *****; 関数 の呼び出し i= 500 での値が決まる ; の値を main へ返す /* 正解の E で波動関数をファイル出力 */ void plotdouble E{ rigtfile = fopen rigt.dat, w ; leftfile = fopen left.dat, w ; /* 左からの計算 */ fori=; i<=500; i++{ 波動関数の計算 行ごとに fprintfleftfile, ; 座標に気を付ける /* 右からの計算 */ fori=; i<=500; i++{ 波動関数の計算 行ごとに fprintfleftfile, ; 座標に気を付ける

19 量子力学固有値問題計算プログラムの大枠 #define error e 8 #define V #define Emin #define Ema main{ do{ 分法で E の範囲を狭めていく関数 diff を呼び出して波動関数計算接続性が良くなるまで再トライ wile エネルギー E での接続が良いか?; 接続性が合格したらその時の E を printf. その E を使って波動関数を最後に計算 ファイル出力 ; 関数 plot /* 左からと右からの ψ を計算し 接続の良さを計算 */ double diffdouble E{ /* 左からの計算 */ fori=; i<=500; i++{ 波動関数の計算 ; tree = *****; 関数 の呼び出し i= 500 での値が決まる ; /* 右からの計算 */ fori=; i<=500; i++{ 波動関数の計算 ; tree = *****; 関数 の呼び出し i= 500 での値が決まる ; の値を main へ返す /* 波動関数計算中に ^ を代入するのに分かりやすくここで場合分けて計算する */ double int i, double E{ if i<500 return なんたら ; /* ポテンシャル外側 */ if i>=500 return なんたら ; /* ポテンシャル内側 */ V0 の符号に注意 /* 正解の E で波動関数をファイル出力 */ void plotdouble E{ rigtfile = fopen rigt.dat, w ; leftfile = fopen left.dat, w ; /* 左からの計算 */ fori=; i<=500; i++{ 波動関数の計算 行ごとに fprintfleftfile, ; 座標に気を付ける /* 右からの計算 */ fori=; i<=500; i++{ 波動関数の計算 行ごとに fprintfleftfile, ; 座標に気を付ける

20 整合性のテスト ポテンシャルの井戸の深さは 0.00 なので 基底状態のエネルギー固有値は 0.00~ の間にある可能性が大きい.. ということで 粗っぽく E= と E= を仮定したときの波動関数を計算してみると E= での Ψ E= での Ψ left rigt left rigt left rigt 両方とも整合条件とれてないが 固有 E 値はこの両者の間にあると思われる 今の問題では 左右の初期値の大きさを同じにとっているので 整合点でΨleft とΨrigt の値が同じになればその微分も同じになる left rigt