偏微分方程式、連立１次方程式、乱数

Transcription

1 数値計算法 2011/7/6 林田清 大阪大学大学院理学研究科

2 モジュール化 長い長いメインプログラムは間違いのもと また 問題を少し変えるたびにデバッグが必要 プログラムを複数の部品 ( モジュール ) にわけて作成することで プログラムの見通しをよくする デバッグの労力を減らせる 複数の人で手分けして作業できる 一部の部品は別のプログラムでも利用できる C 言語では function を利用することでモジュール化できる プログラムに必要な機能を抽象化して モジュール間のインターフェスを明確に定義するのが重要 Functionの引数があまりにも多いようなら インターフェースの切り分けが適当でないと考えるべき モジュール化をより発展させた オブジェクト指向 プログラミング C++ は代表的な オブジェクト指向 言語抽象化 カプセル化 継承 多態 classの概念 private publicの区別など

3 関数の利用 (C) double hayasida(double a,double b); 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; } ポインターについては本で学習してください

4 2 dvx Z Z Aue x 3 dt 4 0m r 2 dvy Z Z Aue y 3 dt 4 0m r dx dt dy dt v v x y v, v, x, yの 4変数の連立一次微分方程式 x y 一般的にn 変数の連立一次常微分方程式をルンゲクッタ法でとく (1stepすすめる) 関数例 #define NMAX 10 void nextstep( int n, double t, double *f, double (*dfdt)(int, double, double *), double dt) { int i; double ft[nmax]; double k1[nmax],k2[nmax],k3[nmax],k4[nmax]; for(i=0;i<n;i++) ft[i]=f[i]; for(i=0;i<n;i++) k1[i]=dfdt(i, t, ft); for(i=0;i<n;i++) ft[i]=f[i]+0.5*k1[i]*dt; for(i=0;i<n;i++) k2[i]=dfdt(i, t+0.5*dt, ft);. }

5 構造体とクラス 構造体 (C 言語でも利用可能 ) struct Particle { double m, q, p, a; }; Particle p1, p2; p1.m = 1.0; p1.q += p1.p/p1.m * 0.001; クラス (C++ で利用可能 ) class Particle { double m, q, p, a; public: void Evolve( double dt ){ q += p/m*dt; p += a*dt; } }; Particle p1, p2; p1.m = 1.0; p1.q += p1.p/p1.m * 0.001; p1.evolve( ); 参照

6 電気力線の書き方 ポアソンの方程式を ( 数値的にor 解析的に ) とい 2 てポテンシャルを求める ポテンシャルの等しい面 ( 線 ) をかく 等ポテンシャル線に常に垂直になるように線を引く タートル法 電場を計算し 電場の方向にそって Turtle をすすめる Q x x Q x x E E x y r1 r2 Q y y Q y y r1 r2 y Q1 Q2 x

7 ポアッソン分布に従う乱数 P x e x x - ( ) /! ただしxは離散値 棄却法はそのままでは使えない Q y dy e y [ y] - ( ) /[ ]! という連続分布関数を定義する ここで [y] はyを超えない整数 Q( y) に従う乱数 yを棄却法により発生させ その整数部をとってxとする

8 シミユレーション (Simulation) = 模擬実験 物理における計算機シミユレーション 現実の ( 複雑な ) 物理現象に対して 物理法則に基づいたモデルをつくり再現を試みる 仮定した物理法則あるいはモデルの検証をする 現象の予想 パラメータの最適化を行う ( 厳密解の計算が困難な場合に ) 近似解を求める モンテカルロ法 ( 乱数を用いた統計的手法 ) を利用したシミュレーション=モンテカルロシミュレーション 微分方程式を解く あるいは連立方程式を解く問題に帰着するシミレーション問題も多い 今回の参考文献 : シミュレーション物理入門 矢部 川田 福田著 ( 朝倉書店 )

9 物理におけるシミュレーション具体例 1 気体 液体の運動 流体力学の方程式を解く方法 個々の粒子の運動を追う方法 ( 分子動力学法 ) 重力多体系 ( 恒星の集団としての銀河等 ) の計算 荷電粒子の運動も同じ扱いが使える このようなシミュレーションでは 通常 空間を格子 ( メッシュ ) 上に分割しそれぞれの場所でのポテンシャル 力等を計算していく 注 ) 初期値を与えるのに乱数を使うことも多いが 乱数を用いた確率論的方法が計算の柱になるわけではない モンテカルロシミュレーションではない

10 物理におけるシミュレーション具体例 2 物質中での放射線の相互作用 検出器の応答のシミュレーション イジングモデル 磁性体のモデル 量子モンテカルロ法 拡散過程の応用 これらは モンテカルロシミュレーション 注 ) モンテカルロ法はシミュレーションだけに使われるわけではない 例 ) モンテカルロ積分

