資料請求番号：SH43 TS41
常微分方程式を使うと、様々な自然現象を記述することができます。
その常微分方程式の数値解法としてルンゲクッタ（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項は強制振動による力。外部からmF0cosωtだけ力が加わっていると考える。力を加えるんだから符号はプラスやな。
それで、両辺mで割って・・・
と、おけば、振動に関する2階の微分方程式が完成する。これをRunge-kutta法で解くたい。
あれ？でも・・・Runge-kuttaで解く方程式って、こんな形の1階常微分方程式じゃなかったかしら？
この方程式、Runge-kuttaで解けるの？
解ける。振動の運動方程式は、以下の2つの常微分方程式の連立と見ればいいだろう？
それで、2つの連立微分方程式はオマエ、こないだ実装したやないか。
あ、ホントだ。これならできそうね。
2つの方程式をt, x, vの関数としてみて・・・
こんな風にRunge-Kutta法を適用すれば出来そう！
おっし。じゃあ、振動の運動方程式を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——– */
ジャジャ～ン！！
（・・・かわいいな。）
とりあえずは、単振動、つまりω=γ=F0=0でやってみたよ。その方が検証しやすいでしょ？
ああ、そうだな。
上記の振動の運動方程式のうち、ω=γ=F0=0とし、抵抗力も外力もないとした場合、つまり
で表される運動を単振動と言います。
初期条件はx(0)=1, v(0)=0でやったのか。
うん。バネを引っ張った状態ね。そこから手を離したら、バネは余弦波の振動をするでしょ？
せや。確かにコサインの波が出来てるな。
それで、ω0を大きくしたら振動数は大きくなって波長が小さくなるところとか、
ω0をπにしたら、この曲線はcos(πt)になって1秒で半周、2秒で1周するところとかちゃんと再現できてるでしょ？
確かに。
ついでだからω0を0.1刻みで0~1まで変化させたときの動画を作ってみたよ！
へぇ～。こんなことができるんやなぁ～。現象の表現がより視覚的になってええな！
この動画の作り方のソースコードはまた今度まとめるね！
γに正の値を与えると、減衰しながら振動します。振動の減衰度合いが強すぎると振動しないまま運動が止まります。
減衰しながら振動する運動を減衰振動、振動せずに停止する運動を過減衰といいます。
減衰振動と過減衰の境目を臨界減衰といいます。臨界減衰を呈しているときγ＝ω0を満たします。
ω0を2にして、減衰係数を上げていく。γの値がω0を超えたら過減衰になるのよ。
この性質もちゃんと表現できてるでしょ？
確かにな。解析的に解くと、減衰振動・臨界減衰・過減衰とでは解曲線が変わってきてしまうからなぁ～。
うっかり、関数を切り替え忘れてて間違ったデータを出すこともあるし、そもそも切り替えがめんどくさかったりするから、数値解析は便利やな。
そうなのよ～。解析的に解けるとは言っても、その解析解が場合分けされているような現象だと、サッと答えを出しにくいのよね。
そういう意味でも、支配方程式を丸ごと数値解析してしまった方が早い場合もあるわよね。もちろん、精度が良ければの話だけどね。
こちらもγを0~3まで変化させたときの動画を作ってみたよ！
おおっ。減衰振動・臨界減衰・過減衰の違いがよく分かってええな！
本来は減衰して振動が止まってしまうところを、何らかの外力（強制振動）を加えることによって、振動を保たせ続けることができます。
本来、γ=0.1の減衰を受けて止まってしまう振動運動に対して・・・
例えば、F0=0.5 ω=1のF0cosωtという外力を加えれば、振動運動を続けることができるのよ。
なるほどな。最初は安定しないかもしれないけど、最終的に減衰と外力が釣り合って、平衡状態になるんやな。
これ、設計を上手くやれば最初から安定した振動が得られるかもな。
う～ん。言われてみればそうね。例えば、F0を強くしたりして0から4に変更すれば・・・・
あ、ホラ、F0=3がちょうどいいね！
おおっ～。これなら外部環境に合わせてどんな外力を加えればいいのかわかるから、振動を使った機械系の設計ができるな！
振動を使った機械系で怖いのが共振だが・・・共振は再現できるのか？
もちろん！今、ω0=2, γ=0.1, F0=3, ω=1で安定している系があるじゃない？
共振は自己の周波数ω0と外部周波数ωが一緒になった時に起きるのよね。
ふぅ～ん。良く知ってんじゃねぇか。
もう共振という危険な現象は、安全会議で耳タコになるくらい出てくるからね～・・・。実際、富士事業所にあるようなデッカイ機械が大きく揺れたらホント怖いし・・。
じゃあ、共振をシミュレーションしてみるね。ωをω0に近づけていくのよ。
おおっ、振幅が7倍になりよった！！
ダッシュポット（減衰力）と外力の比によってはもっと大変なことになりそうたいね。
例えば、γが1/10になってγ＝0.01になったとしたら・・・
※縦軸のスケールに注意
・・・25倍か～。怖いね。
怖い。ウチにはタービンとか攪拌とか回るものが特に多いから気をつけないといかんたいね・・・。
今回は振動の運動方程式をrunge-kutta法で解くための、エクセルのVBAプログラムを作成し、公開しました。
このプログラムを使用して、実際に単振動、減衰振動、強制振動の系について計算し、振動の運動の性質について説明しました。