2.10 原子核の崩壊：問題 2.7 単一核種の崩壊 c

今回も放射性崩壊のシミュレーションの続き。　前回の最後に、²³⁸U の半減期や崩壊定数が桁の極端な値になるねって話がありました。　今回はそれを上手く扱うやり方を見ていきます。

問題

a. プログラム cool を修正し、原子核崩壊で使われる表記に直せ。　そして、 ${ \lambda = 1 }$ について、 ${ t = 1 }$ と ${ t = 2 }$ のときの ${ N(t) / N(0) }$ の解析解の値とオイラー法の差を計算せよ。　さらに、それぞれの ${ t }$ の値では、その差が ${ \Delta t }$ に比例することを示せ。
b. 放射性崩壊の割合は、 ${ \lambda }$ ではなく半減期 ${ T_{1/2} }$ 程度の大きさを単位として与えると都合がよい。　半減期は最初の原子核の半数が崩壊する時間であり、
　　 ${ \displaystyle \begin{align*} T_{1/2} = \frac{\log 2}{\lambda} \end{align*} }$
で与えられる。　修正したプログラムでこの関係を確かめよ。
　　 ${ \displaystyle \begin{align*} {}^{238}\textrm{U} \; \longrightarrow\; {}^{234}\textrm{Th} \end{align*} }$
の半減期が 4.5 × 10⁹ 年なら崩壊定数はいくらか。
c. 半減期を単位として時間を測るようにプログラムを書き換えよ。　そのプログラムで 100μg の²³⁸U が最初の量の20%に減少する時間を求めよ。
d. プログラムを2次のアルゴリズムに書き換えて、設問 a と c をもう一度行え。

今回は問題 c をやります。

半減期を単位とした時間への書き換え

半減期 ${ T_{1/2} }$ を単位として測った時間を ${ \tau }$ （タウ）としましょう。　式で書くと

　　 ${ \displaystyle \begin{align*} \tau = \frac{t}{T_{1/2}} = \frac{\lambda}{\log 2}t \end{align*} }$

です。　この ${ \tau }$ を使って放射性崩壊の法則

　　 ${ \displaystyle \begin{align*} \frac{dN}{dt} &= -\lambda N & \cdots (*) \end{align*} }$

を書き換えます。

${ \tau }$ の定義より

　　 ${ \displaystyle \begin{align*} \frac{d\tau}{dt} = \frac{\lambda}{\log 2} \end{align*} }$

なので、微分方程式 (*) の左辺は

　　 ${ \displaystyle \begin{align*} \frac{dN}{dt} = \frac{dN}{d\tau}\frac{d\tau}{dt} = \frac{\lambda}{\log 2}\frac{dN}{d\tau} \end{align*} }$

となります。　これを用いて (*) 式は

　　 ${ \displaystyle \begin{align*} \frac{\lambda}{\log 2}\frac{dN}{d\tau} = -\lambda N \end{align*} }$

よって

　　 ${ \displaystyle \begin{align*} \frac{dN}{d\tau} = - N \log 2 \end{align*} }$

となります。　つまり、半減期を単位として測った時間では崩壊定数は一定の ${ \log 2 }$ になることが分かります。

問題は初期の ²³⁸U が 100μg あると指定しているので、微分方程式に現れている核数 ${ N }$ を質量 ${ m }$ に書き換えておきましょう。　 ${ N }$ と ${ m }$ は比例するので、書き換えは ${ t }$ から ${ \tau }$ に書き換えるよりも簡単です：

　　 ${ \displaystyle \begin{align*} \frac{dm}{d\tau} = -m\log 2 \end{align*} }$

シミュレーション・スクリプト

今までのシミュレーション・スクリプトと物理量（今の場合は状態変数）の単位が違うので、変数名も t. N の代わりに tau, m を使いましょう。　²³⁸U の質量 ${ m }$ は 100μg を単位に測ることにしましょう（ ${ m = 3 }$ なら 300μg）：

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.2-beta')

import org.waman.groops.builder.SimulationBuilder
import org.waman.groops.simulation.system.*
import java.awt.Color
import static java.lang.Math.*

def m0 = 1.0d, m_end = m0 * 0.2d    // 反復終了の質量 m_end は初期質量 m0 の 20%

new SimulationBuilder().simulation{
    system(new RadioactiveDecaySystem(), dtau:0.01)
    init(tau:0.0, m:m0)

    continueWhile{ it.m > m_end }
    console(outputTiming:FINAL)
    png(fileName:'problem2_7c.png', overwrite:true){
        chart(title:'Problem 2.7 c', domainLabel:'\u03C4 [4.5*10^9 year]', rangeLabel:'mass [\u03BCg]'){
            dot(label:'Simulation', x:'tau', y:{ it.m * 100 })
        }
    }
    //時間を年単位で出力
    //png(fileName:'problem2_7c2.png', overwrite:true){
    //    chart(title:'Problem 2.7 c', domainLabel:'time [year]', rangeLabel:'mass [\u03BCg]'){
    //        dot(label:'Simulation', x:{ it.tau * 4.5E9d }, y:{ it.m * 100 })
    //    }
    //}
}.simulate()

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State BigDecimal tau    // 前の t
    @State double m            // 前の N

    @Parameter BigDecimal dtau

    @Override
    void evolve(Map<String, Object> params) {
        m -= log(2.0d) * m * dtau.doubleValue()
        tau += dtau
    }
}