今回は、少し前に書いた微分方程式に関する記事の続編を書きたいと思います。以前の紹介したのは、一つの微分方程式で与えられる現象の数値解析を IGOR Pro で行うものでした。取り上げた例は、
dy/dx = g(y)f(x)
∴ y = C e－kx C：積分定数
というとても簡単な系で、実は、なにも数値解析しなくても、解析解が求められるものでしたが、今回は、もうちょっとだけ難しい、複数の変数が登場する例（=連立常微分方程式）を考えてみます。
▼ [おさらい] 連立常微分方程式
2種類の物質 A と B を何か溶媒に混ぜると、中間体 C を生成し、その後に生成物 D が得られるという逐次反応（連続反応）を考えてみます。さらに、A と B の反応は可逆反応だとすると、反応式は
Fig.1 反応式
となります。ここで k1、k2、k3 は反応速度定数です。物質A～D の濃度を [A] ～ [D] と表現すると、それぞれの物質の反応速度は次のように表現できますね。
Fig.2 反応速度式
これはまさに、連立常微分方程式となっています。
IGOR Pro では（前回の繰り返しになりますが、）微分方程式の数値解を得るには IntegrateODE という操作関数を使います。
IntegrateODE derivFunc, cwaveName, ywaveSpec
derivFunc は導関数（=ユーザ定義関数）、cwaveName は定数パラメータを与えるウェーブ。ywaveSpec は計算結果を受け取るウェーブで、初期推定値はこのウェーブに与えます。
導関数はユーザ定義関数として次のように定義することができます。ポイントは物質が4種類あるので、濃度も反応速度もマトリクスとして扱っていることです。
Function ChemKinetic(pw, tt, yw, dydt)
// ウェーブと変数の定義
Wave pw // パラメータウェーブで、pw[0] = k1, pw[1] = k2, pw[2] = k3 に
// 対応します。
Variable tt // 導関数のx変数（時刻変数）です。今回は登場しませんが、
// 変数の定義は必要です。
Wave yw // yw[0]～yw[3]は物質 A,B,C,D の濃度に対応します。
Wave dydt // dydt[0]～dydt[3] はそれぞれ d[A]/dt～d[D]/dt に
// 対応します。
// 微分方程式の定義
dydt[0] = -pw[0]*yw[0]*yw[1] + pw[1]*yw[2]
dydt[1] = dydt[0] // AとBは一致します。
dydt[2] = pw[0]*yw[0]*yw[1] - pw[1]*yw[2] - pw[2]*yw[2]
dydt[3] = pw[2]*yw[2]
return 0
End
▼ 次元ラベル
濃度(yw)や反応速度(dydt)を4列のマトリクスとして扱うと、物質A～Dを同時に1つのウェーブで扱うことができて簡単ですが、4列のうちどの列がどの物質に対応しているのか、直感的にはわかりにくいです。そこで役に立つのが「次元ラベル」です。
例えば、濃度を ChemKin というマトリクスウェーブ（100行4列）で扱うことを考えます。
Make/D/O/N=(100,4) ChemKin
1列目は列IDは0、対応する物質はA、2列目の列IDは1、対応する物質はB という置き換えを頭で行うのはわかりにくいので、SetDimLabel 操作関数を使って次元ラベルを設定します。
SetDimLabel dimNumber, dimIndex, label, wavelist
wavelist は設定するウェーブです。dimNumber で指定したオブジェクト（1=列、2=レイヤー、3=チャンク）の dimIndex 番目の要素に、 label というラベルを指定します。
次の例では、1列目（列ID=0）にA, 2列目（列ID=1）に B、3列目（列ID=2）に C、4列目（列ID=3）に D というラベルを指定しています。
SetDimLabel 1,0,A,ChemKin
SetDimLabel 1,1,B,ChemKin
SetDimLabel 1,2,C,ChemKin
SetDimLabel 1,3,D,ChemKin
こうすると、例えば、物質Aの初期値を1 に指定したい場合には、
ChemKin[0][%A] = 1
という風に %次元ラベル でウェーブの列を指定することが可能になります。例えば、物質が10種類（10列）も登場する場合、物質と列IDの置き換えは容易ではありません。そんな時には、きっと便利なんじゃないかと思います。
▼ 解いてみよう
次の条件で解いてみましょう。
- 測定は10秒間隔
- A と B の初期濃度が1
- 反応速度係数がそれぞれ k1 = 0.002、k2=0.0001、k3=0.004
- 反応速度係数を与えるウェーブは kk
次のように操作します。
SetScale/P x 0,10,"s", ChemKin //測定間隔の指定
ChemKin[0][%A] = 1 //初期値の指定
ChemKin[0][%B] = 1 //初期値の指定
ChemKin[0][%C] = 0 //初期値の指定
ChemKin[0][%D] = 0 //初期値の指定
Make/D/O KK={0.002,0.0001,0.004} // 速度パラメータを与えるウェーブの定義
IntegrateODE ChemKinetic, kk, ChemKin //実行
Display ChemKin[][%A], ChemKin[][%B], ChemKin[][%C], ChemKin[][%D] // 次元ラベルを使って、グラフ表示
A と B の濃度が次第に減少し、C の濃度がいったん上昇した後で、生成物 D の濃度が次第に増加する様子が確認できますね。
Fig.3 計算結果
生成物 D の反応速度をもっと早くする（k3 = 0.05）と、ほとんど中間体の C は存在できなくなるので、グラフはもっと顕著になります。
Fig.4 生成物 D の生成速度が速い場合
▼ まとめ
IGOR Pro で連立常微分方程式の数値解析を行ってみました。連立であろうとなかろうと IntegrateODE 操作関数を使うと解析することができます。常微分方程式を解くアプローチは、速度パラメータ既知の条件で、時刻=0 から順に計算する処理でしたが、逆に、複数の式で構成される関数で、実験値から速度パラメータをフィッティングで推計する方法も、そのうち紹介できればと思います。
井上＠技術部でした。
逆に、複数の式で構成される微分を含んだままの関数(反応速度式)で、実験値から速度パラメータをフィッティングで推計する方法をぜひ紹介してください。
k2 さん、コメントありがとうございます。
> 複数の式で構成される微分を含んだままの関数(反応速度式)で、実験値から速度パラメータをフィッティングで推計する方法
→連立常微分方程式の形で表現される関数によるフィッティング方法は、本記事掲載以降、あれやこれやと調べているのですが、できそうに思われて、きわどくできない、そんな状況です。
良い方法が見つかりましたら新しい記事として掲載したいと思っておりますが、現段階では、ご期待にお応えできません。申し訳ありません。
（なお、等高線を塗りつぶす（Filled Contour）に頂いたコメントへのお返事も、このコメントで代えさせてください。）
k2さんと同じ要求は沢山あると思いますが、IGORに限らず、どこかにそのような方法論は紹介されていませんか。あるようで、ない、という感覚を私も持っています。
HH さん、コメントありがとうございます。
> k2さんと同じ要求は沢山あると思いますが、IGORに限らず、どこかにそのような方法論は紹介されていませんか。あるようで、ない、という感覚を私も持っています。
→残念ながら、一般的な方法論の有無を調べたことがありませんのでよくわかりませんが、HH さんのおっしゃるように、案外、できそうでできないテーマなのかもしれませんね。
IGOR Pro での実現性については引き続き調査しています。いい線までは行っているんですが、まだ詰め切れていません。進展があれば、ご報告したいと思います。
先日「微分方程式中のパラメーターをフィッティングで求められないか」
という話題をコメントした張本人です。
たまに訪れておりましたが、思わぬお年玉でした。ありがとうございます。
やはり上記の場合はMathematica ですよね。
理解に欠けていた部分もあり、大変参考になりました。ありがとうございます。
上記の場合をegorで処理する事ができるなら非常に有難いとかねてから思っていました。
p.s.
井上様のegorからのアプローチも楽しみにしております。
守屋様の記事に対して返信する記事の場所を間違えてしまいました。
すいません。
k2 さん、HH さん、長らくお待たせしました、微分方程式で定義される関数のフィッティングを IGOR Pro で行う手順を、ついに公開しました。ぜひご覧ください。
http://blog.hulinks.co.jp/2011/01/ig-ode-3.html