【Excel VBA解説】振動方程式をRunge-Kutta法で数値解析！減衰振動・共振まで再現

Runge-Kutta法を使った振動計算

常微分方程式を使うと、様々な自然現象を記述することができます。
その常微分方程式の数値解法としてルンゲクッタ（Runge-kutta）法があります。

※Runge-kutta法の詳しい解説はこちら
[blogcard url=”https://shimaphoto03.com/science/rk-method/”]
今回は、振動運動の運動方程式をRunge-kutta法で解くためのプログラムをExcelのVBAで作成しました。
その備忘録として、本記事にソースコードを書き残し、計算結果を実現象と照らし合わせてみました。

※ソースコードはこちら

振動の運動方程式

今回、Runge-kutta法で解く方程式は以下の通りだ。振動の運動方程式。

xを変位とすればxの2回微分で加速度になるから、この式はma=Fの運動方程式の形をしているのね。

せや。右辺第一項はフックの法則を表している。バネを引っ張れば、それに対して戻ろうとする力が働くだろう？
だから-kxになってるたい。

第2項はこれ・・・抵抗力かしら？なんでxの微分がついているのかしら？

確かに抵抗力なんやが、これは・・・そうやな・・・
水の中で手をゆっくり動かしたときと早く動かしたときとどっちが手を前に進めにくいか？で考えれば
抵抗力は運動の速度に比例するってことが分かると思うたい。

なるほど・・・。

第3項は強制振動による力。外部からmF₀cosωtだけ力が加わっていると考える。力を加えるんだから符号はプラスやな。

それで、両辺mで割って・・・

と、おけば、振動に関する2階の微分方程式が完成する。これをRunge-kutta法で解くたい。

振動の運動方程式（2階微分方程式）をRunge-Kutta法で解く

2階微分方程式→2個の連立方程式

あれ？でも・・・Runge-kuttaで解く方程式って、こんな形の1階常微分方程式じゃなかったかしら？

この方程式、Runge-kuttaで解けるの？

解ける。振動の運動方程式は、以下の2つの常微分方程式の連立と見ればいいだろう？
それで、2つの連立微分方程式はオマエ、こないだ実装したやないか。

あ、ホントだ。これならできそうね。
2つの方程式をt, x, vの関数としてみて・・・

こんな風にRunge-Kutta法を適用すれば出来そう！

おっし。じゃあ、振動の運動方程式をRunge-kuttaで解くマクロを作ってみろよ。

うん。すぐにできるわよ。前回のソースコードをこの方程式に変えればいいだけでしょう？

[blogcard url=”https://shimaphoto03.com/science/vba-rk/#code”]

・・・・

できた。

おお、早いな。じゃあ、検証していこうか。

ソースコードと検証

ソースコード

こんな感じのワークシートで・・・

このソースコードを走らせると・・・

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

ジャジャ～ン！！

（・・・かわいいな。）

とりあえずは、単振動、つまりω=γ=F₀=0でやってみたよ。その方が検証しやすいでしょ？

ああ、そうだな。

単振動

上記の振動の運動方程式のうち、ω=γ=F₀=0とし、抵抗力も外力もないとした場合、つまり

で表される運動を単振動と言います。

初期条件はx(0)=1, v(0)=0でやったのか。

うん。バネを引っ張った状態ね。そこから手を離したら、バネは余弦波の振動をするでしょ？

せや。確かにコサインの波が出来てるな。

それで、ω₀を大きくしたら振動数は大きくなって波長が小さくなるところとか、
ω₀をπにしたら、この曲線はcos(πt)になって1秒で半周、2秒で1周するところとかちゃんと再現できてるでしょ？

確かに。

ついでだからω₀を0.1刻みで0~1まで変化させたときの動画を作ってみたよ！

へぇ～。こんなことができるんやなぁ～。現象の表現がより視覚的になってええな！

この動画の作り方のソースコードはまた今度まとめるね！

減衰振動

γに正の値を与えると、減衰しながら振動します。振動の減衰度合いが強すぎると振動しないまま運動が止まります。
減衰しながら振動する運動を減衰振動、振動せずに停止する運動を過減衰といいます。
減衰振動と過減衰の境目を臨界減衰といいます。臨界減衰を呈しているときγ＝ω₀を満たします。