CONTENTS

1 Runge-kutta法
2 Runge-kutta法　例題
- 2.1 ロトカ・ヴォルテラの方程式
- 2.2 ソースコード
3 まとめ

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を取得します。

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

らい・ぶらり

数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～

http://shimaphoto03.com/science/rk-method/

資料請求番号：TS33 TS38数値計算による微分方程式解法の基本～runge-kutta法～ルンゲ・クッタ（runge-kutta）法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度...

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

Runge-kutta法　例題

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

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

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

※前回Excelで解いた結果はこちら

らい・ぶらり

数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～

http://shimaphoto03.com/science/rk-method/#Lotka

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

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

こんな風に結果がホラ！

おー。。。

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

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

ソースコード

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

らい・ぶらり

初めてのExcelマクロVBA

http://shimaphoto03.com/program/macro_fund/

資料請求番号：SH45Excelマクロの始め方　まとめ本記事では、まだ一度もExcelマクロに触れたことのない人向けに、Excelマクロの始め方とExcelマクロを使用する上での最低限の知識と文法を記したいと思います。マクロを始めようなぁ、シママ、最近・・・っーつか、研...

/* —-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回だけ出力されるようになっとる！