11 分子動力学法 1 液体のように膨大な数の分子からなる物質のマクロな性質を決める方法 ( 化学反応の再現等にも利用される ) 基本的には 各粒子の位置 速度の初期値を与え 下の運動方程式を差分化した式に従って時間発展させていけばよい 2 粒子間の相互作用のポテンシャルUr () とする du rij j粒子から i粒子に働く力はf, またF F dr r dv i粒子の運動方程式は dt i ij ji ij j F ij, dx dt i ij v i

12 分子動力学法 2 現実的には 計算できる粒子の数は限られている そのために用いられる便法 空間を小さなセルに分割し 周期的境界条件を仮定する近接ポテンシャルと集団ポテンシャルに分けて取り扱う 2 N Ui U( rij ) U( r)4 r 2V r c j(r <r ) ij c 2 dr r c

13 分子動力学のシミュレーション例 Farid F. Abraham, D. Brodbeck, R.A. Rafey and W.E. Rudge, Phys. Rev. Lett. 73, 272 (1994). 張力下での 2 次元固体のひび割れを 100 万個の原子を使ってシミュレーションした シミュレーションにより 観測困難な現象や 実験困難な条件下の現象が理解できる シミュレーション中には 実験で発見された 沢山の現象が見られた より引用 アニメーションは左の web ページで直接みてください

14 重力多体系の計算 粒子 ( 恒星 ) 間に働く重力を計算し 各粒子の運動方程式を差分的に解けば原理的には粒子の配置 運動を追うことができる ポアッソン方程式を解く問題に帰着させることもできる 重力ポテンシャルは Poissonの方程式 4 Gを満たす (c.f. クーロン力静電ポテンシャルは / を満たす ) 1) 空間にメッシュを切り粒子を配置する 粒子の速度の初期値も与える 2) 各メッシュ点でのを計算する 3) 差分化したPoissoの方程式 ( 偏微分方程式 ) を解き 各メッシュ点でのを求める 4) 各粒子に働く力をF -mにより求める i i dvi dxi 5) Fi / mi, v i を差分化した式を使い 各粒子の位置 速度を更新する dt dt 6)2) に戻り操作を繰り返す e 0

15 銀河の衝突 月の起源 Grape Project 重力多体系の2 体相互作用部分を計算する専用計算機を開発 より 例 )Two Disk Galaxies 製作 P.Frisch 例 ) 月の起源 製作 Kokubo et al. kit/index.html に解説あり 参照のこと

16 放射線の物質中での振る舞い 放射線粒子が物質に入射したときの振る舞いをモンテカルロシミュレーションする 放射線粒子が物質中の原子と反応する位置を決定する 考えている反応の平均自由行程を L (L=1/ns とすると p(x)dx=exp(-x/l)/ldx という確率分布に従う乱数 x を発生させる 反応によって粒子が生じる場合 その放出方向 エネルギー等を決定する 微分断面積の式 dwq,に従って放出方向を決定するには dwq,という確率分布に従う乱数 q,を発生させる 特に検出器中での振る舞いをシミュレートする場合は 検出器の信号発生も模擬し直接実験データと比較できる出力を作成する

17 平均自由行程 断面積 到達 ( 衝突 ) 距離の分布 例として光が物質中で吸収される過程を考える 光が物質に入射し微小な厚みdxを通過する前の強度を I 通過したあとの強度 をI diとする diはiとdx に比例することが期待されるので ある定数 kを用い てdI kidxとかける これを積分すると 厚みxを通過したあとの光の強度 I( x) は I( x) I(0)exp( kx) とかける kを吸収係数 l 1/ kを平均自由行程とよぶ 吸収にかかわる物質中の粒子の 個数密度を n それぞれの断面積をsとするとns kの関係がある ( 証明してみよ ) I I di dx 個々の光子が ( 十分厚い ) 物質中の粒子と衝突して吸収されるまでの距離 xの分布は 1 p( x) dx k exp(- kx) dx exp(- x / l) dx l でかける *) ただし (0 x ) で xが有限の範囲に限定される場合は規格化ファクタが異なる 以上は一般の粒子の衝突などにも拡張できる

18 等方散乱のシミュレーション 光が物質中で粒子に衝突し等方に散乱される過程を考える (1) 平均自由行程をとすると 衝突するまでの距離 l xの分布は 1 p( x) dx exp( x / l) l 0,1 の範囲の一様乱数 rを発生させて 逆変換法でx -l log( r ) により求められる 1 1 (2)3次元極座標での角度を極角 q 0,, 方位角 0,2 で指定することにすると 1 1 等方分布は p( q ) dq sin qdq p( ) d dで表すことができる 2 2 それぞれ定義範囲で積分すると1になるように規格化されている 2 a) 方位角 0,2 0,1 の範囲の一様乱数 rを発生させて 2rとして求められる 2 b) 極角 q 0, cos qという変数変換をおこなうとp( q ) dq sinqdq p dから 2 p 1 dq 1 sin 2 q d 2 0,1 の範囲の一様乱数 rを発生させて 2-1 cos として求められる *) 1 3 r3 q q 逆変換法に従うと 累積確率 Pq p( q) dq cosq 1 cosq r r r 1 cos q として求められる と -1 は同じ分布をするので結果は同じ q 0

