倭算数理研究所

科学・数学・学習関連の記事を、「倭マン日記」とは別に書いていくのだ!

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

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

問題

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 × 109 年なら崩壊定数はいくらか。

c. 半減期を単位として時間を測るようにプログラムを書き換えよ。 そのプログラムで 100μg の238U が最初の量の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 } になることが分かります。

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

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

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

今までのシミュレーション・スクリプトと物理量(今の場合は状態変数)の単位が違うので、変数名も t. N の代わりに tau, m を使いましょう。 238U の質量 { 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
    }
}
  • problem2_7c.groovy 直
  • コメントアウトしているチャート出力では、単位の変換を Groovy のクロージャにさせています。 これによって、シミュレーションの実行自体は適度な精度の数値で行えて、データ出力は見やすい単位にする、ということができます。

シミュレーション結果は以下のようになります:


{ \tau = 1 }{ m } の値が 50 μg (初期値の半分)になっているので、時間が半減期で測られていることが分かります。 { m = } 20 μg となるのは { \tau = 2.32 }、年に換算すると

  { \displaystyle
\begin{align*}
    t = 2.32 \times T_{1/2} \sim 1.0 \times 10^{10}
\end{align*}
}

より { 1.0 \times 10^{10} } 年となりました。

計算物理学入門

計算物理学入門

  • 作者: ハーベイゴールド,ジャントボチニク,Harvey Gould,Jan Tobochnik,鈴木増雄,石川正勝,溜渕継博,宮島佐介
  • 出版社/メーカー: ピアソンエデュケーション
  • 発売日: 2000/12
  • メディア: 単行本
  • 購入: 1人 クリック: 28回
  • この商品を含むブログ (45件) を見る
プログラミングGROOVY

プログラミングGROOVY