数値計算法
2011/6/8
林田 清
常微分方程式の応用例１
Rutherford散乱（原子核同士の散乱；金の薄膜にα粒子
をあてる）
2 2 2 2 3 32
cos
,
sin
y
2
2
,
,
0
x y y x x yZe
x
y
f
f
f
f
f
f
r
r
r
dv
dv
Ze
x
Ze
y
dt
m r
dt
m r
dx
dy
v
v
dt
dt
t
0 0 01
クーロン力f=
から
4
粒子のx方向、 方向の速度と座標について
4
4
を差分方程式に変換して
から計算していく。
任意の位置（衝突パラメータ）に対する粒子の軌跡が描ける。
＊）引力に符号を変えればそのまま天体の衝突の式に使える
常微分方程式の応用例２
恒星の内部構造
2 2 3 2 24
3
1
-4
4
4
,
,
)
r r r r r rGM
dp
dr
r
dL
r
dr
L
dT
dr
ac T
r
dM
r
dr
r
p
T
r
M
r
L
静水圧平衡：
熱平衡：
熱輸送の条件：
自己重力の条件：
を独立変数とした、４つの変数(圧力 ,温度 半径 より内側の質量
半径 の球面を通過するエネルギー流量
についての常微分方程式
常微分方程式の境界値問題
2 2 20,
1
2
(0)
(1)
0
Shrodinger
x
x
d
E
m dx
（例）時間に依存しない
方程式
に無限に高いポテンシャルの壁
2 2 1 1 0 1 1
'
/(
2
(
) 2 ( )
(
)
(
)
' ( )
( )
(
)
0
( ),... (
))
j j j j N NE
E
m
x
x
x
x E
x
x
x
x
x
)として差分で近似すると
上の式は(
に対する連立方程式
偏微分方程式
拡散方程式
熱の伝導
時間に依存する
Shrodinger方程式
波動方程式
電磁波の伝播、弦の振動
ラプラス方程式
板の定常温度分布
Poisson方程式
2 2 2 2 2 2 2 2 2 2( , )
u x t
u
u
t
x
u
u
t
x
u
u
x
y
偏微分方程式
拡散方程式
波動方程式
=0 ラプラス方程式
拡散方程式の解法
x 2 2 2 , 1 , 2 1, , 1, 2 , 1 1, , , 1 ( , ) { ( , ) ( , )} 1 ( , ) { ( , ) 2 ( , ) ( , )} ( ) 1 1 ( ) ( 2 ) ( ) ( ) (1 2 ) j n j n j n j n j n j n j n j n j n j n j n j n j n j n j n j x j x t n t u x t u x t t u x t t t u x t u x x t u x t u x x t x x U U U U U t x t x U U U 差分方程式で近似する これから = として 1, ,0 0, , 1 ( ) 0 n j n j j n J n n n U U x U U t t 初期条件 境界条件 前の時刻 での値を使って次の時刻 の値を計算できる陽公式 2 2 (0 1, 0) ( , 0) ( ) (0 1) (0, ) (1, ) 0 ( 0) u u x t t x u x x x u t u t t x
t
0t
0x
x
j nt
U
j n, Jx
拡散方程式の安定性条件
, , 1 1, , 1, 2 2 2 ,( ) exp(
)
(1 2 )
(
1)
{ exp(
) (1 2 )
exp(
)} ( )
(1 4 sin
) ( )
2
(0) 1
( )
(1 4 sin
) ,
(1 4 sin
) exp(
)
2
2
j n j n j n j n j n n n j nU
f n
ikj x
U
U
U
U
k x
f n
ik x
ik x
f n
f n
f
k x
k x
f n
U
ikj x
で表されるような特解を考える
差分方程式
を使うと
とすると
一般解
2 2 21
1 4 sin
1
sin
2
2
2
1
2
(
)
k x
k x
t
k
x
はこの特解の重ね合わせ
発散しない条件は
つまり
任意の に対して成り立つためには
拡散方程式の解法２（陰公式）
, 1 , 2 1, , 1, , 1 , 2 1, 1 , 1 1, 1 1, 1 , 1 1, 1 , 0, , 1, 1, 1, 1 1, 1 1 1 ( ) ( 2 ) ( ) 1 1 ( ) ( 2 ) ( ) (1 2 ) 0 ( ,..., ,..., ) j n j n j n j n j n j n j n j n j n j n j n j n j n j n n J n n J n n J n U U U U U t x U U U U U t x U U U U U U U U U U のかわりに を考える。 )から( を計算するためには 連立一次方程式を解かなければならない（＝陰公式） 2 , (1 4 sin ) exp( ) 2 n j n k x U ikj x 陰公式の場合特解は ＋ になるので がいかなる値であっても発散しない、つまり無条件安定x
t
0t
0x
x
j nt
U
j n, Jx
拡散方程式の差分漸化式を行列形式でかくと。。。
1, 1 , 1 1, 1 , 1, 1 1, 2, 1 2, 2, 1 2, 1, 1 1, (1 2 ) 1 2 0 0 1 2 0 0 1 2 0 0 1 2 j n j n j n j n n n n n J n J n J n J n U U U U U U U U U U U U 陰公式では , 1 1, , 1, 1, 1 1, 2, 1 2, 2, 1 2, 1, 1 1, (1 2 ) 1 2 0 0 1 2 0 0 1 2 0 0 1 2 j n j n j n j n n n n n J n J n J n J n U U U U U U U U U U U U ちなみに陽公式では 0,1,..., j x x j x j J ここで 方向の分点を ととっている波動方程式の解法
, 1 , , 1 1, , 1, 2 2 2 , 1 , , 1 1, , 1, -1 1 ,0 ,1 ,0 1,0 ,0 1 1 ( 2 ) ( 2 ) ( ) ( ) / ) 2 ( 2 ) , ( ) ( ) ( 2 2 j n j n j n j n j n j n j n j n j n j n j n j n n n n j j j j j j j j U U U U U U t x t x U U U U U U t t t U x U U t x U U U 差分方程式で近似する これから =( として これは での値から を決める式 初期条件 1,0 0, , -1 1 ) 0 , / n J n n n n U U t t t t x 境界条件 前の時刻 での値を使って次の時刻 の値を計算できる陽公式 安定性の条件は １ 2 2 2 2 (0 1, 0) ( , 0) ( ), ( , 0) ( ) (0 1) (0, ) (1, ) 0 ( 0) u u x t t x u u x x x x x t u t u t t x
t
0t
0x
x
j nt
U
j n, Jx
ラプラス方程式の解法
, , 1, , 1, , 1 , , 1 2 2 , 1, 1, , 1 , 1 , ,0 1 , 2 , ( ) 1 1 ( 2 ) ( 2 ) 0 1 ( ) 4 ( ), ( ) ( i j i j i j i j i j i j i j i j i j i j i j i j i j i j i j i i i N i x ih y jh u x y U U U U U U U h h U U U U U U U x U x i として に対する差分方程式の解を とする これから これは がその周囲の４個の格子点での値の平均になっているという式 境界条件は 0, 1 , 2 1,1 2,1 ,1 1,2 1, 0,1,..., ) ( ), ( ) ( 0,1,..., ) , ,...., , ,.... ) j j N j j N N N N U y U y j N U U U U U 上の差分式は( に対する連立一次方程式になっている 2 2 2 2 1 2 1 2 (0 1, 0 1) ( , 0) ( ), ( ,1) ( )(0 1) (0, ) ( ), (1, ) ( )(0 1) u u x y x y u x x u x x x u y y u y y y =0連立１次方程式の解法
1,1 1,2 1, 1 1 2,1 2,2 2, 2 2 ,1 ,2 ,...
...
..
...
...
...
.
.
...
N N N N N N N Na
a
a
x
y
a
a
a
x
y
a
a
a
x
y
Ax
y
A
連立一次方程式
は係数行列
det( )
A
解が一意に存在する必要十分条件
は行列式
の値が０でないこと
常微分方程式の応用例１
Rutherford散乱（原子核同士の散乱；金の薄膜にα粒子
をあてる）
2 2 2 2 3 32
cos
,
sin
y
2
2
,
,
0
x y y x x yZe
x
y
f
f
f
f
f
f
r
r
r
dv
dv
Ze
x
Ze
y
dt
m r
dt
m r
dx
dy
v
v
dt
dt
t
0 0 01
クーロン力f=
から
4
粒子のx方向、 方向の速度と座標について
4
4
を差分方程式に変換して
から計算していく。
任意の位置（衝突パラメータ）に対する粒子の軌跡が描ける。
＊）引力に符号を変えればそのまま天体の衝突の式に使える
レポート問題２回目（締め切り７
/６）
金（Z=79)の原子核にα粒子を衝突させるラザフォード散乱の軌跡を図に描 いてみよう。 基礎になるのは、クーロン力による運動方程式（常微分方程式の応用例１）。 金の原子核１個をxy座標の原点におき、x方向、y方向とも-2000fmから +2000fmの範囲(fmは10-15m)でのα粒子の運動を考える。 連立常微分方程式を数値的にとくことで軌跡を求める。できれば、ルンゲクッ タ法を用いるのが望ましいが、オイラー法でもよい。 衝突パラメータとしては、0fm(正面衝突）を含めて、5個から10個程度の値を 自分で設定し、１枚のグラフにそれらの軌跡をかさねてプロットせよ。プロット にはgnuplotを使用することを想定しているが、それ以外でもよい。 入射α粒子は-x方向からx軸に平行に入射するとする。x=-2000fmでの運動 エネルギーは5.3MeVとし、相対論的効果は無視してよい。 物理定数、粒子の質量など必要なら自分で調べること。 提出は2011/7/6までに、khclass@ess.sci.osaka-u.ac.jp までメールするこ と。メイルのタイトルは report2_学籍番号とすること。メールの本文には、自 分で作成したソースプログラムとともに、設定した衝突パラメータの値も記す こと。加えて、軌跡をプロットした結果、気がついたことがあれば一言かきそ えよ。 プロットした結果はepsファイル形式等で保存し、メール添付ファイルと していっしょに提出すること。 以上、授業でも説明します。whileの利用 (条件判断とループ）
eps=1.0e-5; s=10; sold=100; while( fabs(sold-s)>eps) { sold=s; s=sold*0.1; } snew=1; do{ s=snew; snew=s*0.1; }while(fabs(snew-s)>eps); C23456 eps=1.0e-5 s=10 sold=100 10 continue if( dabs(sold-s).le.eps ) * goto 20 sold=s s=sold*0.1 goto 10 20 continue snew=1 110 continue s=snew snew=s*0.1 if( dabs(snew-s).gt.eps ) * goto 110 120 continueC言語ではgoto文は推奨されない。
関数、サブルーチンの利用
(Fortran)
c234567
real*8 a,b,c, func a=2.0 b=3.0 c=func(a,b) write(*,*) a,b,c stop end function func(a,b) real*8 a,b func=a+b; return end c234567
real*8 a,b,c, func a=2.0 b=3.0 call sub(a,b,c) write(*,*) a,b,c stop end subroutine sub(a,b,c) real*8 a,b,c c=a+b; return end
関数の利用
（C）
double hayasida(double a,doubleb); main() { double a,b,c; a=2.0; b=3.0; c=hayasida(a,b); printf(" %lf %lf %lf¥n",a,b,c); } double hayasida( double a, double b) {double y; y=a+b; return(y) ; }
double func(double a,double b, double *y); main() { double a,b,c; a=2.0; b=3.0; func(a,b,&c); printf(" %lf %lf %lf¥n",a,b,c); } double func( double a, double b, double *y) { *y=a+b; } ポインターについては本で学習してください
