数値計算法 011/6/8 林田清 大阪大学大学院理学研究科
常微分方程式の応用例 1 Rutherford 散乱 ( 原子核同士の散乱 ; 金の薄膜に α 粒子をあてる ) 1 クーロン力 f= 4 0 r r r Ze y からf cos, si f f f y f f 粒子の 方向 y方向の速度と座標について dv Ze dvy Ze y, 3 3 dt 40m r dt 40m r d dy v, vy dt dt を差分方程式に変換してt 0から計算していく 任意の位置 ( 衝突パラメータ ) に対する粒子の軌跡が描ける *) 引力に符号を変えればそのまま天体の衝突の式に使える
常微分方程式の応用例 恒星の内部構造 dp GM r 静水圧平衡 : dr r dlr 熱平衡 : 4 r dr dt 3 1 Lr 熱輸送の条件 : - 4 3 4 dr ac T r dm r 自己重力の条件 : 4r dr rを独立変数とした 4つの変数 ( 圧力 p, 温度 T, 半径 rより内側の質量 M, 半径 rの球面を通過するエネルギー流量 L r ) についての常微分方程式 r
常微分方程式の境界値問題 ( 例 ) 時間に依存しない Shrodiger方程式 0, 1に無限に高いポテンシャルの壁 d E m d (0) (1) 0 E' E/( ) として差分で近似すると m ( 1) ( ) ( 1) ( ) E ' ( ) ( ) ( ) 0 0 N 上の式は ( ( ),... ( )) に対する連立方程式 1 N 1
偏微分方程式 偏微分方程式 u(, t) u u 拡散方程式 t u u 波動方程式 t u u =0 ラプラス方程式 y 拡散方程式 熱の伝導時間に依存する Shrodiger 方程式 波動方程式 電磁波の伝播 弦の振動 ラプラス方程式 板の定常温度分布 Poisso 方程式
拡散方程式の解法 u u (0 1, t 0) t u(,0) ( ) (0 1) u(0, t) u(1, t) 0 ( t 0) 差分方程式で近似する, t t u 1 (, t ) { u (, t t ) u (, t )} t t t u 1 (, ) { (, ) (, ) (, )} t u t u t u t ( ) 1 1 ( U, 1 U, ) ( U 1, U, U 1, ) t ( ) t これから = ( ) U U (1 ) U として, 1 1,, 1, 初期条件 U 境界条件 U,0 ( ) U 0, J, 0 U 前の時刻 t での値を使って次の時刻 t 1 の値を計算できる陽公式 t U, t 0 0 J
拡散方程式の安定性条件 U f ( )ep( ik) で表されるような特解を考える, 差分方程式 U U (1 ) U U, 1 1,, 1, を使うと f ik ik f k f f (0) 1とすると ( 1) { ep( ) (1 ) ep( )} ( ) (1 4 si ) ( ) k k f ( ) (1 4si ), U, (1 4si ) ep( ik) 一般解はこの特解の重ね合わせ k k 1 発散しない条件は 1 4si 1つまりsi 任意のkに対して成り立つためには t ( ) 1
拡散方程式の解法 ( 陰公式 ) 1 1 ( U, 1 U, ) ( U 1, U, U 1, ) のかわりに t ( ) 1 1 ( U, 1 U, ) ( U 1, 1 U, 1 U 1, 1 ) t ( ) を考える U (1 ) U U U U 1, 1, 1 1, 1, U 0, J, 0 ( U,..., U ) から ( U,..., U ) を計算するためには 1, J 1, 1, 1 J 1, 1 連立一次方程式を解かなければならない (= 陰公式 ) t t U, t 0 0 J k 陰公式の場合特解はU, (1+ 4 si ) ep( ik) になるのでがいかなる値であっても発散しない つまり無条件安定
拡散方程式の差分漸化式を行列形式でかくと 陰公式では U (1 ) U U U 1, 1, 1 1, 1, 1 0 0 U U 1, 1 1, 1 U, 1 U, 0 0 1 UJ, 1 UJ, 0 0 1 U U J 1, 1 J 1, ちなみに陽公式では U U (1 ) U U, 1 1,, 1, U1, 1 1 0 0 U1, U, 1 1 U, 0 0 UJ, 1 1 UJ, U 0 0 1 U J 1, 1 J 1, ここで 方向の分点を 0,1,..., Jととっている
波動方程式の解法 u u (0 1, t 0) t u u(,0) ( ), (,0) ( ) (0 1) t u(0, t) u(1, t) 0 ( t 0) 差分方程式で近似する 1 1 ( U, 1 U, U, 1 ) ( U 1, U, U 1, ) ( t) ( ) t これから =( t/ ) として U U U ( U U U ), 1,, 1 1,, 1, これはt, t での値からt を決める式 初期条件 U,0-1 1 ( ) U U t ( ) ( U U U 境界条件 U U 0,1,0 1,0,0 1,0 0, J, 前の時刻 t, t での値を使って次の時刻 t の値を計算できる陽公式 -1 1 安定性の条件はt/ 1 ) t U, t 0 0 J
ラプラス方程式の解法 u u =0 (0 1,0 y 1) y u(,0) 1( ), u(,1) ( )(0 1) u(0, y) 1( y), u(1, y) ( y)(0 y 1) ih, y hとしてu( y ) に対する差分方程式の解をU とする i i, i, 1 1 ( U i1, Ui, Ui 1, ) ( U i, 1 Ui, Ui, 1) 0 h h これから 1 Ui, ( Ui 1, Ui 1, Ui, 1 Ui, 1 ) 4 これはU がその周囲の4 個の格子点での値の平均になっているという式 i, 境界条件は Ui,0 1 ( i ), Ui, N ( i ) ( i 0,1,..., N) U0, 1 ( y ), U N, ( y ) ( 0,1,..., N) 上の差分式は ( U, U,..., U, U,... U ) に対する連立一次方程式になっている 1,1,1 N,1 1, N 1, N
連立 1 次方程式の解法 連立一次方程式 a1,1 a1,... a1, N 1 y1 a,1 a,... a, N y............. an,1 an,... a N, N N y N A y Aは係数行列 解が一意に存在する必要十分条件は行列式 det( A) の値が0でないこと
常微分方程式の応用例 1 Rutherford 散乱 ( 原子核同士の散乱 ; 金の薄膜に α 粒子をあてる ) 1 クーロン力 f= 4 0 r r r Ze y からf cos, si f f f y f f 粒子の 方向 y方向の速度と座標について dv Ze dvy Ze y, 3 3 dt 40m r dt 40m r d dy v, vy dt dt を差分方程式に変換してt 0から計算していく 任意の位置 ( 衝突パラメータ ) に対する粒子の軌跡が描ける *) 引力に符号を変えればそのまま天体の衝突の式に使える
レポート問題 回目 ( 締め切り 7/6) 金 (Z=79) の原子核に α 粒子を衝突させるラザフォード散乱の軌跡を図に描いてみよう 基礎になるのは クーロン力による運動方程式 ( 常微分方程式の応用例 1) 金の原子核 1 個を y 座標の原点におき 方向 y 方向とも -000fm から +000fm の範囲 (fm は 10-15 m) での α 粒子の運動を考える 連立常微分方程式を数値的にとくことで軌跡を求める できれば ルンゲクッタ法を用いるのが望ましいが オイラー法でもよい 衝突パラメータとしては 0fm( 正面衝突 ) を含めて 5 個から 10 個程度の値を自分で設定し 1 枚のグラフにそれらの軌跡をかさねてプロットせよ プロットには guplot を使用することを想定しているが それ以外でもよい 入射 α 粒子は - 方向から 軸に平行に入射するとする =-000fm での運動エネルギーは 5.3MeV とし 相対論的効果は無視してよい 物理定数 粒子の質量など必要なら自分で調べること 提出は 011/7/6 までに khclass@ess.sci.osaka-u.ac.p までメールすること メイルのタイトルは report_ 学籍番号とすること メールの本文には 自分で作成したソースプログラムとともに 設定した衝突パラメータの値も記すこと 加えて 軌跡をプロットした結果 気がついたことがあれば一言かきそえよ プロットした結果は eps ファイル形式等で保存し メール添付ファイルとしていっしょに提出すること 以上 授業でも説明します
while の利用 ( 条件判断とループ ) eps=1.0e-5; s=10; sold=100; while( fabs(sold-s)>eps) { } sold=s; s=sold*0.1; sew=1; do{ s=sew; sew=s*0.1; }while(fabs(sew-s)>eps); C 言語では goto 文は推奨されない C3456 eps=1.0e-5 s=10 sold=100 10 cotiue if( dabs(sold-s).le.eps ) * goto 0 sold=s s=sold*0.1 goto 10 0 cotiue sew=1 110 cotiue s=sew sew=s*0.1 if( dabs(sew-s).gt.eps ) * goto 110 10 cotiue
関数 サブルーチンの利用 (Fortra) c34567 real*8 a,b,c, fuc a=.0 b=3.0 c=fuc(a,b) write(*,*) a,b,c stop ed c34567 real*8 a,b,c, fuc a=.0 b=3.0 call sub(a,b,c) write(*,*) a,b,c stop ed fuctio fuc(a,b) real*8 a,b fuc=a+b; retur ed subroutie sub(a,b,c) real*8 a,b,c c=a+b; retur ed
関数の利用 (C) double hayasida(double a,double b); mai() { double a,b,c; a=.0; b=3.0; c=hayasida(a,b); pritf(" %lf %lf %lf ",a,b,c); } double hayasida( double a, double b) {double y; y=a+b; retur(y) ; } double fuc(double a,double b, double *y); mai() { double a,b,c; a=.0; b=3.0; fuc(a,b,&c); pritf(" %lf %lf %lf ",a,b,c); } double fuc( double a, double b, double *y) { *y=a+b; } ポインターについては本で学習してください