IGOR Pro にはいろんな解析機能とそのサンプルが付属していますが、今回は微分方程式に挑戦してみます。今の教育指導要領では高校の数学で扱われなくなったそうですが、微分方程式とは、その名の通り方程式に微分が入っている方程式のことですね(説明になっていないなぁ)。

▼ [おさらい] 変数分離による微分方程式の解法

すっかりにぶってしまっている記憶に刺激を与えるつもりで、ちょっと教科書を紐解いてみました。

dy/dx = g(y)f(x)

という方程式を考えます。この方程式は、左辺に y を x で微分した項が含まれているので「微分方程式」の仲間ですね。g(y) ≠ 0 であれば両辺を g(y) で割って、両辺に dx を掛けると

dy/g(y) = f(x)dx

となり、x を右辺に、y を左辺に分離できます。あとは、初期条件と g(y) と f(x) の複雑さ(面倒さ)次第ですが、両辺を積分するときっと解ける...気がしますよね。

例えば、次のような簡単な例を考えます。

dy/dx = -ky

変形し、積分定数を C0 として積分すると、

dy/y= -kdx
ln y = -kx + C0
∴ y = e(-kx+C0)

となります。ec0 は定数ですので改めて C とおくと、

y = C e-kx

と、解くことができます。解けると、気持ちイイです。数学も、このくらいまでなら好きでいられる気がします。

▼ でも、数値積分

こうやって、解析解がいつでも簡単に求まればよいのですが、そうは行かないのが世の常。そういう時には、初期条件からどんどん計算していって、微分方程式が表現している現象を確認する方法、数値解析(数値積分)が利用されます。化学反応の速度について次の式を考えてみましょう。

dC/dt = -kC

......上の例と同じ式ですね。解析的に解けるのですが、簡単なので数値解析についてもこの例を使います。化学反応としては、よく混合された反応槽で、物質 C が大過剰に存在する他の物質と反応して減少していく1次反応のケースに相当しますね。もちろん、C は濃度、k は反応速度定数です。

さて、IGOR Pro で微分方程式の数値解を得るには IntegrateODE という操作関数を使います(やっと IGOR Pro の話になりました)。

IntegrateODE derivFunc, cwaveName, ywaveSpec

derivFunc は導関数(ユーザ定義関数)、cwaveName は定数パラメータを与えるウェーブ。ywaveSpec は計算結果を受け取るウェーブで、初期推定値はこのウェーブに与えます。

導関数はユーザ定義関数として、都度定義してあげる必要があります。今回の例では、次のような関数の定義が必要です。

Function FirstOrder(pw, xx, yw, dydx)

Wave pw       //パラメータ(ウェーブ)
Wave yw       //yデータ(ウェーブ)
Variable xx   //xデータ(変数)
Wave dydx     //導関数(ウェーブ)

dydx[0] = -pw[0]*yw[0]

End

pw はパラーメータのウェーブで pw[0] には k が含まれます。xx は導関数を計算する x を与える変数ですが、この例では関数に x の項がありませんので使用されません。ただし、導関数としては定義に含める必要があります。 yw は y データ(濃度データ)、dydx は導関数です。

この関数は、データレンジ全体に対して実行されるのではなく、1点ずつ実行されます。今回の微分方程式は式が1つだけなので、実際にはyw[0] と dydx[0] だけが利用されるわけです。(導関数の定義方法については、また別の機会でもご紹介しますので、今回はここまでで。)

さて、定義された関数で数値解を得るには、例えば次のようなコマンドを実行します。

Make/D/O/N=101 YY      // 結果を受け取るウェーブを作成 
YY[0] = 1              // 初期値を設定 
Display YY             // グラフを作成
Make/D/O PP = {0.05}   // パラメータ k を指定
IntegrateODE FirstOrder, PP, YY   // 実行!

最初に、結果を受け取るウェーブ(YY)を作成し、次に初期値を YY[0] に与え、グラフに表示します。PP はパラメータ(k)を保持する 1 ポイントだけのウェーブを作成し、その値として 0.05 を与えます。そして IntegrateODE を実行です。

ig_ODE_01.pngFig.1 IntegrateODE の実行結果

きれいな減衰曲線が得られます。k = 0.01 すなわち PP = 0.01 として再度実行するともっと緩やかな減衰、0.1 とするともっとシャープな減衰を確認できます。

ig_ODE_02.pngFig2. パラメータ k による変化

▼ まとめ

微分方程式の数値解析を IntegrateODE という操作関数で行いました。例に上げた微分方程式が解析的に解ける系なので、あまりメリットは見当たりませんが、IntegrateODE を使うと、こんなことができるんだということは、お伝えできたかと思います。今回のポイントは、導関数をユーザ定義関数として定義するところです。
近いうちに、もうちょっと複雑な系、連立微分方程式の数値解析を IntegrateODE で行う例を紹介することにします。

 井上@技術開発部でした。

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

トラックバック

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

コメントの投稿

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

HULINKS サイトの新着情報