* 本記事はスクリプトや操作関数に関する記載が多くなっています。あらかじめご了承ください。

2010年9月の微分方程式2という記事の末尾に、

常微分方程式を解くアプローチは、速度パラメータ既知の条件で、時刻=0 から順に計算する処理でしたが、逆に、複数の式で構成される関数で、実験値から速度パラメータをフィッティングで推計する方法も、そのうち紹介できればと思います。

と(気軽に)書いたところ、たくさん反響をいただきました(改めて、感謝)。記事にしようと調査したのですが、これが(気軽な発言を悔いるに十分なほど)困難で、なかなか手順をご紹介できませんでしたが、ついに、本日御開帳ー!!

▼ おさらい

いくつかの物質が関与する反応系として、前回同様に、

ig-ode-2_1.png

を考えましょう。2種類の物質 A と B を何か溶媒に混ぜると、中間体 C を生成し、その後に生成物 D が得られるという逐次反応(連続反応)で、A と B の反応は可逆反応、というものです。もちろん、k1、k2、k3 は反応速度定数です。物質A~D の濃度を [A] ~ [D] と表現すると、それぞれの物質の反応速度は次のように表現できますね。

ig-ode-2_2.png

この(教科書的な)反応速度式(=連立微分方程式)に対して、

「反応速度パラメータが既知であれば、この連立方程式に従って計算することで、濃度データが得られますよ」

というのが前回の内容でした。しかし、実際にはそういうことは少なくて、むしろ、先日の守谷さんの記事にもありましたように、

「それぞれの初期濃度と、生成物の濃度だけが観測できるとき、それぞれの速度パラメータの推計したい」

というケースの方が現実的でしょう。

求めたいパラメータを含んだ関数が、平易な関数であれば、いわゆるフィッティングで良いわけですね。この手順については、これまでもご紹介してきました(例えば、2010/1/21 の記事)ので割愛。

しかし、この反応速度式のように関数系は微分を含む連立方程式で与えられ、一方、フィッティングできるのはその解(=積分)というような場合、非常に難儀します。

▼ なぜ、難儀してしまうのか?

上述の関数でフィッティングするための基本的なアプローチとしては、Mathematica の記事 のように、フィッティングカーブにこの連立微分方程式を含めて定義する方法を考えることになりますが、ここで、IGOR Pro の標準的なユーザ定義関数の定義方法が問題になってきます。

Function Normal_Function(pw, x) : FitFunc
    WAVE pw
    Variable x

    return pw[0] + pw[1]*exp(-x/pw[2])
End

関数 Normal_Function はとても一般的なカーブフィッティングのユーザ定義関数で、

y = a + b*Exp(-x/b)

という関数が定義されています。さて引数に注目してください。pw は Wave で宣言されているのですが、x は Variable で宣言されています。Variable で宣言できるのは値を1つだけ保持する変数。つまり、FuncFit というフィッティングを制御する操作関数から、1ポイントずつ呼び出して処理する ようにユーザ定義関数は作られています。これが大問題。

何が問題って、呼び出された関数に1ポイント分のデータしか渡されないと、関数の内部で微分とか積分とか計算できません。。

▼ All-At-Once

さんざん調べて、ついに解法を見つけました。All-At-Once というユーザ定義関数の定義方法を利用すれば OK でした!

Function AllAtOnce_Function(pw, yw, xw) : FitFunc
WAVE pw, yw, xw

yw = pw[0] + pw[1]*exp(-xw/pw[2])
End

上記のユーザ定義関数は、先ほどの関数を All-At-Once の手法で定義したものです。通常の定義方法との違いは、

  1. 引数として x データ を Wave として与えている
  2. 引数として y データ も Wave として与えている
  3. Return 文がない(Return 文を記述しても無視されてしまいます)

の3点です。この違いが非常に大きいのです。x データも y データも、ウェーブとして与えています。1ポイントずつではなく、ウェーブ全体がこの関数に与えられますので、微分のように複数のポイントの関与する計算も可能になります。

また、Return 文ではなく、yw への代入文が書かれていることにも、注目してください。この関数では、引数 xwpw はいわば関数への入力ですが、yw は出力として扱われます。y データ(yw)に直接出力されますので、Return 文は不要になるわけです。

なお、このタイプのユーザ定義関数は、ダイアログから作成することはできず、プロシージャウィンドウに直接定義する必要がありますが、回帰分析ダイアログから利用することはできます。

▼ やってみよう。

前置きが長くなりました。論より証拠。すべての手順をノーカットでお届けします。

1) 必要な関数をプロシージャウィンドウに定義します。次の二つの関数をプロシージャウィンドウに貼り付けて、プロシージャウィンドウをコンパイルしてください。

Function ODE_Function(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]
    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

関数 ODE_Function は、名称こそ異なりますが、内容は前回の記事と同じですので、今回は説明を省きますね。


Function AllAtOnce_Function(pw, yw, xw) : FitFunc
    Wave pw  // パラメータウェーブ
