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

数値計算法

2011/7/6

林田 清

モジュール化



長い長いメインプログラムは間違いのもと。また、問題を少し変

えるたびにデバッグが必要。



プログラムを複数の部品（モジュール）にわけて作成することで

、プログラムの見通しをよくする。

 デバッグの労力を減らせる。  複数の人で手分けして作業できる。  一部の部品は別のプログラムでも利用できる。 

C言語ではfunctionを利用することでモジュール化できる。

 プログラムに必要な機能を抽象化して、モジュール間のインターフェス を明確に定義するのが重要。  Functionの引数があまりにも多いようなら、インターフェースの切り分 けが適当でないと考えるべき。 

モジュール化をより発展させた“オブジェクト指向”プログラミン

グ

 C++は代表的な“オブジェクト指向”言語  抽象化、カプセル化、継承、多態  classの概念、private、publicの区別など

関数の利用

_（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; } ポインターについては本で学習してください

一般的に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); …. } 2 3 0 2 3 0

4

4

x Au y _Au x y

dv

Z Z e

x

dt

m r

dv

_{Z Z e}

_y

dt

m r

dx

v

dt

dy

v

dt

   













, , , 4 x y v v x yの 変数の連立一次微分方程式

構造体とクラス



構造体(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( 0.001 ); http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/cprogram4.html参照

電気力線の書き方



ポアソンの方程式を（数値的にor解析的に）とい

てポテンシャルを求める。



ポテンシャルの等しい面（線）をかく。



等ポテンシャル線に常に垂直になるように線を引く。



タートル法



電場を計算し、電場の方向にそって“

Turtle”をすすめ

る。

2     

















1 1 2 2 3 3 1 2 1 1 2 2 3 3 1 2 x y Q x x Q x x E r r Q y y Q y y E r r         x y 1 Q Q₂

ポアッソン分布に従う乱数

[ ]

-( )

/ !

( )

/[ ]!

( )

x y

P x

e

x

x

Q y dy

e

y

Q y

y

x

 









ただし は離散値、棄却法はそのままでは使えない

　という連続分布関数を定義する

ここで[y]はyを超えない整数

に従う乱数 を棄却法により発生させ、その整数部を

とって とする

シミユレーション

_(Simulation)



＝模擬実験



物理における計算機シミユレーション



現実の（複雑な）物理現象に対して、物理法則に基づいたモデ

ルをつくり再現を試みる。

 仮定した物理法則あるいはモデルの検証をする。  現象の予想、パラメータの最適化を行う。  （厳密解の計算が困難な場合に）近似解を求める。 

モンテカルロ法（乱数を用いた統計的手法）を利用したシミュ

レーション＝モンテカルロシミュレーション

 微分方程式を解く、あるいは連立方程式を解く問題に帰着するシ ミレーション問題も多い。 今回の参考文献：シミュレーション物理入門、矢部・川田・福田著（朝倉書店）

物理におけるシミュレーション具体例１



気体、液体の運動



流体力学の方程式を解く方法



個々の粒子の運動を追う方法（分子動力学法）



重力多体系（恒星の集団としての銀河等）の計算



荷電粒子の運動も同じ扱いが使える



このようなシミュレーションでは、通常、空間を格子（メッ

シュ）上に分割しそれぞれの場所でのポテンシャル、力等

を計算していく。



注）初期値を与えるのに乱数を使うことも多いが、乱数を用いた確

率論的方法が計算の柱になるわけではない

_{→モンテカルロシミュ}

レーションではない。

物理におけるシミュレーション具体例２



物質中での放射線の相互作用



検出器の応答のシミュレーション



イジングモデル



磁性体のモデル



量子モンテカルロ法



拡散過程の応用



これらは、モンテカルロシミュレーション。



注）モンテカルロ法はシミュレーションだけに使われるわ

けではない。 例）モンテカルロ積分。

分子動力学法１



液体のように膨大な数の分子からなる物質のマクロな性質

を決める方法（化学反応の再現等にも利用される）



基本的には、



各粒子の位置、速度の初期値を与え、



下の運動方程式を差分化した式に従って時間発展させていけば

よい。

( ) , , ij ij ji ij ij i ij j U r r dU j i F F F dr r dx i F dt dt     



 i i ２粒子間の相互作用のポテンシャル とする 粒子から 粒子に働く力は また dv 粒子の運動方程式は　 v

分子動力学法２



現実的には、計算できる

粒子の数は限られている。

そのために用いられる便

法



空間を小さなセルに分割し、

周期的境界条件を仮定す

る



近接ポテンシャルと集団ポ

テンシャルに分けて取り扱

う。

2 2 ) ( ) ( )4 2 c i ij r N U U r U r r V



 







