コンピュータ物理学２
第13回
(2016.1.22)
第1回
10/ 2(金) ガイダンス
第2回
10/ 9(金) 数値表現と誤差
第3回
10/16(金) “
第4回
10/23(金) 数値微分・積分
第5回
10/30(木) “
第6回
11/13(金) “
第7回
11/20(金) 常微分方程式
第8回
11/27(金)
“
第 回
12/ 4(金) 休講
第 9回 12/11(金) 常微分方程式
第10回 12/18(金) 偏微分方程式
第11回 12/25(金) “
第12回
1/ 8(金) “
第13回
1/ 22(金) 量子力学
第14回
1/ 29(金) “
（最終レポート）
第15回
2/ 5(金) モンテカルロ法
20
max
N
N
max
100
200
max
N
20
max
N
N
max
100
200
max
N
収束の具合
課題4‐(1) 解答例
逐次 469回 : δmax= 9.996688e‐06
361points
9801points
39601points
逐次 8427回: δmax=9.995611e‐06
逐次28045回: δmax=9.998112e‐06
x
y
x
y
x
y
Mesh: 100 x 100
iter=360 errorMax=9.983922e‐06
課題4‐(2) 解答例
y=2.5cm （中心面）
y=中心+1cm
課題4‐(3) 解答例
U(x,y)
x
y
x
y
x (cm)
U(x,y)
iter=7392 deltaMax=9.985374e‐07
iter=8645 deltaMax=9.984241e‐08
iter=9898 deltaMax=9.983107e‐09
課題4‐(3) 補足
deltaMax=1e‐06
deltaMax=1e‐07
deltaMax=1e‐08
y=+1cm
中心面(A) での電位は 0 になるはずだが…
実際は計算誤差が生じる
→収束条件を厳しくすれば誤差は減らせる
無限に広い電極なら
同一平面上では同電位だが…
電極間内部
電極内部でも電極端付近 では一様性が悪い 電極中心付近 では一様性が良いx (cm)
x (cm)
+15kV ‐15kV +15kV ‐15kV 3 2 1
1
2
3有限サイズの電極
では….
while (errorMax > 1e‐6){ errorMax = 0; for(j=1; j<=N_max‐1; j++){ for(i=1; i<=N_max‐1; i++){ if((i>=20 && i<=80)&&((j==20)||(j==80))) continue; delta = U[i][j]‐0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1]); U[i][j] = 0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1]); if(fabs(delta)>errorMax) errorMax = fabs(delta); } } printf("iter=%d errorMax=%e ¥n", iter, errorMax); iter++; }
(1) 電極上の点では電位計算を行わない
for(iter=0;iter<=10000;iter++){ delta_max=0.; for(i=1;i<=Nmax‐1;i++){ for(j=1;j<=Nmax‐1;j++){ delta=U[i][j]‐0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1]); U[i][j]=0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1]); if(i>=20 && i<=80) { U[i][20]=‐15000; U[i][80]=15000; } if(((j!=20)&&(j!=80)) || ( ((j==20)||(j==80))&&((i<20)||(i>80)))){ if((fabs(delta)> delta_max)) delta_max=fabs(delta); } } } printf("%d %e¥n", iter, delta_max); if(delta_max < 1.0e‐6) break; }(2) 電極上の電位値を設定値にリセットする
課題4‐(3) プログラム補足
電極線上の電位 (+15000V, ‐15000V)は値が書き換え
られてはならない
電極
値のリセット 収束判断は電極以外の場所のみで。/* calculation */
for (iter=0; iter<10000; iter++){ errorMax = 0;
for(j=1; j<=N_max‐1; j++){ for(i=1; i<=N_max‐1; i++){
corr = U[i][j] ‐ 0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1])‐rho[i][j]*h*h/e0; U[i][j] = 0.25*(U[i+1][j]+U[i‐1][j]+U[i][j+1]+U[i][j‐1])+rho[i][j]*h*h/e0; if(fabs(corr)>errorMax) errorMax = fabs(corr); } } printf("iter=%d errorMax=%e ¥n", iter, errorMax); /* convergence check */ if(errorMax < 1e‐5) break; } for(i=0; i<=N_max; i++){ for(j=0; j<=N_max; j++){ fprintf(outfile, "%f %f %f ¥n", i*h, j*h, U[i][j]); } } } /* Partial Differential equation */ /* Sample problem; Poisson */ #include <stdio.h> #include <math.h> main(){ int i, j, iter; int N_max=100; double h, e0=8.85e‐12, q=1.0e‐9; double errorMax, corr; double U[N_max+1][N_max+1], rho[N_max+1][N_max+1]; FILE *outfile; char filename[20]; h=1.0/N_max; /* bin width (x‐range=0‐1) */ /* set output file */ printf("Output file ?"); scanf("%s", &filename); outfile = fopen(filename, "w"); /* Initialize and boundary condition */ for(j=0; j<=N_max; j++){ for(i=1; i<=N_max; i++){ U[i][j]=0; rho[i][j]=0.0; } } rho[N_max/2][N_max/2]=q;
例題2 の program sample
課題4‐(4)
/* y[i][j] ： i=x座標, j=時刻(=0,1,2) */ /* 初期変位 y[i][0], y[i][1] を代入 */ for(i=0; i<81; i++){ y[i][0] = 0.00125*i; /* y=0.1 @ i=80 */ } for(i=81; i<101; i++){ y[i][0] = 0.1 ‐ 0.005*(i‐80); /* y=0.1 @ i=80 */ } for(i=0; i<100; i++){ y[i][1] = y[i][0] + 0.5*(63.0/100)*(63.0/100)* (y[i+1][0]+y[i‐1][0]‐2*y[i][0]); } /* 最初の２つの時刻は初期値として書き出せる */ for(i=0; i<100; i++) fprintf(outfile, "0 %d %f¥n", i, y[i][0]); for(i=0; i<100; i++) fprintf(outfile, "1 %d %f¥n", i, y[i][1]); /* 次々と時刻を進めて変位を計算していく */ for(k=2; k<kmax; k++){ for(i=1; i<100; i++){ y[i][2] = 2.0*y[i][1]‐y[i][0]+(63.0/100)*(63.0/100)* (y[i+1][1]+y[i‐1][1]‐2.0*y[i][1]); }/* 最新の y[i][2] が計算されたので y[i][1] を y[i][0] に、 y[i][1] を y[i][2] に代入し直す → 次のループで y[i][2] を 再び計算する */ for(i=1; i<100; i++){ y[i][0] = y[i][1]; y[i][1] = y[i][2]; } if(k%3==0){ for(i=0; i<101; i++){ fprintf(outfile, "%d %d %f¥n", k, i, y[i][2]); if(k==75) fprintf(outfile2, "%d %f ¥n", i, y[i][2]); } } }
量子力学における数値計算； 1次元井戸型ポテンシャル問題
・ 原子レベルまたはそれ以下のサイズでの世界を扱う理論（特殊な系ではマクロサイズもあり得る）
・ 物理量について確率的に記述される
ある粒子が
x ～ x+dx の間に発見される確率は粒子の状態（波動関数）ψ(x) を用いて
dx
x
x
P
(
)
(
)
2)
(
)
(
)
(
)
(
2
2 2 2x
E
x
x
V
dx
x
d
m
粒子が定常状態であるとき、波動関数は時間に依存しない。たとえば
・ 中心にある原子核からのクーロンポテンシャル内に束縛された電子
・ 原子核の一体場ポテンシャルに束縛された陽子・中性子
・ 光電場にトラップされたイオン
このとき閉じ込められた粒子の波動関数は Schrodinger 方程式に従う：
ここで m は粒子の質量、V(x) はポテンシャル、E は粒子のエネル
ギー。ある与えられた V(x) のもとで、許される ψ(x) と E を求めるの
が解くべき課題である。粒子のエネルギーが量子化されているか
を確認する。ここでは粒子が井戸型ポテンシャルに存在するとする
（陽子や中性子が原子核内に束縛されている場合等）
x
a
x
E
m
dx
x
d
a
x
x
V
E
m
dx
x
d
)
(
2
)
(
)
(
2
)
(
2 2 2 0 2 2 2
a
x
a
x
V
x
V
0
)
(
0a
a
V
0V
x
と表せる。
:ポテンシャルの井戸の端, x=±a において、「粒子を発見する確率が連続で、粒子の流れも連続」である
→ 「波動関数とその微分が連続」である必要がある （接続条件）
従って、
a aCe
a
B
Ce
a
B
sin
cos
0
tan
a
が必要となる
この両辺の比をとれば
a
a
,
とおくと、
tan
0
を解くことに帰着する。
2
0
.
04829
83
2
216
.
032
2 2 0 2 2 2 2 2
a
mV
a
一方、ξ と η は
ξ, η 平面上で、この２つのグラフを描いて
交点を求めれば良い。おおまかに求めると：
ここで
,
1
.
25
,
3
.
80
,
3
.
60
,
1
.
74
解が２つ存在する。（反対称な波動関数は考慮
していない）
tan
03
.
16
2 2
a
x
a
x
a
a
x
Ce
x
B
Ce
x
x x
(
)
cos
と仮定して
2 -1 -2 2 fm MeV 04829 . 0 fm 197MeV MeV 940 2 2 mMeV
83
0
V
数値計算で解く； ヌメロフ法
微分方程式の数値解は４次のルンゲクッタ法が推奨できるが、Shrodinger方程式のように1次導関
数項を含まない方程式の場合、方程式そのものを利用して極めて高精度に計算することができる。
ヌメロフ法と呼ばれ、実際に幅広く使われている。（Numerov 法の他に、Fox‐Goodwin 法と呼ばれる
こともある。）
0
)
(
)
(
)
(
2 2 2
k
x
x
dx
x
d
x
a
a
x
E
m
a
x
a
V
E
m
k
,
2
)
(
2
2 0 2 2
まず Schrodinger 方程式を変形：
)
(
!
4
)
(
!
3
)
(
!
2
)
(
)
(
)
(
x
h
x
h
(1)x
h
2
(2)x
h
3
(3)x
h
4
(4)x
波動関数をテーラー展開する。
この式について、１つの微小区間 x → x+h について積分を行う際にポテンシャルを一定と
みなす。
)
(
!
4
)
(
!
3
)
(
!
2
)
(
)
(
)
(
x
h
x
h
(1)x
h
2
(2)x
h
3
(3)x
h
4
(4)x
)
(
)
(
12
)
(
)
(
2
)
(
)
(
x
h
x
h
x
h
2
(2)x
h
4
(4)x
O
h
6
x‐h の周りでの展開も考える
両者を足すと奇数次の項がキャンセルされる：
: :数値的にSchrodinger 方程式を解く； ヌメロフ法
: :これから２次の導関数は
)
(
)
(
12
)
(
2
)
(
)
(
)
(
2 (4) 4 2 ) 2 (h
x
O
h
h
x
h
x
h
x
x
2階微分を中心差分で表現
誤差：～
h
40
)
(
)
(
)
(
2 ) 2 (x
k
x
x
に上式を代入すると、
0
)
(
)
(
)
(
12
)
(
2
)
(
)
(
2 (4) 4 2 2
x
k
h
O
x
h
h
x
h
x
h
x
しかし、4階微分項はどう表せるのか？
(
)
(
)
0
)
(
2 2 2 ) 4 (
k
x
x
dx
d
x
そこで Shrodinger方程式に２階微分を作用させる：
2 2 2 2 2 2 2 2 ) 4 ()
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
h
h
x
h
x
k
x
x
k
x
x
k
h
x
h
x
k
x
x
k
dx
d
x
もとの方程式
:
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
0
12
)
(
2
)
(
)
(
4 2 2 2 2 2 2 2 2
x
k
h
O
h
h
x
h
x
k
x
x
k
x
x
k
h
x
h
x
k
h
h
x
h
x
h
x
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
0
12
)
(
2
)
(
)
(
x
h
x
h
x
h
2k
2x
h
x
h
k
2x
x
k
2x
x
k
2x
h
x
h
O
h
6
h
2k
2
x
)
(
)
(
)
(
12
1
)
(
)
(
)
(
6
2
)
(
)
(
12
1
h
2k
2x
h
x
h
h
2k
2x
h
2k
2x
x
h
2k
2x
h
x
h
O
h
6
)
(
)
(
)
(
12
1
)
(
)
(
6
5
2
)
(
)
(
12
1
h
2k
2x
h
x
h
h
2k
2x
x
h
2k
2x
h
x
h
O
h
6
)
(
)
(
12
1
)
(
)
(
12
1
1
)
(
)
(
12
5
1
2
)
(
6 2 2 2 2 2 2h
O
h
x
k
h
h
x
h
x
k
h
x
x
k
h
h
x
この４階微分の表現を先ほどの Schrodinger 方程式に代入すると
・
x+h
での
ψ
を計算するのに
x
と
x-h
での値を使えば良い（前２ステップでの
ψ
の値）。
・ プログラムでは、離散的な空間座標
x = ih
を使って(h=刻み幅, i=0,1,2,・・・) 計算を
進めていけば良い。
・ ルンゲクッタに比べても精度良く計算できる。しかも式の数も少なく、ψ(x+h/2) は計算不要。
→ 適用できる方程式の形が限定されているとはいえ、量子力学等において重要な計算方法である。
)
(
12
1
)
(
)
(
12
1
1
)
(
)
(
12
5
1
2
)
(
2 2 2 2 2 2h
x
k
h
h
x
h
x
k
h
x
x
k
h
h
x
・ この式は漸化式と見なせるので、Runge=Kutta 法と同様に、変数領域の一方の端点で与えられた
初期値から次々と方程式を積分していって解を求めることができる。
・ 実際の量子力学の固有値問題では、無限遠で波動関数が 0 に収束するような漸近条件が課される
ことが多い。
・ このような場合は、まず適当なエネルギー固有値を与え、一方の端点から波動関数を計算していき、
他方の端点での値が条件を満たすように、
エネルギー固有値を調整して再び波動関数を計算し直す
、
というプロセスを繰り返す方法が良く使われる。
0
)
(
x
x
において
Numerov法において波動関数の満たすべき関係式
井戸型ポテンシャル問題のヌメロフ法での解き方
解析的な手法では：
ポテンシャル内とポテンシャル外との関数形を仮定して、接続条件を
課してパラメータを決定する。
数値的な解法では：
解析的手法と同様に上記の条件を利用する。
(1) エネルギー E を推定してある値に仮定する。
(2) 大きな正の x で ψ(x)= exp(‐βx) であるとして出発し、中心に向かって h だけ
ステップを進めてヌメロフ法で波動関数を計算する。
(3) (2)と並行して同じ微分方程式を、今度は大きな負の x においてψ(x)= exp(βx)
として出発し、中心に向かって h だけステップを進める。
(4) (2), (3) を交互に繰り返して進めていけば適当な位置（接続条件を判断する
位置）において、左から来た解と右から来た解に整合性があるかどうかを判
断する。
(4) 整合性が悪ければ２分法を適用して整合性を満たす E の値を追い詰めて
いく。（ (1) からやり直す。）
(5) 左からの解と右からの解の値と傾きが許容範囲内で一致する場合に、その
ときのE の値が 解（固有値）となり、。
※ ただし今の場合、値が等しいという条件だけで十分。
main関数： 2分法（後述）により 推定 E 値を追跡していく。計算が終わったら E と波動関数を出力する。
diff 関数：
main関数から受け取った E を用いて、波動関数を左右両側からヌメロフ法により計算
していき、整合性のチェックを行う。計算に出てくる k
2(x) は k2 関数から持ってくるのが便利。
k2 関数：
diff 関数内の計算から呼び出される。x の値に応じて2種類の値を返す。
plot 関数： 固有値 E が決定されたらその E を使って全領域の波動関数を計算しなおしてファイル出力。
main 関数の最後で呼び出す。
例題：実際にプログラムを作成する
・ 波動関数（のひとつ）データを作成せよ
・ 固有エネルギーのひとつを求めよ
・ xmin= ‐1000 ~ xmax= 1000 の区間(2000ステップ)を考える。 従って波動関数は 2000 個の数値の並びとして表わされる。 ・ i=0, 1,2,・・ と左から出発して波動関数ψleftを計算する。ヌメロフ 法では、左からψ の値を one, two, three とおき、one, two から three を求める（one, two には初期値が必要（例えば one = 0.0; two = 0.00001 とおいて良い）。 同様に右から i=0, 1, 2 ・・ と波動関数ψrightを計算する。こちらも one, two, three を今度は右から定義する。 ・ ポテンシャルの深さは V0= ‐0.001 とする。 ・ 整合性の良さを確認するための関数を用意する。 右側（プラス側）のポテンシャル壁の位置で整合性のチェックを 行うことにする。 左側からx = 0,1,2, … 1500. 右側から x = 0,1, 2, … 500. と計算して いき、重なる位置での両者の違い f を見る。 （ここでは f が 1e‐8 以下かどうかで判別する。） ・ まずは E = ‐0.0009 と仮定して、波動関数を計算してみる。 （両側からの波動関数が接続しないことを確認する）。 ・ V0=‐0.001 としたので、エネルギー固有値の推定値として Emin= ‐0.001 と Emax= ‐0.00085 を2分法の両端の初期値とし、 整合性チェック時に左右の波動関数の差が 1e‐8 以下になるまで 2分法を繰り返す。500
500
500
0
500
)
(
0x
x
V
x
V
1000
1000
i = 0
i = 1500
i = 0
i = 500
)
(
)
(
)
(
E
E
E
f
left
right簡単のため
2
1
2
m
0
)
(
)
(
)
(
2 2 2
k
x
x
dx
x
d
x
a
a
x
E
a
x
a
V
E
k
,
0 2２分法によるゼロ点の探索
ある関数
f(ξ)
がゼロになる点 ξ = ξ
0を効率的に探したい。
ある区間
ξ
-< ξ < ξ
+において関数
f(ξ)
が符号を変えた場合
0
)
(
)
(
f
f
この条件を満たしながら
ξ
-と
ξ
+を徐々に接近させる
ξ‐ と ξ+ の差がある許容限界ε 以下になったとき、解に到達
したとする
実際に行う作業は
各区間の中間点を考える
2
1
c
c
区間1
区間2
c
区間3
区間2
0
)
(
)
(
f
f
c
cそうでないなら
cなら
としてあらたな区間 [ξ‐, ξ+] を設定する。
このように区間次々と半分にし、その都度どちらの側に解が含まれるかを調べ、
| ξ
-- ξ
+| < ε
となるまで繰り返す。
を探す
となる
例えば
f
(
)
tan
0
ここを効率的に 見つけたい)
(
f
)
(
f
1
i
0
i
i
2
) ( 12 1 two ) ( 12 1 1 one ) ( 12 5 1 2 three 2 2 2 2 2 2 h x k h h x k h x k h 500
500
1000
1000
i = 0
i = 1500
i = 0
i = 500
x
・ ポテンシャルから十分遠方では粒子に力が働かず、
束縛状態においては、|x| の増加に伴い波動関数は
exp(‐β|x|) に従って減少する。
・ 大きな正の x では Ψ = exp(‐βx) であるとして中心に向
かって（左向きに） h だけステップを進める
・ 大きな負の x では Ψ = exp(βx) として出発して中心に
向けて（右向き） h だけステップを進める。
・ 上記の２つのステップを繰り返していけば、いづれは
適当な位置で左から来た解と右から来た解に整合性
があるかを確かめられる
・ 整合性が悪ければ2分法によりやり直し、整合性が十分
良ければそれを解とする（波動関数およびβ(つまりE)）。
・ 今、波動関数の大きさは問題にしてないので、満たす
べき整合条件はただ一つで良い。（波動関数の値自体の連続性）
0
.
0
0
0.00001
1
2
1 i i0 2 i0
.
0
0
0.00001
1
2
left
right ※ 整合性のある解を探してエネルギー固有値を求めるだけならばわざわざ波動関数の 値を配列にに記憶する必要はない。しかし整合点でのΨleft とΨright はそれぞれ独立 に記憶させる必要はある。波動関数の計算
[one] [two] [three]整合点
逆
main(){ E= ‐0.0009; ← まず E を固定してテストプログラム。 その E を使って波動関数を計算・ファイル 出力;関数plot } /* 左からと右からのψを計算し、接続 の良さを計算*/ double diff(double E){ /* 左からの計算 */ for(i=1; i<=1500; i++){ 波動関数の計算； three = *****; 関数k2の呼び出し } i = 1500 での値が決まる; ① /* 右からの計算 */ for(i=1; i<=500; i++){ 波動関数の計算; three = *****; 関数k2の呼び出し } i = 500 での値が決まる；② ①‐② の値を main へ返す } /* 波動関数計算中に k^2 を代入するのに分かりやすく ここで場合分けて計算する */
double k2(int i, double E){
if (i<500) return(なんたら); /* ポテンシャル外側*/ if (i>=500) return(なんたら); /* ポテンシャル内側*/ (V0 の符号に注意) } /* 正解のE で波動関数をファイル出力 */ void plot(double E){ rightfile = fopen(“right.dat”, “w”); leftfile = fopen(“left.dat”, “w”); /* 左からの計算 */ for(i=1; i<=1500; i++){ 波動関数の計算 1行ごとにfprintf(leftfile, ) ; (x座標に気を付ける） } /* 右からの計算 */ for(i=1; i<=500; i++){ 波動関数の計算 1行ごとにfprintf(leftfile, ); (x座標に気を付ける） } }
量子力学固有値問題計算プログラムの大枠； E を固定バージョン
#define error 1e‐8 #define V0 ‐0.001 #define Emin‐0.0010 #define Emax‐0.00085main(){ do{ ２分法でE の範囲を狭めていく 関数 diff を呼び出して波動関数計算 接続性が良くなるまで再トライ }while(エネルギーE での接続が良いか？); 接続性が合格したらその時の E を printf. その E を使って波動関数を最後に計算・ファイル 出力;関数plot } /* 左からと右からのψを計算し、接続 の良さを計算*/ double diff(double E){ /* 左からの計算 */ for(i=1; i<=1500; i++){ 波動関数の計算； three = *****; 関数k2の呼び出し } i = 1500 での値が決まる; ① /* 右からの計算 */ for(i=1; i<=500; i++){ 波動関数の計算; three = *****; 関数k2の呼び出し } i = 500 での値が決まる；② ①‐② の値を main へ返す } /* 波動関数計算中に k^2 を代入するのに分かりやすく ここで場合分けて計算する */
double k2(int i, double E){
if (i<500) return(なんたら); /* ポテンシャル外側*/ if (i>=500) return(なんたら); /* ポテンシャル内側*/ (V0 の符号に注意) } /* 正解のE で波動関数をファイル出力 */ void plot(double E){ rightfile = fopen(“right.dat”, “w”); leftfile = fopen(“left.dat”, “w”); /* 左からの計算 */ for(i=1; i<=1500; i++){ 波動関数の計算 1行ごとにfprintf(leftfile, ) ; (x座標に気を付ける） } /* 右からの計算 */ for(i=1; i<=500; i++){ 波動関数の計算 1行ごとにfprintf(leftfile, ); (x座標に気を付ける） } }
