資料請求番号：SH43 SH44
VBAでルンゲ・クッタ法のような数値計算を行う中で、刻み幅を細かくしていくと、その分だけ計算結果の数が多くなります。
例えば、0から10秒まで計算したいとき、刻み幅を0.1にすれば100個、0.01にすれば1000個、0.001にすれば1万個の計算結果を得ます。
刻み幅をな～細かくしていくと、それだけ精度の良い計算ができるし、イマドキのパソコンだと10万回や100万回くらい簡単にループできる。
でもな～。10万個のプロットでグラフを書くのは、なんつーか、だるい。100万個とかありえへん。
あ～わかるな～その気持ち。100万個のプロットとか描いたことないけど。
グラフ重いし、刻み幅を変更したらいちいちグラフの仕様を変更しなくちゃいけないのがメンドクサイ。
こういう作業をいちいちやるのも嫌だけん、何とかならん？
なるわよ！要は、n回に1回だけ出力するプログラムを書けばいいんでしょ？
ワタシが書いたやつ見せてあげる！
というわけで、本記事では、VBAにおいて、n回に1回だけ出力するプログラムの備忘録を残したいと思います。
今回は、振動の運動方程式を例題にしてみるね。
まず、この状態でワークシートを作って・・・
プログラムを実行すると・・・
こんな風に結果が出るの！スゴイでしょ？
これは・・・。50秒間積分するために刻み幅を0.01秒にしているから5000回ループしている。
でも、10回に1回しか出力しないように設定しているから、グラフのプロット数は500個だけなんやな。
そういうこと！
ソースコードは以下の通り。このコードは前回の振動計算の記事と同じものよ。
/* –vibration formula rk-method ——– */
Public omega0, gamma, force, omega As Double
Sub rkvib()
Dim x0, v0, init, ed, h, h2, kx(4), kv(4) As Double
Dim i, j As Integer
‘Calculate propertires
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3) ‘h2回に1回出力
‘Constants
omega0 = Cells(5, 3)
gamma = Cells(6, 3)
force = Cells(7, 3)
omega = Cells(8, 3)
‘Initial Condition
x0 = Cells(4, 7)
v0 = Cells(4, 8)
x = x0
v = v0
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) = v
End If
kx(1) = h * F1(t, x, v)
kv(1) = h * F2(t, x, v)
kx(2) = h * F1(t + h / 2, x + kx(1) / 2, v + kv(1) / 2)
kv(2) = h * F2(t + h / 2, x + kx(1) / 2, v + kv(1) / 2)
kx(3) = h * F1(t + h / 2, x + kx(2) / 2, v + kv(2) / 2)
kv(3) = h * F2(t + h / 2, x + kx(2) / 2, v + kv(2) / 2)
kx(4) = h * F1(t + h, x + kx(3), v + kv(3))
kv(4) = h * F2(t + h, x + kx(3), v + kv(3))
nx = x + (kx(1) + 2 * kx(2) + 2 * kx(3) + kx(4)) / 6
nv = v + (kv(1) + 2 * kv(2) + 2 * kv(3) + kv(4)) / 6
nt = t + h
t = nt
x = nx
v = nv
Next i
End Sub
Function F1(ByVal t As Double, ByVal x As Double, ByVal v As Double) As Double
F1 = v
End Function
Function F2(ByVal t As Double, ByVal x As Double, ByVal v As Double) As Double
F2 = -omega0 ^ 2 * x – 2 * gamma * v + force * Cos(omega * t)
End Function
/* –vibration formula rk-method End——– */
太字の部分見てみて！
刻み幅hに対して、新たに「h2」という変数を用意して、
h2回に1回結果を出力するようにしたの。
太字の部分の・・・・j = 6 + i / h2ってなんだ？
6は位置合わせ。計算結果が6行目から出てるでしょ？
それで、i/h2っていうがポイント！
えっと、h2=10なら、iが0のとき、jは6,
iが1のとき、jは6.1
iが2のとき、jは6.2
・・・
iが10のとき、jは7。
あ、これでi=10の時の結果を7行目に、i=20の時の結果を8行目に・・・ってことか。
そうそう！そして、jって整数型（integer）だから結局i=1~9のとき6か7になるんだけど、
i=1~9のときに出力しないでね！って言うのが次のif文。
If i Mod h2 = 0 Then・・・か。
確かにi=1~9のとき、このifはFALSEになるたいね。
i=0,10,20・・・のとき、iがh2で割り切れてTrueになって、
i=0のときj=6になって6行目、i=10のときj=7になって7行目、i=20のときj=8になって8行目に出力されるのか～。
あ～このちょっと、こそばゆいところに手が届く感じええなぁ～。
いいでしょ～。
ラクできるとこはラクしなきゃね～。マクロってそういうものだからね～。
これなら、新しい方程式を100万回くらい計算させて、10万回計算させた結果と照合して問題なければ、「その方程式に関しては刻み幅を大きくしてもよい」という判断ができるから、ムダがないな。非線形性の大きい問題にもチャレンジできる。
100万回計算させている間にタバコ吸いに行ける（※）からさらにムダがない。
※Corei5-2400, 8GBで100万回計算させるのに必要な時間は3分4秒でした。
ひせんけい？
2次関数とかexpとか非線形やろ？直線でないってことたい。
それで、こういうのって、変化が激しい時は激しいから、線形に近似する数値計算、オイラー法もルンゲクッタもそうやが・・・
注意が必要たい。4次関数とかタチ悪いたいね。
※ちなみに、生物生存闘争モデル（例：ウサギとキツネの個体数の時間変化）のLotka-Volterra方程式も非線形の微分方程式です。
右辺のxyの部分が2次になっています。解曲線も極大を超えたら急激に減少するような挙動を示します。
[blogcard url=”http://shimaphoto03.com/science/rk-method/#Lotka”]
そっか・・・。4次関数は精度が荒いと2次関数と間違える可能性もあるよね・・。
そうそう。そういう危険も孕んでいるから、数値計算結果の過信はよくないんやが・・・。
100万回の計算結果と10万回の計算結果が同じなら、計算の信頼性が強くなるたいね。
例えばな、今回、n回に1回出力するプログラムが出来たから、
100万秒の現象を刻み幅1秒で計算して1万回に1回出力したときのグラフと
100万秒の現象を刻み幅10秒で計算して1,000回に1回出力したときのグラフを比較して一致すれば
刻み幅10秒で信頼たる計算が出来とると判断できるやろ？
一応言うが、支配方程式が間違っとるっていうのはナシやで。
確かに。両方ともプロットの数は100個だもんね。グラフのフォーマットも参照範囲も変えなくていいよね。
（そうか・・・プログラミングを自然現象の解析に使っている人はそういう事を考えるのか・・・・。）
いやぁ～。これは助かるよ！ありがとう！
いえいえ。いつも色々教えてもらってるから、こういうことくらいは・・・ね。
専門分野を異にするダチってええなぁ。
そうね。
Runge-kutta法で微分方程式を解いていて、
n回に1回出力できるようにならないかな～と思って、試行錯誤した結果をまとめました。
VBAではCなどとは違って、一工夫が必要で、Cと同じようにMODだけ使っていると、こうなります。
結果を見れば「そりゃ、そうか。」って感じなのですけどね。
VBAはあくまでツールで、それ専門じゃないので、またVBAでプログラムを組むのは数か月後ということは、あり得ることです。
そういった意味で、同じところで苦戦しないためにも備忘録を残しました。