ij c j(r <r dr

r

_c

分子動力学のシミュレーション例

 Farid F. Abraham, D. Brodbeck,

R.A. Rafey and W.E. Rudge, Phys. Rev. Lett. 73, 272 (1994).

 張力下での2次元固体のひび割 れを、100万個の原子を使ってシ ミュレーションした。シミュレーショ ンにより、観測困難な現象や、実 験困難な条件下の現象が理解で きる。 シミュレーション中には、 実験で発見された、沢山の現象 が見られた。  http://www.human.nagoya- u.ac.jp/~yasudak/sample3/MD-simulation.htmlより引用 アニメーションは左 のwebページで直接 みてください

重力多体系の計算



粒子（恒星）間に働く重力を計算し、各粒子の運動方程式を

差分的に解けば原理的には粒子の配置、運動を追うことが

できる



ポアッソン方程式を解く問題に帰着させることもできる。

4 -e i i G F m dt                    0 i 重力ポテンシャル は、Poissonの方程式 を満たす （c.f.クーロン力静電ポテンシャル は / を満たす） １）空間にメッシュを切り粒子を配置する。　　粒子の速度の初期値も与える。 ２）各メッシュ点での を計算する。 ３）差分化したPoissoの方程式（偏微分方程式）を解き、各メッシュ点での を求める ４）各粒子に働く力を により求める。　 dv ５）　 F m_i / _i, dxi dt  v を差分化した式を使い、各粒子の位置、速度を更新する。i ６）２）に戻り操作を繰り返す

銀河の衝突、月の起源



Grape Project



重力多体系の２体相互作用部分を計算する専用

計算機を開発

 http://grape.astron.s.u-tokyo.ac.jp/grape/より

 例）Two Disk Galaxies



製作P.Frisch



例）

月の起源



製作 Kokubo et al.



http://yso.mtk.nao.ac.jp/~kokubo/moon/

放射線の物質中での振る舞い

放射線粒子が物質に入射したときの振る舞いをモンテカル

ロシミュレーションする。



放射線粒子が物質中の原子と反応する位置を決定する。

 考えている反応の平均自由行程をL (L=1/nsとすると p(x)dx=exp(-x/L)/Ldxという確率分布に従う乱数xを発生させる。 

反応によって粒子が生じる場合、その放出方向、エネルギー等

を決定する。

 微分断面積の式dWq,に従って放出方向を決定するには、 dWq,という確率分布に従う乱数q,を発生させる。 

特に検出器中での振る舞いをシミュレートする場合は、検出器

の信号発生も模擬し直接実験データと比較できる出力を作成

する。

平均自由行程、断面積、到達（衝突）距離の分布



以上は一般の粒子の衝突などにも拡張できる。

( ) ( ) (0) exp( ) 1 / dx I I dI dI I dx k dI kIdx x I x I x I kx k l k       例として光が物質中で吸収される過程を考える。 光が物質に入射し微小な厚み を通過する前の強度を 、通過したあとの強度 を とする。 は と に比例することが期待されるので、ある定数 を用い て とかける。 これを積分すると、厚み を通過したあとの光の強度 は とかける。 を吸収係数、 を平均自由行程とよぶ。吸収にかかわる物質中の粒子の n s ns k 個数密度を 、それぞれの断面積を とすると の関係がある。（証明してみよ） 1 ( ) exp(- ) exp(- / ) ) x p x dx k kx dx x l dx l x x      個々の光子が（十分厚い）物質中の粒子と衝突して吸収されるまでの距離 の分布は でかける。 ＊）ただし(0 で、 が有限の範囲に限定される場合は規格化ファクタが異なる。 I I  dI dx

等方散乱のシミュレーション

      1 1 1 ( ) exp( / ) 0,1 - log( ) (2)3 0, , 0, 2 1 1 ( ) sin ( ) 2 2 l x p x dx x l l r x l r p d d p d d q    q q q q          光が物質中で粒子に衝突し等方に散乱される過程を考える。 (1)平均自由行程を とすると、衝突するまでの距離 の分布は の範囲の一様乱数 を発生させて、逆変換法で により求められる。 次元極座標での角度を極角 方位角 で指定することにすると 等方分布は 、 で表すことが             2 2 1 3 3 2 ) 0, 2 0,1 2 2 ) 0, 1 cos ( ) sin 2 1 1 sin 2 2 0,1 2 cos *) a r r b p d d p d d p d r r     q   q q q q q   q  q   q           できる。 それぞれ定義範囲で積分すると１になるように規格化されている。 方位角 の範囲の一様乱数 を発生させて、 として求められる。 極角 という変数変換をおこなうと から の範囲の一様乱数 を発生させて、 -1、 として求められる。 逆変換法に従う           0 0 1 3 3 3 1 1 ( ) cos 1 cos 2 2 cos 1 2 1 2 2 P p d r r r q q q q q q q q         _{ }       



