原子核の崩壊のシミュレーション 〜 『計算物理学入門』読み(その8)
ここで、λは崩壊定数である。
この方程式の形は、ニュートンの冷却の法則のシミュレーションとほぼ同じ形をしている。
それ故、放射性崩壊のシミュレーションのプログラムは以下の通りの短いものになる。
(define (decay Lamda time N dt)
(euler 0.0 N time (lambda (t N) (- (* Lambda N))) dt))
λを直接指定するよりも、半減期を指定できた方が便利である。半減期は以下で与えられる。
(define (decay-half-reduce T_1/2 time N dt)
(decay (/ 2 (exp T_1/2)) time N dt))
オイラー法は1次の精度しかないので、2次の精度で計算してみる。2次の方程式は以下の通り。
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で、解析解、オイラー法と比較してみると
あれ? なんか値が解析解に収束していかないぞ。誤差が溜まってきてしまうのかなあ。まあ、途中までの結果をみると、2次の精度で求められているようだ。
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)
