エクセルVBAを用いた微分方程式の計算(ルンゲクッタ法)

資料請求番号:SH43

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

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

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

Runge-kutta法

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

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

シママ

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

シママ

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

Runge-kutta法 例題

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

ストーク

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

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

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

シママ

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

シママ

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

シママ

こんな風に結果がホラ!

ストーク

おー。。。

ストーク

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

シママ

(もうちょっとこう、なんかないの・・・?なんか・・・!!?)

ソースコード

シママ

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

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

シママ

気づいた?そうなのよ~。
1000日を刻み幅0.1日で計算するから計算回数は1万回。
でも、h2というパラメータを10にすれば、計算の刻み幅は0.1日のまま、出力を10回の計算に1回、つまり1日ごとに出来てプロットは1000個に調整できるのよ。
この仕組みの詳しい解説はこっちの記事にまとめといた。

ストーク

いや、何て言うのかな・・・。イマドキのマシンだったら1万回計算しても10万回計算しても時間かからんけど、
10万行の計算結果が出るのが、ホンマせからしくてな。

シママ

ワタシもそう思う!グラフは重いし参照範囲は変えないといけないし・・・。
だからこんな風に書いてみたのよ~!!

ストーク

いやぁ~。これは助かるよ!ありがとう!

シママ

いえいえ~!!

シママ

このVBAのソースコードを使って、色々な方程式が解けると思うから、ぜひコピペして色々試してみてね!

まとめ

今回はRunge-kutta法による微分方程式の計算をVBAでやってみましたという内容の備忘録でした。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です