資料請求番号：SH43 TS41

CONTENTS

1 振動の運動方程式
2 振動の運動方程式（2階微分方程式）をRunge-Kutta法で解く
- 2.1 2階微分方程式→2個の連立方程式
3 ソースコードと検証
4 まとめ

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

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

※Runge-kutta法の詳しい解説はこちら

らい・ぶらり

数値計算を使って常微分方程式を解く～ルンゲクッタ法の解説～

http://shimaphoto03.com/science/rk-method/

資料請求番号：TS33 TS38数値計算による微分方程式解法の基本～runge-kutta法～ルンゲ・クッタ（runge-kutta）法は、コンピュータを使った微分方程式解法の基本として工学部の基礎教育に登場します。中でも4次のルンゲ・クッタ法はアルゴリズムが簡単な上に高い精度...

今回は、振動運動の運動方程式を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で解くマクロを作ってみろよ。

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

らい・ぶらり

エクセルVBAを用いた微分方程式の計算（ルンゲクッタ法）

http://shimaphoto03.com/science/vba-rk/#code

資料請求番号：SH43Excel VBAマクロでRunge-kutta法を使い、常微分方程式を解くルンゲクッタ法とは、数値計算による常微分方程式の解法の一つであり、特に4次精度Runge-kutta法はアルゴリズムが簡便であり、精度も高いことから、工学者から高い支持を受けており、現...

・・・・

できた。

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

ソースコードと検証

ソースコード

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

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

/* –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まで変化させたときの動画を作ってみたよ！

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

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

減衰振動

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