計算物理学２第 13 回：微分方程式の解法その３

計算物理学２第 13 回：微分方程式の解法その３

ver. 2018/7/20

偏微分方程式

(partial differential equation)

一般的に二階二変数の偏微分方程式は

A(x, y) ∂ ² U

∂x ² + 2B(x, y) ∂ ² U

∂x∂y + C(x, y) ∂ ² U

∂y ² + D(x, y) ∂U

∂x + E(x, y) ∂U

∂y + F (x, y) = 0 (1)

B ² (x, y) − A(x, y)C(x, y)

3

すべての

x, y

B ² (x, y) − A(x, y)C(x, y) < 0

(elliptic partial differential equation)

Poisson

∂ ² U

∂x ² + ∂ ² U

∂y ² = − 4πρ(x, y) (2)

があります。

Poisson

U(x, y)

B ² (x, y) − A(x, y)C(x, y) > 0

(hyperbolic partial differential equation)

1

∂ ² U

∂t ² − c ² ∂ ² U

∂x ² = 0 (3)

が該当します。

B ² (x, y) − A(x, y )C(x, y) = 0

(parabolic partial differential equation)

1

∂U

∂t = λ C V

∂ ² U

∂x ² (4)

が該当します。

1 ^{楕円型偏微分方程式}

Poisson

法によって解くことができます。

例えば三点公式を用いると

∂ ² U

∂x ² (x i , y j ) = U (x i+1 , y j ) − 2U(x i , y j ) + U (x i − 1 , y j )

(∆x) ² , (5)

∂ ² U

∂y ² (x i , y j ) = U (x i , y j+1 ) − 2U(x i , y j ) + U (x i , y j − 1 )

(∆y) ² , (6)

などのように

U(x i , y j ), (i = 1, · · · , N x , j = 1, · · · , N y )

これら

N x N y

U

Poisson

AU = b (7)

の形の連立一次方程式となります。

b

ρ(x i , y j )

Jacobi, Gauss-Seidel, SOR,

2 ^{双曲型偏微分方程式}

次回時間があれば取り扱う予定です。

3 ^{放物型偏微分方程式}

放物型偏微分方程式として

1

Schr¨ odinger

i¯ h ∂

∂t ψ(x, t) = [

− ¯ h ² 2m

∂ ²

∂x ² + V (x) ]

ψ(x, t) (8)

の形で、初期値としては時刻

t = 0

ψ(x, t = 0)

例えば初期値としてはガウス波束

ψ(x, 0) = Ce ⁻

2

a2

e ^ik

^x (9)

を考えます。係数

C

a

x 0

p 0 = ¯ hk 0

空間

0 < x < L

(

ます

)

N

x i = − L + i∆x, (

∆x = 2L

N + 1 , i = 0, 1, 2, · · · , N + 1 )

(10)

とします。時間についても適当なタイムステップ

∆t

t j = j∆t, (j = 0, 1, 2, · · · ) (11)

とします。

波動関数の時間発展は形式的に

ψ(x, t) = exp [

− i

¯ h

Ht ˆ ]

ψ(x, t = 0) (12)

とかけます。ここで

Hamiltonian

H ˆ = − ¯ h ² 2m

∂ ²

∂x ² + V (x) (13)

で与えられます。

3.1

式

(12)

t = t j

t j+1

∆t

Taylor

1

ψ(x i , t j+1 ) = (

1 − i

¯ h

H∆t ˆ )

ψ(x i , t j ) + O ((∆t) ² ) (14)

とかけます。これは、時間依存

Schr¨ odinger

i¯ h ψ(x _i , t _j+1 ) − ψ(x _i , t _j )

∆t + O (∆t) = − ¯ h ² 2m

ψ(x _i+1 , t _j ) − 2ψ(x _i , t _j ) + ψ(x _{i − 1} , t _j )

(∆x) ² + V (x _i )ψ(x _i , t _j ) + O ((∆x) ² ), (15)

t = t _j

t = t _j+1

ψ(x i , t j+1 ) = ψ(x i , t j ) − i ∆t

¯ h

[

− ¯ h ² 2m

ψ(x i+1 , t j ) − 2ψ(x i , t j ) + ψ(x i − 1 , t j )

(∆x) ² + V (x i )ψ(x i , t j ) ]

+ O ((∆t) ² , ∆t(∆x) ² ) (16)

これは陽解法

(explicit method)

Euler

Runge-Kutta

(

FTCS(Forward-Time Central-Space)

)

(16)

i¯ h ψ(x i , t j+1 ) − ψ(x i , t j − 1 )

2∆t + O ((∆t) ² ) = − ¯ h ² 2m

ψ(x i+1 , t j ) − 2ψ(x i , t j ) + ψ(x i − 1 , t j )

(∆x) ² + V (x i )ψ(x i , t j ) + O ((∆x) ² ), (17)

ψ(x i , t j+1 ) =ψ(x i , t j − 1 ) − i 2∆t

¯ h

[

− ¯ h ² 2m

ψ(x _i+1 , t _j ) − 2ψ(x _i , t _j ) + ψ(x _{i − 1} , t _j )

(∆x) ² + V (x i )ψ(x i , t j ) ]