Wave yw // y ウェーブ(出力)
Wave xw // x ウェーブ(入力)
//濃度データと時刻データを一時ウェーブとして作成 Make/D/O/N=(numpnts(yw),4) tmpConc //初期濃度の指定 tmpConc[0][0] = 1 //A の初期濃度 tmpConc[0][1] = 1 //B の初期濃度 tmpConc[0][2] = 0 //C の初期濃度 tmpConc[0][3] = 0 //D の初期濃度 //連立微分方程式の処理 IntegrateODE/x=xw ODE_Function, pw, tmpConc //D の濃度データを出力 yw = tmpConc[p][3] End

関数 AllAtOnce_Function はフィッティングカーブを定義したユーザ定義関数。今回の目玉です。与えられた速度パラメータ(pw)と、決まった初期濃度から積分計算(IntegrateODE)を実行して濃度を計算し、y データ(yw)として返します。返されたデータは、フィッティングを担当する操作関数(FuncFit)が受け取って、当てはまり具合を検討し、十分当てはまっていれば計算は終了。そうでなければ次の速度パラメータを与えて、再度積分計算を実行するという計算が繰り返されることになります。

なお、連立微分方程式の結果を受け取る一時ウェーブである tmpConc は、入力と同じポイント数になるように numpnts 関数でポイント数を指定しています

2) 次の初期値と速度パラメータで、20 ポイントの濃度データを作ります。

初期濃度:
  • [A]0 = 1
  • [B]0 = 1
  • [C]0 = 0
  • [D]0 = 0
速度パラメータ:
  • k1 = 0.03
  • k2 = 0.002
  • k3 = 0.04

次のコマンドを順に実行します。

Make/D/O/N=(20,4) ChemKin    //濃度データ作成
SetScale/P x 0,10,"s", ChemKin    //測定間隔の指定 10 s
ChemKin[0][0] = 1    //初期値の指定
ChemKin[0][1] = 1    //初期値の指定
ChemKin[0][2] = 0    //初期値の指定
ChemKin[0][3] = 0    //初期値の指定

Make/D/O KK={0.03,0.002,0.04}    // 速度パラメータを与えるウェーブの定義
IntegrateODE ODE_Function, kk, ChemKin    //実行

// 以下3行は、グラフ表示と編集。
Display ChemKin[][0], ChemKin[][1], ChemKin[][2], ChemKin[][3]
ModifyGraph rgb(ChemKin)=(0,0,65280),rgb(ChemKin#1)=(0,0,65280),rgb(ChemKin#2)=(0,26112,0)
Legend/C/N=text0/J/F=0/A=RC/X=-0.69/Y=5.45 "\\F'Calibri'\\s(ChemKin) A, B\r\\s(ChemKin#2) C\r\\s(ChemKin#3) D"

ig-ode3-1.png


実験データっぽくなるように、濃度 D のデータに少しノイズを乗せたウェーブ Dconc を作成して、このウェーブについてフィッティングを考えることにします。

Make/D/O/N=20 Dconc = ChemKin[p][3]+gnoise(0.02)    // Y データの作成
Make/D/O/N=20 Dtime = p*10    // X データの作成
// グラフ表示と編集
Display Dconc vs Dtime
ModifyGraph mode(Dconc)=3, marker(Dconc)=19, rgb(Dconc)=(0,0,65280), msize(Dconc)=2
SetAxis/A/N=1 left

ig-ode3-2.png


3) これから求めるフィッティングパラメータを保持するウェーブ Dconf を作成し、初期値を与えます。

Make/D/O/N=3 Dconf = {0.01, 0.001, 0.01}

上記のとおり、ノイズを乗せてはいますが、「答え」が分かっていますので、当たらずも遠からずな初期値を与えることができます。実際にはこの作業が非常に重要になってくるはずです。

4) フィッティングを実行します。次のコマンドを実行!!

FuncFit/L=20/NTHR=0 AllAtOnce_Function Dconf Dconc /x=Dtime /D

このコマンドは、回帰分析ダイアログ(分析メニュー>回帰分析)を開いて、

(1) 関数とデータのタブで、関数リストから AllAtOnce_Function、Y データとして Dconc、X データとして Dtime を選択し

(2) 係数のタブで、係数ウェーブのリストで Dconf を選択し、

(3) さらに出力オプションのタブのポイントの数のボックスに 20 と入力して、

(4) 実行ボタンをクリック

しても同じ操作を実行できます。

ポイントは /L=20 です。出力ウェーブを自動設定にしています(/D)ですので、fit_Dconc というウェーブがフィッティングで得られたカーブとして生成されますが、既定の設定ではこのポイント数は200ポイントになります。しかし、All-At-Once フィッティングでは、入力ウェーブを超えるポイント数の出力ウェーブを作成しようとするとうまくいかないみたいです。なので、ここで /L=20 を指定することは、かなり重要です。

実行ボタンをクリックして、計算を開始すると、繰り返し IntegrateODE が実行されますので、ご利用の環境によっては、次のダイアログで少々待たされるかもしれません。

ig-ode3-3.png


OK ボタンがクリックできるようになれば、フィッティングは終了です。 

