化学反応シミュレーション


実験で振動反応を間近で観察することができたと思います。これを計算機の中でシミュレーションしてみましょう。

1次反応

いきなり振動反応に飛び込んでいくのはしんどいので、まずはもっとも単純な反応系から考えていきます。次のような化学反応を考えましょう。

A --> P

つまり、"A"(反応物)という物質が反応して"P"(生成物)という物質に変化していく反応です。1つの化合物が変化するときの化学反応を"1次反応"と呼びます。今、"A"という物質がどういう割合で"P"に変わっていくか、が興味があります。つまり"A"の時間に対する濃度変化が知りたいわけです。一般に1次反応の場合、反応物の濃度変化は反応物の濃度に比例します。

d[A]/dt = - k[A]

ここで、[A]は"A"の濃度、d[A]/dtは"[A]"の濃度変化、kは比例定数を表します。右辺にマイナスがついているのは、時間とともに"A"の濃度は減少するからです。定数kを通常速度定数と呼びます。数学の言葉でいうと、上の式は1階の微分方程式です。

上の式をもう少し変形してみましょう。ここで、左辺は濃度変化ですので、ある時間tの時の濃度を[A]t, 適当な時間たったときの濃度を[A]t'としましょう。dtとは時間間隔をあらわしますので、上の式は

([A]t' - [A]t)/(t'-t) = -k [A]t

[A]t' = [A]t - k [A]t (t'-t)

と変形できます。つまり、時刻tの時の濃度がわかれば、適当な時間間隔後(あまり大きすぎるとよろしくない)の濃度がわかるわけです。このような操作のことを"微分方程式の差分化"と呼びます。計算機で微分方程式を解くための1つのテクニックです。

実は、上の微分方程式は比較的簡単に解くことができます。差分化して解いた答えとあわせるために、その解を示しておきます。

[A] = [A]0 exp(-kt)

ここで[A]0は"A"の初期濃度を表します。

それでは計算機で解いてみましょう。以下にFORTRANのプログラムを示します。

program first
implicit real*8(a-h,o-z)

a0 = 1.0
dk = 1.0
dt = 0.01

t = 0.0
a_sabun = a0

do i = 1 , 1000
t = t + dt
a_exact = a0 * exp( - dk * t )
a_sabun = a_sabun - dk * a_sabun * dt
write(6,'(f10.5,3f20.10)') t, a_exact, a_sabun, a_exact - a_sabun
end do

end

これを仮に"first.f"という名前で保存し、コンパイルしましょう。

f77 first.f
a.out > first.txt

計算出力を"first.txt"に保存するようにしましょう。そうしないと、たくさんの数字が画面に流れていってしまいます。

計算出力には正確な解と差分で求めた解との差も出力するようにしています。まず今の条件でどの程度正確に計算できたか、確認しましょう。

less first.txt

次に、出力ファイルをwindowsにftpで持ってきて、gnuplotを使ってグラフ化してみましょう。


逐次反応

”ちくじ”と読みます。次のようなもう少し複雑な反応を考えましょう。

A --> B --> C

このとき、"A", "B", "C"の各濃度変化は次のように書き表すことができます。

d[A]/dt = - k1[A]

d[B]/dt = k1[A] - k2[B]

d[C]/dt = k2[B]

1次反応の時と同様にして上の方程式を差分化して解いてみましょう。また、この反応系についても正確な解を求めることができます。解法については大学レベルの数学が必要なのでここでは省略します。これを仮に"cons.f"という名前で保存し、コンパイルしましょう。

f77 cons.f
a.out > cons.txt

今度は3つの成分があるので、正確な解と差分による解の差は出力していません。gnuplotでグラフ化するところでどの程度正確に計算できたかを確認してください。

1次反応のときと異なり、速度定数の比を変えることによってグラフの様子を変えてやることができます。いろいろ変えてみましょう。


いよいよ振動反応

いよいよ振動反応にチャレンジしてみましょう。いくつか反応系がありますが、ここではロトカ−ヴォルテラ機構(Lotka-Volterra mechanism)についてシミュレーションしてみましょう。

A --> 2A

A + B --> 2B

B --> C

化学反応系として考えるのが難しいのであれば、"A"をうさぎ、"B"をおおかみ、"C"をおおかみの死と考えてもらってよいでしょう。つまり"A"は食べられる側、"B"は食べる側で、食物連鎖の関係にあります。"A", "B"の速度式は次のようになります。

d[A]/dt = k1[A]2 - k2[A][B]

d[B]/dt = k2[A][B] - k3[B]

今までと同様、差分化して問題を解きましょう。ただし、今回は正確な解はありません。"vib.f"という名前で保存し、コンパイルしましょう。

f77 vib.f
a.out > vib.txt

先ほどと同様、速度定数の比をいろいろ変えながら、グラフの様子がどう変わっていくか、調べてみましょう。