So-net無料ブログ作成
検索選択

原子核の崩壊のシミュレーション 〜 『計算物理学入門』読み(その8) [計算物理学入門]

放射性原子核の崩壊の基本法則は以下のように書ける。

AtomDecayEquation.png

ここで、λは崩壊定数である。
この方程式の形は、ニュートンの冷却の法則のシミュレーションとほぼ同じ形をしている。 それ故、放射性崩壊のシミュレーションのプログラムは以下の通りの短いものになる。
(define (decay Lamda time N dt)
  (euler 0.0 N time (lambda (t N) (- (* Lambda N))) dt))

λを直接指定するよりも、半減期を指定できた方が便利である。半減期は以下で与えられる。

HalfLifeEquation.png
半減期を指定できるようにしたプログラムは以下の通り。単にラッパーの関数を作っただけ。
(define (decay-half-reduce T_1/2 time N dt)
  (decay (/ 2 (exp T_1/2)) time N dt))

オイラー法は1次の精度しかないので、2次の精度で計算してみる。2次の方程式は以下の通り。

SecondPrecisionEquation.png

2次の精度で求めるプログラムは以下の通り。

(define (decay-2nd-precision Lambda time N dt)
  (let loop ((x-i 0)
             (y-i N))
    (if (< x-i time)
        (loop (+ x-i dt) (- y-i (* (/ Lambda 2)
                                   (- (* 2 y-i) (* Lambda y-i dt))
                                   dt)))
        y-i)))

様々なdtで、解析解、オイラー法と比較してみると

gosh> (* 100 (exp (- (* 1 1))))
36.787944117144235
gosh> (map (lambda (dt)
	     (decay 1 1 100 dt))
	   '(0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001))
(31.381059608999998
 36.603234127322985
 36.76954247709651
 36.78242603282879
 36.78739229905556
 36.78792572316515
 36.787938598960544)
gosh> (map (lambda (dt)
	     (decay-2nd-precision 1 1 100 dt))
	   '(0.1 0.01 0.001 0.0001 0.00001 0.000001 0.0000001))
(33.352959127436435
 36.78856187161925
 36.78795025306916
 36.78426556798378
 36.7875762401565
 36.78794411715042
 36.78794043835247)
あれ? なんか値が解析解に収束していかないぞ。誤差が溜まってきてしまうのかなあ。まあ、途中までの結果をみると、2次の精度で求められているようだ。


nice!(0)  コメント(0)  トラックバック(0) 
共通テーマ:学問

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。

×

この広告は1年以上新しい記事の更新がないブログに表示されております。