読者です 読者をやめる 読者になる 読者になる

倭算数理研究所

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

放射性崩壊の法則に従う系に対する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*}

が得られます。 

計算物理学入門

計算物理学入門

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