ig-ode3-4.png

じゃーん。こんなグラフが得られたら、大成功です。もし、上記のとおり操作してもフィッティングに失敗するようでしたら、ノイズを追加しないデータで試してみると良いかもしれません。



▼ まとめ

連立微分方程式で定義されるカーブフィッティングの方法を(ついに)ご紹介しました。これまでご紹介した記事の中では、かなり、難易度の高い内容になってしまった上に、説明はかなりは省略しています。 IGOR Pro に付属する Curve Fitting.ihf あるいは Analysis.ihf というヘルプファイルには、今回紹介した内容が詳しく記載されていますので、ぜひご参照ください。

井上@技術部でした。


注:もし、積分の精度を上げるために、IntegrateODE にもっとたくさんのデータポイントを渡す必要があるときには、別のX ウェーブ(例えば tmpX)を定義して、tmpConc と tmpX のポイント数を増やして IntegrateODE を実行させて、最後に yw へ値を返すところで間引いてあげても OK です。

※この記事の内容は執筆者の個人的見解で、ヒューリンクスによる公式情報ではありません。[免責事項]

トラックバック

この記事へのトラックバックURL
http://blog.hulinks.co.jp/cgi/mt/mt-tb.cgi/453
内容に対しての関連性がみられないものは削除する場合があります

コメント一覧

実際にやってみたのですが、
プロシージャウィンドウをコンパイルする部分で、
2つめの関数のコンパイルが出来ません。

Make/D/O/N=(numpnts(yw),4) tmpConc

の部分にエラーと出てしまいます。
どうしたらよいでしょうか?

fitman さん

本記事の内容をお試しいただきありがとうございます。関数のコンパイル中、Make 文でエラーが出るとのことですが、こちらのいろんな環境(IGOR Pro 5.05J、6.06AJ、6.12AJ、6.2E)で試してみましたが、同じ場所で失敗するケースは見当たりませんでした。

(本文にすっかり書き忘れてしまいましたが、)ご紹介しました一連の操作は、バージョン 6.1 以降での動作を確認しています。以前のバージョンではこのままの内容では動作しない可能性もございます。

残念ながら当ブログのコメント欄では、これ以上の対応はいたしかねますので、ご利用のOS 環境や IGOR Pro のバージョン、エラーメッセージ詳細などを添えて、ヒューリンクスのお問い合わせフォームをご利用ください。
https://www.hulinks.co.jp/forms/support/common.php

お手数をおかけしますが、よろしくお願いいたします。

fitmanさん

私も同じでした(igor pro 5.05j)が、コピペをすると
Wave pw // パラメータウェーブ
Wave yw // y ウェーブ(出力)
Wave xw // x ウェーブ(入力)
が改行されなかった為だと思われます。
そのために下の段の
Make/D/O/N=(numpnts(yw),4) tmpConc
にエラーが出るかと…。なんにせよ改行すれば大丈夫だと思います。

井上さん

ご無沙汰しております。
ありがとうございます。
色々試してみましたが、現実的には
「初期値を与えるよりもそのフィッティングで初期値を通るのか?」
と疑問に思われる事が多いような気がします。(私が実験系の化学屋だからでしょうか…)
ですので初期値もpwの変数に変えてしまい、
3(速度乗数)+4(初期値)の7変数でフィッティングにしてみました。
私としてはこちらの方がしっくりきました。

何にせよigorを用いて微分方程式のみでフィッティングするのは
プロの方からすれば当たり前なのかも知れませんが、
詳細に(日本語で)解説したのはwebを探す限り
ここしか見当たらないのでありがたいです。

個人としてはとても勉強になるので
このような濃い内容のこともたまには紹介していただければうれしいです。
ありがとうございました。


k2 さん:

コメントありがとうございます。なるほど、コードの改行に問題があったのですね。手元の環境では、偶然なのか、問題なく貼り付けできていたので気づきませんでした。

fitman さん:

改行位置を修正した場合、うまく動作するでしょうか? ご確認いただけると嬉しいです。


また、k2 さん:

> 色々試してみましたが、現実的には
> 「初期値を与えるよりもそのフィッティングで初期値を通るのか?」
> と疑問に思われる事が多いような気がします。(私が実験系の化学屋だからでしょうか…)

→こういった疑問はあって当然だと思います。フィッティングは、非常に便利ですが、数学的な当てはめにすぎませんので、物理的な現象にマッチしているかどうかを判断するのは、あくまでユーザ側にゆだねられているものだと思います。


> このような濃い内容のこともたまには紹介していただければうれしいです。

→励ましのお言葉、ありがとうございます。
最近はライトな話題を続けていますが、またじっくり準備して、面白い話題を提供できればと思います。何か良いネタがあれば教えてください。(私に理解できることが条件ですが。。)

コメントの投稿

Emailアドレスは表示されません。は必須項目です。
ヒューリンクス取り扱い製品の内容や購入に関するお問い合わせはヒューリンクスサイト連絡先へお願いいたします。投稿前にその他の注意事項もご覧ください。

HULINKS サイトの新着情報