Categories: PC/プログラミング化学/物理

【Excel VBA】微分方程式パラメータを動的変更！リアルタイム描画で挙動を迅速確認





【Excel VBA】微分方程式パラメータを動的変更！リアルタイム描画で挙動を迅速確認



資料請求番号：SH43 SH44 TS41

微分方程式パラメータが解曲線に与える影響を考察するために

大学や会社で研究開発活動をしていると、微分方程式における様々なパラメータが変化すると解曲線はどのように変化するのかをストレスなく確認したいという場面があると思います。
また、パラメータを求める実験を行う際には、ターゲットとしている微分方程式がパラメータによってどのように変動するのかを予め知っておけば、実験計画をスムーズに立てられると思います。

今回は、微分方程式のパラメータの変化が解曲線にどのような影響を与えるのかを理解するためのアニメーションの作成方法をまとめたいと思います。

今回作成したいアニメーションは以下の通りです。こちらは振動現象の微分方程式

をルンゲ・クッタ法で解いた曲線で、「γ=0.1の粘性によって減衰してしまうところ、どのくらいの強制振動F0を与えると振動を安定に続けることができるのか」を考察しています。

※こちらのページは振動に関する上記の微分方程式をルンゲクッタ法で解くプログラムがあることを前提にしています。ルンゲクッタ法のプログラムは以下の記事で解説しています。
[blogcard url=”https://shimaphoto03.com/program/rk-vib/”] ※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”ってのは・・・ルンゲクッタループか？

※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——– */

シママ

そうそう！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で重ね合わせ特有の不安定な振動がなくなるのもわかるたいね！

ストーク

いやぁ～。よかよか！ありがとう！

シママ

じゃ、テンイチ行きましょ！

ストーク

はいはい・・・。

shimakei8364

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
      の構文を書けばよろしいかと思います。

Leave a Comment
Share
Published by
shimakei8364
6年 ago

Recent Posts

【大分観光モデルコース】2泊3日でめぐる絶景と温泉満喫ドライブプラン

【大分観光モデルコース】2泊3…

1週間 ago

Pythonで始めるPID制御シミュレーション：2槽タンクモデルで直感的理解

Pythonで始めるPID制御…

3週間 ago

移動現象論とは？物質・熱・運動量が「多いところから少ないところへ」移動するしくみを身近な例で解説

移動現象論とは？物質・熱・運動…

3週間 ago

【高校物理】ドップラー効果の公式と図解による簡潔な解法～場合分け不要～

【高校物理】ドップラー効果の公…

3週間 ago

【ラプラス変換】RLC回路・質量ばねダンパ系の複雑な微分方程式を解く方法

【ラプラス変換】RLC回路・質…

3週間 ago

浦島太郎の起源と変遷：御伽草子版と現代版に見る「異類婚姻譚」と「善悪応報譚」のギャップ

浦島太郎の起源と変遷：御伽草子…

3週間 ago