資料請求番号：SH43
ルンゲクッタ法とは、数値計算による常微分方程式の解法の一つであり、
特に4次精度Runge-kutta法はアルゴリズムが簡便であり、精度も高いことから、
工学者から高い支持を受けており、現在の研究活動および企業のR＆Dの現場でも良く用いられています。
本記事では、Runge-kutta法をVBAマクロに実装し、連立常微分方程式を解いてみました。
その備忘録を、実際のソースコードと共に記します。
以下に示される常微分方程式を解くことを考えます。
Runge-kutta法は、xi周りで4種類の傾きを取得し、その傾きによるyの増加量を重み付き平均することにより、yiからyi+1を取得します。
ルンゲ・クッタ法の詳細はこちら！
数値計算の基本から詳しく説明しています！
| 資料請求番号：TS33 TS38数値計算による微分方程式解法の基本～runge-kutta法～ルンゲ・クッタ（runge-kutta）法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度... 数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～ - らい・ぶらり |
今回は、ソースコードが主役なので、Runge-kutta法の理論に関する詳しい勉強は上のページでお願いします。
何？前回Excelで解いたロトカ・ヴォルテラの方程式のVBAマクロを作ったんだって？
※ロトカ・ヴォルテラの方程式（生態系における被食者と捕食者の生存闘争モデル）
x：ウサギ（被食者）の個体数 y：キツネ(捕食者)の個体数 t：時間
※前回Excelで解いた結果はこちら
| 資料請求番号：TS33 TS38数値計算による微分方程式解法の基本～runge-kutta法～ルンゲ・クッタ（runge-kutta）法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度... 数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～ - らい・ぶらり |
そうよ！これでExcelシートがかなりスッキリしたし、解きたい方程式が変わった時のプログラム変更も比較的ラクにできるようになったよ！
この状態でRunge-kuttaのボタンを押すと・・・
こんな風に結果がホラ！
おー。。。
で、中身はどうなってるんだ？
（もうちょっとこう、なんかないの・・・？なんか・・・！！？）
ソースコードはこんな感じ！割と簡単でしょ？
初心者の方はVBAマクロの基本の記事を読みながら使ってみてね！
| 資料請求番号：SH45Excelマクロの始め方 まとめ本記事では、まだ一度もExcelマクロに触れたことのない人向けに、Excelマクロの始め方とExcelマクロを使用する上での最低限の知識と文法を記したいと思います。マクロを始めようなぁ、シママ、最近・・・っーつか、研... 初めてのExcelマクロ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個に調整できるのよ。
この仕組みの詳しい解説はこっちの記事にまとめといた。
| 資料請求番号：SH43 SH44n回の計算に1回だけ結果を取り出す方法VBAでルンゲ・クッタ法のような数値計算を行う中で、刻み幅を細かくしていくと、その分だけ計算結果の数が多くなります。例えば、0から10秒まで計算したいとき、刻み幅を0.1にすれば100個、0.01にすれば... エクセルVBAによる計算結果を間引きする - らい・ぶらり |
いや、何て言うのかな・・・。イマドキのマシンだったら1万回計算しても10万回計算しても時間かからんけど、
10万行の計算結果が出るのが、ホンマせからしくてな。
ワタシもそう思う！グラフは重いし参照範囲は変えないといけないし・・・。
だからこんな風に書いてみたのよ～！！
いやぁ～。これは助かるよ！ありがとう！
いえいえ～！！
このVBAのソースコードを使って、色々な方程式が解けると思うから、ぜひコピペして色々試してみてね！
今回はRunge-kutta法による微分方程式の計算をVBAでやってみましたという内容の備忘録でした。