«前の日記(2012年05月03日) 最新 次の日記(2012年05月16日)» 編集

【小波の京女日記】


2012年05月09日 放射性物質の基礎的な振る舞いを記述する (2) [長年日記]

_ Euler法で減衰のシミュレーション

放射性物質というのは,不安定な原子核をもつ原子たちの集団のことだ。 もっとも,たとえば放射性セシウムとかいう場合には,純粋のセシウムだけが あるのではなく,何らかの化合物となっているので,他の原子もいっしょに 存在しているのだけど,単に放射線のことを考えるときには,それは考えなくてもよい。

それではいよいよ,放射線の量が減っているようすをシミュレーションで確認していくことにしよう。

繰り返し計算で攻める

いま時刻$t$を 0 として,ある放射性原子核の集まりがあったとして,その量が $x_0$であったとす る。その中で $\lambda$ の割合が単位時間あたりで崩壊していくとすると, 短い時間 $\Delta t$ の間に減る量は $\lambda x_0 \Delta t$ となるだろう。これを $x_0$ から引けば残る量が求められる。 式の上でこの作業を繰り返してやると,次のようにどんどん式はややこしく 混み合ってしまう。


時刻 $t$ 量 $x$ 減少量
0 $x_0$ $\lambda x_0 \Delta t$
$\Delta t$ $x_0 - \lambda x_0 \Delta t$ $\lambda (x_0 - \lambda x_0 \Delta t) \Delta t$
$2\Delta t$ $x_0 - \lambda x_0 \Delta t -\lambda (x_0 - \lambda x_0 \Delta t) \Delta t $ めんどくさ!
$3 \Delta t$ ... ...

が,よく見ると,これは次のような単純な作業の繰り返しにすぎない。

  1. 時間の進み幅 $\Delta t$, 量の初期値 $x_0$ を適当な値に設定する
  2. 時刻$t$をゼロに,量$x$を初期値 $x_0$をセットする
  3. 量 $x$に $\lambda\Delta t$を掛けて減少量を得る
  4. 時刻を $\Delta t$だけ進める
  5. $x$から減少量を引いて,現時刻での新しい $x$とする
  6. 3 に戻る

勘違いを避けるために注意しておくと,数学の時間にやらされるような数式の操作を実行することは,コンピュータにはむつかしい(因数分解とか式の展開とかできるようなソフトはある)。そういう高級なテクニックは使わずに,具体的な値を使って「どん臭く」計算を進めるのが,この種の数値シミュレーションである。

Ruby のソースで表現する

上の繰り返しをRubyのソースで書くとどうなるだろうか? 変数にギリシャ文字は使えないので,$\lambda$ は lambda, $\Delta t$は dt というふうに置き換えて書いてみる。もちろん添字も使えないので $x_0$ は x0 のようにする。

lambda = 0.1
dt = 0.05
x0 = 10.0
x = x0
t = 0.0
#繰り返しはじめ
dx = lambda * x * dt
t += dt
x -= dx
#繰り返しおわり

繰り返しには for ループや while ループがあるが,ここではどちらを使えばよいだろうか。比較すると,

  • for ループは回数を指定できる
  • while ループは終了条件を指定できる

ということはよく知っているだろう。回数指定のほうが頭を使わないですむので,まずそれでやってみよう。上の繰り返しの部分を次のように for ループで囲む。n_repeat は繰り返し回数であり,最初のほうで適当に 100 とか 1000 とかを代入しておく。

for i in 1 .. n_repeat #繰り返しはじめ
  dx = lambda * x * dt
  t += dt
  x -= dx
end #繰り返しおわり

一応これでプログラムは走るはずだが,計算結果を外に吐き出す部分がないので,ただ黙ってひたすら走って終わる。つまり,

本人は分かっているのだろうけど,それじゃ他人はわからないよね。

って感じだよね。

せっかく計算して得た数値を出力するには,puts なり printf なりを使えばよいが,この種の計算では printf を使って,出力を精密に制御するほうがよい。

出力されたデータの形は次のように GNUPLOT でそのまま使えるスペース区切りのフォーマットにする。

0.000000  10.0000000
0.010000   9.9800000

課題: それではソースを補って,上のような出力を得られるようにして下さい。プログラム名は radiation_decay01.rb とし,適当にバージョンを上げながら改良するとよいでしょう。

得られたグラフをみる

結果を GNUPLOT で描画してみると次のようになる。

画像の説明


課題

  • パラメータ lamda をいろいろ変更して上の計算を行い,その結果を同一のグラフの上にプロットしなさい。グラフ右上の線の説明にはファイル名でなく,パラメータを表示するようにすること。
  • この種の減衰では,いわゆる半減期が現れる。量が半減する時間が,どんな量から始めても同じであることをこのシミュレーションで確認しなさい。
  • 半減期と $\lambda$ の関係はどうなっているか。log(2) の値を参考にしてシミュレーションの結果から推測しなさい。
  • 結果のグラフと共に簡単な説明を付けて,何らかの文書にしてメールに添付して提出しなさい。文書作成は Word, OpenOffice, TeX などなんでもよいが,もらうのはPDF の形がうれしい(OpenOffice には PDF でエクスポートする機能がある)。

" log(2) は2の自然対数である。irb で include Math としてからこの式を打ち込むと得られる。また Windows などの関数電卓でも容易に計算できる。