常微分方程式の近似計算について

常微分方程式を近似計算するプログラムをしばらく書いていなかったので、復習がてらまとめてみました。

近似計算の方法はいくつかありますが、基本的にテイラー展開した式の何番目の項まで使うかという問題だと思います。

それなりの精度を必要とするシミュレーションの場合は、ルンゲクッタ法を使うことがほとんどだと思いますが、大雑把な挙動だけ分かればいいという場合にはオイラー法などを用いることもあるようです。

今回は神経科学の研究者には馴染みのある、Hodgkin-Huxley方程式の数値計算を、オイラー法、修正オイラー法、ルンゲクッタ法を用いて行ってみました。

自作のGUIに結果を表示しており、軸が表示されていませんが、横軸は時間で全体が100msec、縦軸が電位で上限が150mV、下限が100ｍV（ただし細胞外の電位ではなく細胞内の静止電位を0Vとしています）です。

印加電流を10μAにすると規則的に発火することが知られているため、その条件で数値計算を行いました。

まずは刻み幅（dt）を0.1msecにした場合の、オイラー法の計算結果。

f:id:m6chin9h96d:20150226185942j:plain

途中で電圧の値が発散してしまいました。

次に刻み幅（dt）を0.1msecにした場合の、修正オイラー法の計算結果。

f:id:m6chin9h96d:20150226190027j:plain

こちらもオイラー法とほぼ同様の結果に。

次に刻み幅（dt）を0.1msecにした場合の、ルンゲクッタ法の計算結果。

f:id:m6chin9h96d:20150226190158j:plain

こちらでは実験データの即した規則的な発火が起きています。

刻み幅が同じであれば、やはりオイラー法（修正オイラー法）よりもルンゲクッタ法の方が精度のよい結果が得られることが分かります。

最後に刻み幅を（dt）を0.01msecにした場合の、オイラー法、修正オイラー法、ルンゲクッタ法の計算結果。

f:id:m6chin9h96d:20150226190229j:plain

f:id:m6chin9h96d:20150226190251j:plain

f:id:m6chin9h96d:20150226190318j:plain

刻み幅を小さくすれば、Hodgkin-Huxley方程式の場合には、オイラー法でも実験結果を再現することができるようです。

〈参考資料〉

ほかに修正オイラー法やルンゲクッタ法の記事も。

ルンゲクッタ法による常微分方程式の数値解法

Euler法とその改良

講義参考資料　井本敬二

電脳ラボ