資料請求番号：SH43 SH44 TS41
CONTENTS
大学や会社で研究開発活動をしていると、微分方程式における様々なパラメータが変化すると解曲線はどのように変化するのかをストレスなく確認したいという場面があると思います。
また、パラメータを求める実験を行う際には、ターゲットとしている微分方程式がパラメータによってどのように変動するのかを予め知っておけば、実験計画をスムーズに立てられると思います。
今回は、微分方程式のパラメータの変化が解曲線にどのような影響を与えるのかを理解するためのアニメーションの作成方法をまとめたいと思います。
今回作成したいアニメーションは以下の通りです。こちらは振動現象の微分方程式
をルンゲ・クッタ法で解いた曲線で、「γ=0.1の粘性によって減衰してしまうところ、どのくらいの強制振動F0を与えると振動を安定に続けることができるのか」を考察しています。
※こちらのページは振動に関する上記の微分方程式をルンゲクッタ法で解くプログラムがあることを前提にしています。ルンゲクッタ法のプログラムは以下の記事で解説しています。
| 資料請求番号：SH43 TS41Runge-Kutta法を使った振動計算常微分方程式を使うと、様々な自然現象を記述することができます。その常微分方程式の数値解法としてルンゲクッタ（Runge-kutta）法があります。※Runge-kutta法の詳しい解説はこちら今回は、振動運動の運動方程... 振動の運動方程式をエクセルVBAを使って解く - らい・ぶらり |
/* –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 properties
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)
‘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——– */
なんでそんなケンカ腰なのよ！？イライラしてるからってワタシに当たらないでちょうだい！！
まぁ、確かにないけど。ワタシはプログラム書くことが仕事だら？一応微分方程式を解くプログラムは書けるけど、その微分方程式の中身はよく分かんない～なんてよくあることだし？
ただね、例えば・・・この微分方程式のF0ってパラメータあるでしょ？
ああ。
これを0.1ずつ変動させながら、繰り返しグラフを描画させられる方法なら知ってるわよ。
おおっ、ってことは、いちいち0.1ずつ変化させてルンゲクッタループしなくてもええってこと？
まぁね。1回アニメーション見て、気になったパラメータで刻み幅を小さくするなり、他のパラメータ考察するなりすれば、効率よさそうじゃない？
せやな！教えてくれよ！
お腹すいたな～。今日、テンイチ行きたいな～
・・・しょうがねぇな。ええよ。行こう。
やった！
それで、今回やりたいことっていうのは、例えばF0を0.5から4まで動かしたときの解曲線の挙動が見たいっていう目的なら、
0.5でルンゲクッタループ回して、グラフ見て、0.6でルンゲクッタループ回して・・・を続けて4までたどり着くようにすればいいのよね。
せやな。
なら、基本形として、C7セルにあるF0の値を代入しては計算させてプロットさせての繰り返しになるわね。
/* –vibration constant simulation——– */
For force = 0.5 To 1 Step 0.1
Cells(7, 3) = force
Application.Run “rkvib”
Next force
/* –Rose curve Realtime plot ——End– */
このApplication.Run “rkvib”ってのは・・・ルンゲクッタループか？
/* –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 properties
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)
‘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——– */
そうそう！0.5でルンゲクッタループ回して、Step0.1だから次は0.6で回して・・ってこと。
ほう。繰り返しの制御変数に整数型以外も使えるのか。便利やな。
まぁ、一応、テストとして、とりあえず1までにしてやってみたら？
わかった。
・・・・。
数字が変化するだけで、リアルタイムでグラフが変化せんのやけど。
もうテンイチ行かんぞ！！
まぁまぁ待って待って。
あのね、VBAのプログラム回している間はExcelの制御を完全にVBAに渡している状態になるの。
だから、時間のかかるマクロを回している間はExcelがフリーズしたみたいな感じになるのよね。
で？
それで一旦制御をOSに返しますっていう命令を書かなきゃいけないの。それがDo Events関数。
これを入れてやってみなよ。制御がOSに返ったらグラフを描くことができるから。
/* –vibration constant simulation——– */
For force = 0.5 To 1 Step 0.1
Cells(7, 3) = force
Application.Run “rkvib”
DoEvents: DoEvents: DoEvents
Next force
/* –Rose curve Realtime plot ——End– */
※Sub rkvib()
/* –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 properties
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)
‘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——– */
おおっ！グラフがうごきよった！！最初からこれ教えろよ～（笑）
でしょ？スゴイでしょ？そしたら、F0 = 4までまわしてごらん？
これでF0=3がちょうどいいってことがわかるし、F0 = 2.1~2.2で重ね合わせ特有の不安定な振動がなくなるのもわかるたいね！
いやぁ～。よかよか！ありがとう！
じゃ、テンイチ行きましょ！
はいはい・・・。
View Comments
こんにちは、
貴殿のサイト読ましてもらっています。
とても詳しくわかりやすく説明されていて、VBA初心者の私にとってはとても学びやすいです。
ところで、以下のコードなんですが、これをどこに書けば最後の動画のようになるのでしょうか？
rkvib内のどこかなのかモジュール（Module1 やModule2）内に導入すればいいのでしょうか？
For force = 0.5 To 1 Step 0.1
Cells(7, 3) = force
Application.Run “rkvib”
DoEvents: DoEvents: DoEvents
Next force
こんにちは。
質問の回答ですが、同一モジュール内にrkvibとは別の新しいSubプロシージャを作成して
For force = 0.5 To 1 Step 0.1
Cells(7, 3) = force
Application.Run “rkvib”
DoEvents: DoEvents: DoEvents
Next force
の構文を書けばよろしいかと思います。