【Excel VBA入門】ルンゲクッタ法で微分方程式を解く！サンプルコード付き解説

Excel VBAマクロでRunge-kutta法を使い、常微分方程式を解く

ルンゲクッタ法とは、数値計算による常微分方程式の解法の一つであり、
特に4次精度Runge-kutta法はアルゴリズムが簡便であり、精度も高いことから、
工学者から高い支持を受けており、現在の研究活動および企業のR＆Dの現場でも良く用いられています。

本記事では、Runge-kutta法をVBAマクロに実装し、連立常微分方程式を解いてみました。
その備忘録を、実際のソースコードと共に記します。

Runge-kutta法

以下に示される常微分方程式を解くことを考えます。

Runge-kutta法は、x_i周りで4種類の傾きを取得し、その傾きによるyの増加量を重み付き平均することにより、y_iからy_i+1を取得します。

ルンゲ・クッタ法の詳細はこちら！
数値計算の基本から詳しく説明しています！

[blogcard url=”http://shimaphoto03.com/science/rk-method/”]

今回は、ソースコードが主役なので、Runge-kutta法の理論に関する詳しい勉強は上のページでお願いします。

Runge-kutta法　例題

ロトカ・ヴォルテラの方程式

何？前回Excelで解いたロトカ・ヴォルテラの方程式のVBAマクロを作ったんだって？

※ロトカ・ヴォルテラの方程式（生態系における被食者と捕食者の生存闘争モデル）
x：ウサギ（被食者）の個体数　y：キツネ(捕食者)の個体数　t：時間

※前回Excelで解いた結果はこちら
[blogcard url=”http://shimaphoto03.com/science/rk-method/#Lotka”]

そうよ！これでExcelシートがかなりスッキリしたし、解きたい方程式が変わった時のプログラム変更も比較的ラクにできるようになったよ！

この状態でRunge-kuttaのボタンを押すと・・・

こんな風に結果がホラ！

おー。。。

で、中身はどうなってるんだ？

（もうちょっとこう、なんかないの・・・？なんか・・・！！？）

ソースコード

ソースコードはこんな感じ！割と簡単でしょ？
初心者の方はVBAマクロの基本の記事を読みながら使ってみてね！

[blogcard url=”http://shimaphoto03.com/program/macro_fund/”]

/* —-Lotka-Volterra-rk———- */

Public a, b, c As Double
Sub rklotka()
Dim x0, y0, init, ed, h, h2, kx(4), ky(4) As Double
Dim i, j As Integer

‘Calculate property
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
a = Cells(5, 3)
b = Cells(6, 3)
c = Cells(7, 3)

‘Initial Condition
x0 = Cells(4, 7)
y0 = Cells(4, 8)
x = x0
y = y0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = x
Cells(j, 8) = y

End If

kx(1) = h * F1(t, x, y)
ky(1) = h * F2(t, x, y)
kx(2) = h * F1(t + h / 2, x + kx(1) / 2, y + ky(1) / 2)
ky(2) = h * F2(t + h / 2, x + kx(1) / 2, y + ky(1) / 2)
kx(3) = h * F1(t + h / 2, x + kx(2) / 2, y + ky(2) / 2)
ky(3) = h * F2(t + h / 2, x + kx(2) / 2, y + ky(2) / 2)
kx(4) = h * F1(t + h, x + kx(3), y + ky(3))
ky(4) = h * F2(t + h, x + kx(3), y + ky(3))
nx = x + (kx(1) + 2 * kx(2) + 2 * kx(3) + kx(4)) / 6
ny = y + (ky(1) + 2 * ky(2) + 2 * ky(3) + ky(4)) / 6
nt = t + h

t = nt
x = nx
y = ny

Next i

End Sub

Function F1(ByVal t As Double, ByVal x As Double, ByVal y As Double) As Double
F1 = a * x – c * x * y
End Function

Function F2(ByVal t As Double, ByVal x As Double, ByVal y As Double) As Double
F2 = -b * y + c * x * y
End Function

/* —————————– */

とりあえず、方程式が変わったら、下にあるF1やF2だけを変化させればいいからわかりやすくなったな。
自分で考えた方程式がサッと試せてええな。

そうなのよ～。VBAを使うとこうやって自分で関数を定義することができる（※）から、単にExcelで計算する場合に比べれば可読性が増すと思うわ。

※このことを「Functionプロシージャを使う」と言います。

あ、これ、10回に1回だけ出力されるようになっとる！