+ O ((∆t) ³ , ∆t(∆x) ² ) (18)

この式では

t = t j+1

t = t j

t j − 1

t = ∆t

t j − 1

(16)

(18)

時間依存

Schr¨ odinger

ことです。

t = 0

∫ L

− L

dx | ψ(x, t) | ² = 1 (19)

は満たされるように規格化定数

C

∆x

∆t

3.2

陰解法

(implicit method)

ψ(x i , t j+1 ) = (

1 + i

¯ h

H ˆ ∆t ) ₋ 1

ψ(x i , t j ) (20)

となり、時間に関して後退差分を用います。

i¯ h ψ(x i , t j+1 ) − ψ(x i , t j )

∆t = − ¯ h ² 2m

ψ(x i+1 , t j+1 ) − 2ψ(x i , t j+1 ) + ψ(x i − 1 , t j+1 )

(∆x) ² + V (x i )ψ(x i , t j+1 ) (21)

この式を解く際に既知の波動関数は

ψ(x i , t j )

¯ h ²

2m(∆x) ² ψ(x _i+1 , t _j+1 ) + (

− ¯ h ²

m(∆x) ² + V (x _i ) + i ¯ h

∆t )

ψ(x _i , t _j+1 ) + ¯ h ²

2m(∆x) ² ψ(x _{i − 1} , t _j+1 ) = i ¯ h

∆t ψ(x _i , t _j ) (22)

ψ(x i , t j+1 ) (i = 1, 2, · · · N )

t j

3.3 Crank-Nicolson

Crank-Nicolson

exp ( i

¯ h ∆t H ˆ

)

= 1 − 1

2 i

¯ h

H∆t ˆ

1 + 1 2 i

¯ h

H∆t ˆ

(23)

と近似します。これを用いると

t _j

t _j+1

(

1 + 1 2 i

¯ h

H ˆ ∆t )

ψ(x _i , t _j+1 ) = (

1 − 1 2 i

¯ h

H ˆ ∆t )

ψ(x _i , t _j ) (24)

となります。

Crank-Nicolson

t j

t j+1

Hamiltonian

t j , t j+1

i¯ h ψ(x i , t j+1 ) − ψ(x i , t j )

∆t = 1

2 [

− ¯ h ² 2m

ψ(x i+1 , t j+1 ) − 2ψ(x i , t j+1 ) + ψ(x i − 1 , t j+1 )

(∆x) ² + V (x _i )ψ(x _i , t _j+1 ) ]

+ 1 2

[

− ¯ h ² 2m

ψ(x i+1 , t j ) − 2ψ(x i , t j ) + ψ(x i − 1 , t j )

(∆x) ² + V (x i )ψ(x i , t j ) ]

(25)

この式も陰解法のときと同様に連立方程式の形にして時刻

t j

(Crank-Nicolson

) Crank-Nicolson

4 ^可視化

時間依存する波動関数

ψ(x, t)

(46)

(46)

mkdir

4.1

まず時刻

t

1

(46)

wavefunction 0000000.dat

(0

) x

Reψ(x, t), Imψ(x, t)

3

時刻ごとに出力するデータファイルの名前を変更するにはプログラムの中で少しテクニックが必要です。

Fortran

ti

7

(I7.7)

tindex

WRITE

文で変換しています。これを用いてファイル名を文字列の連結で作成しています。

Fortran

//

(

50

)

TRIM

C

言語では

sprintf

4.2 gnuplot

データファイルのセットが出来上がったら

gnuplot

gif

plot.gpl

data.gpl

plot.gpl

set term gif animate

とします。出力ファイル名は拡張子

”.gif”

x

y

set xrange

set yrange

tmin, tmax, dt(

)

data.gpl

data.gpl

t

tmin

定します。データファイルの内容をプロットし、

t

dt

tmax

reread

data.gpl

tmin

tmax

ターミナルで

% gnuplot plot.gpl

とすると

gif

gif

ファイルマネージャーでディレクトリを開き、ダブルクリックしてください。

4.3 linux

多くのデータファイルを作成しますのでディレクトリの中に大量のファイルが作成されます。容量等にも注 意を払ってください。

% ls -lh

と

-h

% du -h .

とします。

”.”

複数のファイルをまとめて操作する場合はワイルドカード

*

*

% ls -lh wavefunction_00003*.dat

とすると、時刻が

0000300

0000399

% rm wavefunction_*.dat

とすると、

wavefunction

.dat

% rm *

とするとディレクトリにある全てのファイルが削除されますので気をつけてください。特に

wavefunc- tion *.dat

wavefunction * .dat

5 演習問題

(46) ¯ h ² /2m = 20 (MeV)

∆t = 10 ^{− 4} (¯ h/MeV),

L = 100 (fm),

N = 1000

t = 0

(a = 5 fm, x ₀ = 20 (fm), k = 5 (fm ^{− 1} ))

Schr¨ odinger

t = 5000(¯ h/MeV)

(47) x = L/2

b (fm)

V ₀ (MeV)

b

V 0

(48)

(49)

ψ(x + L, t) = ψ(x, t)

(50) Crank-Nicolson