19 素粒子実験での応用例 コンピューターシミュレーションによる BELLE 検出器で観測される B 中間子の崩壊事象 中央飛跡測定器の電子信号 ( 印 ) から発生した粒子の飛跡を求め 運動量を割り出します ( openhouse96/b_bell. html を転載 )

20 イジングモデル 1 磁性体を記述するもっとも簡単なモデル 相転移も再現できる 格子状に粒子を配置し そのスピンが上向き (1) か下向き (-1) しかないとする スピンは隣あう格子上のスピンとしか相互作用しない H - JSiS j i, j <i,j> は隣あう格子を意味する Jが正なら強磁性 Jが負なら反強磁性 格子点のスピンの組み合わせで指定される状態の確率分布をモンテカルロ法によりシミュレートする

21 イジングモデル 2 1) 各格子点にスピン (+1か-1) をランダムに与える ( これを状態 m1とする ) 2) ある 1個の格子点のスピンを逆転させる ( これを状態 m2とする ) 3) 状態 m1 状態 m2の実現確率の比 rを計算する Hm2 Hm 1 r fm2 / fm 1 exp( ) kt 4) 0から1までの一様乱数 sを発生させ s rのときこのスピンを逆転させる ( 状態 m2を採用する ) そうでないときは状態 m1のまま 5)4) で採用された状態を新たな状態 m1として 2) へ戻り 次のステップを行う 十分平衡状態に達するように全格子点のスピンについて何度か 2)-5) の操作 を繰り返す そのあと 各ステップの 4) で採択された状態は 状態の実現確率 を再現したものになる 参考 ) イジングモデルのアプレット

22 量子力学の固有値問題と変分法 2 2 d H x x V x x E x 2 2m dx 0 1 固有値を低い方からE, E,... と並べ 対応する固有関数を u, u,... とする Hu n n n x try n n n0 x 2 2 try n n n n0 n0 try 0 E u a u 0 1 ここで規格化された試行関数 E a E a E E 0 0 try 試行関数が規格化されていない場合はどういう表式になるか H=H 0 +V とおいた摂動法との関係は? x をパラメータを含む形で仮定し E なるようにパラメータを決める このときのE がE に対する近似解 を考える try が最小に try の最小値

23 量子モンテカルロ法 1 時間依存のシュレディンガー方程式は虚時間の拡散方程式に帰着できる 拡散方程式をモンテカルロ方で解くことで エネルギー固有値決定する 2 2 V -i t 2m で時間を it / 2 2 V 2m これは拡散定数 D= に変換すれば 2 /2mの拡散方程式と同じ形 =exp( iet / ) ( x) の形にかけるとき=exp( E ) ( x) つまりEは虚時間拡散方程式の減衰率として求まる

24 量子モンテカルロ法 2( 箱型ポテンシャルへの応用例 ) 0(0 x L) V ( x 0, x L) 1) まず区間 0からLまでに粒子をばらまく 1/ 2 2) 拡散方程式は時間 tの間に平均して x 2 (2 Dt) 1/ 2 のステップで 酔歩運動する過程を表している 各粒子を時間 tごとに P=exp(-x 2 / Dt) で与えられる確率 xで動かす 3) 箱の外に出た粒子は消滅させる 4) 粒子数 N の減衰率によって エネルギー固有値 E=-(dN/dt)/N は決定できる

25 昨年度のシミュレーション問題 2 次元 外部磁場なしのイジングモデルをモンテカルロシミュレーションするプログラムをつくってみよう 格子の大きさは20X20 周期的境界条件を用いる printf, write 文などを使って20 文字 X20 行でスピンの上下を表示する ( 例えば上向きに1 下向きに0という数字を使う + と-でもよい ) sleep() という関数を使用することでプログラムの実行を一時休止することができる kt/jの値を変えて計算してみる 特に相転移に対応するkT/Jを探ってみよ 参考 : 年度は自由課題とする 締切りは

26 レポートには使用する必要ないが curses の使用例 compile するときは -l curses が必要 #include <stdlib.h> #include <stdio.h> #include <curses.h> #include <unistd.h> int main() { int i,j,l; initscr(); for(l=0;l<10;l++) { for(i=0;i<10;i++) { move(i,0); clrtoeol(); for(j=0;j<10;j++){ } mvprintw(i,j+l,"%1d",(i+j+l)%10); refresh(); } sleep(1); } endwin(); }