と、累積確率 として求められる。 と -1 は同じ分布をするので結果は同じ。

素粒子実験での応用例



コンピューターシミュ

レーションによる

BELLE検出器で観測

されるB中間子の崩壊

事象。中央飛跡測定器

の電子信号（●印）か

ら発生した粒子の飛跡

を求め、運動量を割り

出します



（

http://ccwww.kek.jp/

openhouse96/B_Bell.

html

を転載）

イジングモデル１



磁性体を記述するもっとも簡単なモデル



相転移も再現できる

, - _{i j} i j H JS S   



格子状に粒子を配置し、そのスピンが上向き(1)か下向き(-1)しかないとする。 スピンは隣あう格子上のスピンとしか相互作用しない。 <i,j>は隣あう格子を意味する。　Jが正なら強磁性、Ｊが負なら反強磁性。 格子点のスピンの組み合わせで指定される状態の確率分布をモンテカルロ法により シミュレートする。

イジングモデル２

2 1 2 1 1 1 2 1 2 / exp( ) 0 2 1 m m m m m m m m r H H r f f kT s s r m m      １）各格子点にスピン(+1か-1)をランダムに与える（これを状態 とする） ２）ある 個の格子点のスピンを逆転させる（これを状態 とする） ３）状態 、状態 の実現確率の比 を計算する。 　　　 ４）　から１までの一様乱数 を発生させ、 のときこのスピンを逆転させる （状態 を採用する）。　そうでないときは状態 のまま。 ５）４）で採用された状態を新たな状態 として、２）へ戻り、次のステップを行う。m1 十分平衡状態に達するように全格子点のスピンについて何度か２）－５）の操作 を繰り返す。　そのあと、各ステップの４）で採択された状態は、状態の実現確率 を再現したものになる。

参考）イジングモデルのアプレット

http://jc.maxwell.jp/statisticalmechanics/ising2D/index.html http://www2.truman.edu/~velasco/ising.html

量子力学の固有値問題と変分法



試行関数が規格化されていない場合はどういう表式になるか



H=H

₀

+



Vとおいた摂動法との関係は？

 

2 2₂

 

   

 

2 d H x x V x x E x m dx               0 1 0 1 0 2 2 0 0 0 0 0 , ,... , ,... n n n try try n n n try n n n n n try try try E E u u Hu E u x x a u E a E a E E x E E E                







固有値を低い方から と並べ、対応する固有関数を とする。 ここで規格化された試行関数 を考える をパラメータ を含む形で仮定し、 が最小に なるようにパラメータ を決める。このときの の最小値 が に対する近似解。

量子モンテカルロ法１



時間依存のシュレディンガー方程式は虚時間の拡散方

程式に帰着できる。



拡散方程式をモンテカルロ方で解くことで、エネルギー

固有値決定する。

2 2 2 2 2 / 2 / 2 exp( / ) ( ) exp( ) ( ) V m it V m m iEt x E x E       _{    }    _{    }     2 -i t で時間を に変換すれば これは拡散定数D= の拡散方程式と同じ形 ＝ の形にかけるとき ＝ つまり は虚時間拡散方程式の減衰率として求まる。

量子モンテカルロ法２（箱型ポテンシャル

への応用例）





1/ 2 2 1/ 2 2 0(0 ) ( 0, ) (2 ) / ) x L V x x L t x D t x D t x     _{  }          １）まず区間0からLまでに粒子をばらまく ２）拡散方程式は時間 の間に平均して のステップで 酔歩運動する過程を表している。　各粒子を時間 tごとに P=exp(-で与えられる確率 で動かす ３）箱の外に出た粒子は消滅させる ４）粒子数Ｎの減衰率によって、エネルギー固有値E=-(dN/dt)/Nは決定できる

昨年度のシミュレーション問題

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

レポートには使用する必要ないが

_{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(); }