放射性崩壊の法則に従う系に対する2次の精度のアルゴリズム

放射性崩壊の法則

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

• \( t \) ： 時間
• \( N \) : 核の個数
• \( \lambda \) ： 崩壊定数

に従う系に対して、数値計算法を2次の精度で行う方法を導きます。　数値計算法の精度が \( n \) 次とは

　　\( t \) の値を一定にしたときに解析解と数値解の差が \( (\Delta t)_n \) に比例するなら、その数値計算法は \( n \) 次であるという。

で定められます。　ここで導く方法は崩壊の法則の微分方程式 (1) を使うため、以前に導いた Runge-Kutta （ルンゲ=クッタ）法や高次テイラー法に比べて適用範囲がかなり限定されます。

導出

時刻 \( t + \Delta t \) での核数 \( N(t+\Delta t) \) を \( \Delta t \) の冪に展開すると
\begin{align*}
N(t+\Delta t) = N(t) + \frac{dN(t)}{dt}\Delta t + \frac{1}{2}\frac{d^2N(t)}{dt^2}(\Delta t)^2 + \cdots
\end{align*}
となります。　この式を \( n \) 次まで残して計算すると \( n \) 次の精度の数値計算法が得られます。　これは高次テイラー法そのものですね。　ここまでは \( N(t) \) がどんな微分方程式に従う場合でも使えます。

次に \( \Delta t \) の2次までを残して、かつ崩壊の法則 (1) 式を使って \( N(t) \) の微分を書き換えると
\begin{align*}
N(t+\Delta t) &= N(t) - \lambda N(t) \Delta t + \frac{\lambda^2}{2} N(t) (\Delta t)^2 & \cdots (2)
\end{align*}
となります。　崩壊の法則 (1) を使っているので、これ以降は微分方程式が (1) の形の系にしか適用できません。

一方、オイラー法で計算した \( \Delta t \) 秒後の \( N \) の値を \( N_e \) とすると
\begin{align*}
N_e &= N(t) - \lambda N(t) \Delta t & \cdots (3)
\end{align*}
です。　これは (2) 式で \( \Delta t \) の1次まで残したものですね。

(3) 式を使って (2) 式の \( \Delta t \) の2次の項を消去しましょう。　(3) 式の両辺に \( \lambda \Delta t/2 \) をかけて
\begin{align*}
\frac{\lambda}{2} N_e = \frac{\lambda}{2}N(t) - \frac{\lambda^2}{2}N(t)(\Delta t)^2
\end{align*}
この式と (2) 式の辺々を加えて、\( N(t+\Delta t) \) について解くと

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

が得られます。　

計算物理学入門