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

今回は問題 2.7 の続き。　放射性原子核を扱う場合、崩壊定数 ${ \lambda }$ よりも半減期 ${ T_{1/2} }$ の方がよく使われるので、こちらを元にしたコードに書き換える問題です。

問題

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} \qquad \cdots(*) \end{align*} }$
で与えられる。　修正したプログラムでこの関係を確かめよ。
　　 ${ \displaystyle \begin{align*} {}^{238}\textrm{U} \; \longrightarrow\; {}^{234}\textrm{Th} \end{align*} }$
の半減期が ${ 4.5 \times 10^9}$ 年なら崩壊定数はいくらか。
c. 半減期を単位として時間を測るようにプログラムを書き換えよ。　そのプログラムで 100μg の [\tex:{ {}^{238}U }] が最初の量の20%に減少する時間を求めよ。
d. プログラムを2次のアルゴリズムに書き換えて、設問 a と c をもう一度行え。

半減期を崩壊定数で表す (*) 式はこちらで導いています。

b. 半減期をシミュレーションで求める

まずは前回に書いた、崩壊定数を用いたシミュレーション・スクリプトを変更して半減期 ${ T_{1/2} }$ を求めてみます。　このためには、核数 N の値が、初期値（サンプル・コードでは N0 = 1.0d）の半分（サンプル・コードでは N_end = N0*0.5d = 0.5d）になるまで反復を繰り返します。

@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 N0 = 1.0d, N_end = N0*0.5d
def lambda = 1.0d

new SimulationBuilder().simulation{
    system(new RadioactiveDecaySystem(), dt:0.001, lambda:lambda)
    init(t:0.0, N:N0)

    continueWhile{ it.N > N_end }
    console(outputTiming:FINAL)
    png(fileName:'problem2_7b.png', overwrite:true){
        chart(title:'Problem 2.7 b', domainLabel:'t', rangeLabel:'N/N(0)'){
            dot(label:'Simulation', x:'t', y:'N')
            marker(label:'N / N0 = 1/2', y:0.5d, paint:Color.BLUE)
            marker(label:'T_1/2', x:log(2.0d)/lambda, paint:Color.GREEN)
        }
    }
}.simulate()

class RadioactiveDecaySystem extends AnnotationObservableSystem{

    @State BigDecimal t
    @State double N

    @Parameter BigDecimal dt
    @Parameter double lambda

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