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

今回はオイラー法よりも精度の高い計算法を用いたシミュレーション。　

問題

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 をもう一度行え。

問題 d にある精度が2次のアルゴリズムは、オイラー法に代わって

　　 ${ \displaystyle \begin{align*} N(t+\Delta t) = N(t) - \frac{\lambda}{2}\left(N(t) + N_e\right) \Delta t \end{align*} }$
ただし
　　 ${ \displaystyle \begin{align*} N_e =N(t) - \lambda N(t) \Delta t \end{align*} }$

という式を用いて ${ N(t+\Delta t) }$ の計算を行います。　導出は「放射性崩壊の法則に従う系に対する2次の精度のアルゴリズム」の記事参照。　設問 c をもう一度やるのは面倒なのでパス。　設問 a だけを2次のアルゴリズムで行います。

2次の精度のアルゴリズムを実装する

2次の精度のアルゴリズムを用いるには、物理系の更新方法を変更するだけです。　問題 2.7 では RadioactiveDecaySystem の evolve() メソッドを変更します。　まぁ、式で表されるとおりに計算させればいいだけです。　RadioactiveDecaySystem のみを抜き出して実装してみると以下のようになります：

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State BigDecimal t
    @State double N

    @Parameter BigDecimal dt
    @Parameter double lambda

    // 2次の精度のアルゴリズム
    @Override
    void evolve(Map<String, Object> params) {
        double deltat = dt.doubleValue()
        double Ne = N - lambda * N * deltat    // オイラー法での値
        N -= lambda / 2.0d * (N + Ne) * deltat    // 2次の精度の計算法
        t += dt
    }
}

dt を BigDecimal 型にしているので double 値を保持しておくためちょっと行数が増えてますが、やってることはたいしたことないですね。

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

上記この RadioactiveDecaySystem の実装を踏まえて、設問 a と同様のシミュレーション・スクリプトは以下のようになります：

@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.*

def N0 = 1.0d, lambda = 1.0d

// 解析解
def analytic = { double t -> N0 * Math.exp(-lambda * t) }
def N1 = analytic(1.0), N2 = analytic(2.0)
println "N(1) = $N1, N(2) = $N2"

// dt のリスト
def dtList = [0.0125, 0.025, 0.05, 0.1, 0.2]

new SimulationBuilder().parallel{

    def pngout = png(fileName:'problem2_7d1.png', overwrite:true){
        chart(title:'Problem 2.7 d', domainLabel:'dt', rangeLabel:'N/N(0)'){
            point(label:'t = 1', x:'dt', y:{ it.N - N1 })
            point(label:'t = 2', x:'dt', y:{ it.N - N2 })
        }
    }

    for(dt in dtList){
        sequence{
            simulation{
                system(new RadioactiveDecaySystem(), dt:dt, lambda:lambda)
                init(t:0.0, N:N0)
                parameterObservableSet()

                continueWhile{ it.t <= 1.0 }
                console(dataEntries:['t', 'N'], outputTiming:OBSERVE_FINAL)
                output(pngout['t = 1'], outputTiming:OBSERVE_FINAL)
            }

            simulation{
                continueWhile{ it.t <= 2.0 }
                console(dataEntries:['t', 'N'], outputTiming:OBSERVE_FINAL)
                output(pngout['t = 2'], outputTiming:OBSERVE_FINAL)
            }
        }
    }
}.simulate()

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State BigDecimal t
    @State double N

    @Parameter BigDecimal dt
    @Parameter double lambda

    @Override
    void evolve(Map<String, Object> params) {
        double deltat = dt.doubleValue()
        double Ne = N - lambda * N * deltat
        N -= lambda * (N + Ne) * deltat / 2.0d
        t += dt
    }
}

problem2_7d.groovy

実行結果は以下のようになります：

解析解との誤差が、直線ではなく（おそらく）放物線に沿って小さくなっている感じがしますね。　それ以上に、誤差の絶対値が 0.003 以下で、オイラー法の場合の 0.04 以下に比べて1桁以上小さくなっていることに注目。　2次の精度のアルゴリズムでは、オイラー法の値を計算した後に、さらに同じくらいの計算をして新しい ${ N }$ の値を得るので一見すると倍くらいの計算量がいるように思いますが、それによって精度が1桁以上よくなるのでリターンは大きいですね。