Pythonによる数値計算-半減期(Euler法)

はじめに

　入門を卒業していよいよ数値計算らしい問題を扱っていこうと思っています。今回のテーマは半減期です。不安定な放射性物質は放射線を出しながら崩壊を繰り返して安定な物質へと変化します。放射性物質の安定化に伴って放射能の強さは減少するわけですが、この放射能の強さが半分になる時間を半減期と呼びます。この半減期を数値計算で求めてみましょう。
　Pythonによる数値計算入門-Part1でも書きましたが、まずは数値計算の手法で解くことのできる数式を用意します。放射性物質の個数Nは時定数$\tau$の割合で減少すると仮定し、次の一階の常微分方程式を得ます。

\begin{equation}
\frac{dN}{dt} = -\frac{N}{\tau} \tag{1}
\end{equation}

これは式変形なしで数値計算の手法を適用できます。今回はEuler法という数値計算の手法を用いて、この一階の常微分方程式を解いてみましょう。ちなみに人間的な方法でこれを解くと次のようになります。グラフに書き出してみるとわかるでしょう。

N(t) = N_{0}e^{-t/\tau} \tag{2}

#Euler法
時刻tにおける放射性物質の$N(t)$が時間$\Delta t$後に$N(t+\Delta t)$になっているとすると、これは$\Delta t$で展開できて

N(t+\Delta t) = N(t) + \frac{dN(t)}{dt}\Delta t + O({\Delta t}^{2}) \tag{3}

式(1)を代入すると

N(t+\Delta t) = N(t) - \frac{N(t)}{\tau}\Delta t + O({\Delta t}^{2}) \tag{4}

式(4)の第三項$O({\Delta t}^{2})$は$\Delta t$が十分小さければ無視できます。このようにして時間$\Delta t$後の放射性物質の個数を時刻$t$での放射性物質の個数から計算することができるようになりました。これをEuler法と言います。
　Euler法を使って数値計算するわけですが、プログラム的には次のようになります。

for i in range(steps-1):
    t[i+1] = t[i] + dt
    n[i+1] = n[i]-n[i]/tau*dt

変数stepsは前回も説明した分割数です。分割して計算してその結果を逐次リストに格納します。時間リストtには一つ前の時刻からdtを足すことで次の時刻を定義します。粒子数リストnにはEuler法によって一つ前の時刻での粒子数から求めることができるので、その値を計算して格納します。

計算

　では全体のプログラムをみてみましょう。

decay.py

import math, pylab

#①
n_0 = 100
tau = 1
dt = 0.01
t_max = 5
steps = int(t_max/dt)
n = [0.0]*steps
t = [0.0]*steps

#②
t[0] = 0.0
n[0] = n_0

#③
for i in range(steps-1):
    t[i+1] = t[i] + dt
    n[i+1] = n[i]-n[i]/tau*dt


#④

pylab.plot(t, n)
pylab.show()

　プログラム内ではまず①で初期値を設定します。Pythonによる数値計算入門ではグラフを描くことを目標としていたので定義域を設定するという表現を用いていましたが、一般にはこれはパラメータの初期値を設定することに対応しています。②では生成した時間リストtと粒子数リストnに初期値を格納しています。③は先ほども説明したようにEuler法を使った計算です。④でグラフを書き出しています。
　半減期を出力したい場合は③のしたに次のようなプログラムを入れると標準出力されるはずです(コマンドプロンプト上に表示されること)。

for i in range(steps-1):
    if int(n[i]) == n_0/2:
        print(t[i])

粒子数リストをint型に変換しておおよそ半分になっている時間を求めて出力しています。

#最後に
　Euler法が有効なのは$\Delta t$が十分に小さいときのみであることに注意してください。そうでない場合は収束性が保証されません。Euler法は学習目的には最適ですが、本格的な数値計算にはあまり適していません。
　しかし、Euler法は数学的にもプログラム的にも理解しやすいと思いますので、このプログラムを拡張して遊んでください。標準入力から数値を代入できるようにしたり二粒子系にするなどが考えられます。あるいは別の一階線形微分方程式を持ってきてそれを解いてみるのもいいと思います。