【ルンゲクッタ】連立の二階微分方程式をシミュレーションする方法をご紹介

私は別のページでコマのシミュレーションを紹介していますが、回るコマの動きを知るには結構複雑な連立の2次微分方程式を解く必要がありました。

微分方程式をシミュレーションで解こうとした際、「ルンゲクッタ」という手法をよく耳にしますが、連立で高次（２次）の微分方程式に対してどのようにすればルンゲクッタを利用できるかはじめはわかりませんでした。

しかし、少し調べて取り組んでみると、驚くほど簡単に連立の2次微分方程式をルンゲクッタで解くことができたので、このページではプログラミングに落とし込めるよう具体的にご紹介します。

あわせて読みたい

【倒れない】コマのシミュレーションは意外に簡単だった なぜ回っているコマは倒れないの？という疑問は多くの方が持っていると思います。 私が「コマが倒れない理由」でネット検索しても、「ジャイロ効果があるから」「歳差運…

ルンゲクッタは微分方程式を精度良くシミュレーションできる

まず、一番簡単な例として微分方程式$\frac{\text{d}y}{\text{d}x}=xy$をシミュレーションしてみましょう。（ただし、$y(0)=4$とします。）

なお、微分方程式の理論解は$y=4\exp(x^2/2)$です。

オイラー法で解いてみる

この微分方程式をオイラー法で解いてみましょう。

オイラー法のアルゴリズムはある時刻$i$での$x$の値、$x_i$だけを使って$x$を時間発展させます。微分方程式$\frac{\text{d}y}{\text{d}x}=xy$を離散化させると、$\frac{y_{i+1}-y_i}{\Delta{x}}=x_i\times{y_i}$なので、$y_{i+1}=y_i+{x_i}{y_i}\Delta{x}$とできます。

この式を使うと、例えば$\Delta{x}=0.02$ととると、$y_2=4+0*4*0.02=4$、$y_3=4+0.02*4*0.02=4.0016$と順に時間発展させることができます。

オイラー法は簡単にシミュレーションできますが、オイラー法では、$x$が大きい領域で解析解との差が大きくなってしまい、精度が高くありません。

そこで、登場するのがルンゲクッタの方法です。

ルンゲクッタの方法で解いてみる

では、この微分方程式をルンゲクッタの方法で解いてみましょう。

ルンゲクッタの方法は簡単に言うと、「未来の傾きを先読みする」シミュレーションです。

ルンゲクッタの方法では、$y$を$y_{i+1}=y_i+(k_1+2k_2+2k_3+k_4)*\Delta{x}$で時間発展させます。この式において、$k_1$は$(x,y)=(x_i,y_i)$での傾き、$k_2$は$(x,y)=(x_i+\frac{\Delta{x}}{2},y_i+\frac{k_1\Delta{x}}{2})$での傾き、$k_3$は$(x,y)=(x_i+\frac{\Delta{x}}{2},y_i+\frac{k_2\Delta{x}}{2})$での傾き、$k_4$は$(x,y)=(x_i+{\Delta{x}},y_i+k_3\Delta{x})$での傾きです。

つまり、微分方程式$\frac{\text{d}y}{\text{d}x}=xy$の場合、$k_1=x_i{y_i}$、$k_2=(x_i+\frac{\Delta{x}}{2})(y_i+\frac{k_1\Delta{x}}{2})$、$k_3=(x_i+\frac{\Delta{x}}{2})(y_i+\frac{k_2\Delta{x}}{2})$、$k_4=(x_i+\frac{\Delta{x}}{2})(y_i+k_3\Delta{x})$になります。

オイラー法に比べるとややこしいね

ややこしい分、オイラー法より精度がかなり上がります

ルンゲクッタ法は$(x_i,y_i)$から$(x_{i+1},y_{i+1})$の値に更新する際、$x=x_i$での傾きだけでなく、$x=x_i+\frac{\Delta{x}}{2}$などの値も使うため、オイラー法よりもかなり精度が上がります。

高次の微分方程式の場合

難易度を少し上げて、次は二階微分方程式$\frac{\text{d}^2{x}}{\text{d}t^2}+\omega_0^2{x}-f_0\cos\omega_ft=0$をルンゲクッタの方法で解いてみましょう。

ルンゲクッタの方法を二階微分方程式に適用するためには、少し工夫が必要です。具体的には、この式を$\frac{\text{d}{x}}{\text{d}{t}}=y$と置くことで、（無理やり）二元一階微分方程式に書き換えます。

$ \left\{ \begin{array}{l} \frac{\text{d}{x}}{\text{d}{t}}&=y\\ \frac{\text{d}{y}}{\text{d}{t}}&=-\omega_0^2x+f_0\cos\omega_ft\end{array} \right. $

このようにすることで、二階微分方程式であっても一階微分方程式であった$\frac{\text{d}y}{\text{d}x}=xy$の時と同じようにして解くことができます。

微分方程式を書き換えた後は、ルンゲクッタは次のような計算式によって値を更新していくことができます。

$x_{i+1}=x_i+\frac{(k_{x1}+2k_{x2}+2k_{x3}+k_{x4})}{6}$、$y_{i+1}=y_i+\frac{(k_{y1}+2k_{y2}+2k_{y3}+k_{y4})}{6}$

$k_{x1}=y_i\Delta{t}$、$k_{y1}=\{-\omega_0^2x_i+f_0\cos\omega_ft_i\}\Delta{t}$

$k_{x2}=(y_i+\frac{k_{y1}}{2})\Delta{t}$、$k_{y2}=\{-\omega_0^2(x_i+\frac{k_{x1}}{2})+f_0\cos\omega_f(t_i+\frac{\Delta{t}}{2})\}\Delta{t}$

$k_{x3}=(y_i+\frac{k_{y2}}{2})\Delta{t}$、$k_{y3}=\{-\omega_0^2(x_i+\frac{k_{x2}}{2})+f_0\cos\omega_f(t_i+\frac{\Delta{t}}{2})\}\Delta{t}$

$k_{x4}=(y_i+{k_{y3}})\Delta{t}$、$k_{y3}=\{-\omega_0^2(x_i+{k_{x3}})+f_0\cos\omega_f(t_i+{\Delta{t}})\})\Delta{t}$

定数を具体的な値に置き換えて実際に解いてみましょう。$\omega_0=2$、$\omega_f=1$、$f_0=3$とすると、二階微分方程式$\frac{\text{d}^2{x}}{\text{d}t^2}+4{x}-3\cos{t}=0$の解析解は$x(t)=\sqrt{2}\cos(2t-\frac{\pi}{4})+\cos(t)$で与えられます。

ルンゲクッタの方は、上の計算式に$\omega_0=2$、$\omega_f=1$、$f_0=3$をそれぞれ代入したものをpythonでプログラミングしてシミュレーションしました。（時間ステップは$\Delta{t}=0.02$を利用）

import math
import numpy as np
import numpy.linalg as LA
np.set_printoptions(precision=6,floatmode='fixed')

omega_0  = 2
omega_f  = 1
f_0      = 3
Dt   = 0.02

x  = 2
y  = 2

f = open('x.txt','w')


ite = 0

for ite in range(400):

 t = Dt * ite

 k_x1     =  y  * Dt
 k_y1     = (-(omega_0**2) * x + f_0*math.cos(omega_f*t)) * Dt
 
 k_x2     =  (y+k_y1/2)  * Dt
 k_y2     = (-(omega_0**2) * (x+k_x1/2) + f_0*math.cos(omega_f*(t+Dt/2))) * Dt

 k_x3     =  (y+k_y2/2)  * Dt
 k_y3     = (-(omega_0**2) * (x+k_x2/2) + f_0*math.cos(omega_f*(t+Dt/2))) * Dt

 k_x4     =  (y+k_y3)  * Dt
 k_y4     = (-(omega_0**2) * (x+k_x3) + f_0*math.cos(omega_f*(t+Dt))) * Dt

 x   =   x + (k_x1 + 2*k_x2  + 2*k_x3  + k_x4)/6
 y   =   y + (k_y1 + 2*k_y2  + 2*k_y3  + k_y4)/6


 if (ite % 1) == 0:
    print(ite,t,x,y,file=f)

f.close()

ルンゲクッタの結果を見ると、解析解と非常に近い計算結果を与えていることが分かります。

このように、ルンゲクッタでは二階微分方程式も問題なくシミュレーションができます。

非常に複雑な微分方程式の場合

最後に、「非常に複雑な微分方程式」として、回るコマの動きをシミュレーションするための微分方程式をルンゲクッタで解いてみましょう。

回るコマの運動方程式は次のようになります。

最後にふさわしいくらい複雑な微分方程式だね。。

$\begin{eqnarray}
\left\{
\begin{array}{l}
A\ddot{\theta}-A\dot{\phi}^2\sin\theta\cos\theta+Cn\dot{\phi}\sin\theta=Mgh\sin\theta\\
A\ddot{\phi}\sin\theta+2A\dot{\phi}\dot{\theta}\cos\theta-Cn\dot{\theta}=0\\
\dot{\phi}\cos\theta+\dot{\psi}=n
\end{array}
\right.
\end{eqnarray}$

あわせて読みたい

【倒れない】コマのシミュレーションは意外に簡単だった なぜ回っているコマは倒れないの？という疑問は多くの方が持っていると思います。 私が「コマが倒れない理由」でネット検索しても、「ジャイロ効果があるから」「歳差運…

複雑なこの微分方程式に関しても、前と同じように、$\dot{\theta}=\eta$、$\dot{\phi}=\epsilon$と置くことで、二階微分の部分を一階微分に書き換えてみましょう。

$\begin{eqnarray}
\left\{
\begin{array}{l}
\dot{\theta}=\eta\\
\dot{\phi}=\epsilon\\
\dot{\eta}=\frac{(-A{\epsilon}^2\cos\theta-Cn\dot{\phi}+Mgh)\sin\theta}{A}\\
\dot{\epsilon}=\frac{2A{\eta}{\epsilon}\cos\theta+Cn{\eta}}{A\sin\theta}\\
\dot{\psi}=n-{\epsilon}\cos\theta
\end{array}
\right.
\end{eqnarray}$

このように、一階の微分方程式に書き直すことで、これまでと同じようにルンゲクッタの方法が使えるようになります。

$\begin{eqnarray}
\left\{
\begin{array}{l}
k_{\theta1}=\eta_i*\Delta{t}\\
k_{\phi1}=\epsilon_i*\Delta{t}\\
k_{\eta1}=\frac{\{-A\epsilon_i^2\cos\theta_i-Cn\epsilon_i+Mgh\}\sin\theta_i}{A}*\Delta{t}\\
k_{\epsilon1}=\frac{-2A\eta_i\epsilon_i\cos\theta_i+Cn\eta_i}{A\sin\theta_i}*\Delta{t}\\
k_{\psi1}=(n-\epsilon_i\cos\theta_i)*\Delta{t}\\
k_{\theta2}=(\eta_i+\frac{k_{\theta1}}{2})*\Delta{t}\\
k_{\phi2}=(\epsilon_i+\frac{k_{\epsilon1}}{2})*\Delta{t}\\
k_{\eta2}=\frac{\{-A(\epsilon_i+\frac{k_{\epsilon2}}{2})^2\cos(\theta_i+\frac{k_{\theta1}}{2})-Cn(\epsilon_i+\frac{k_{\epsilon1}}{2})+Mgh\}\sin(\theta_i+\frac{k_{\theta1}}{2})}{A}*\Delta{t}\\
\vdots\\
k_{\epsilon3}=\frac{-2A(\eta_i+\frac{k_{\eta2}}{2})(\epsilon_i+\frac{k_{\epsilon2}}{2})\cos(\theta_i+\frac{k_{\theta2}}{2})+Cn(\eta_i+\frac{k_{\eta2}}{2})}{A(\sin\theta_i+\frac{k_{\theta2}}{2})}*\Delta{t}\\
\vdots\\
k_{\psi4}=\{n-(\epsilon_i+\frac{k_{\epsilon3}}{2})\cos(\theta_i+\frac{k_{\theta3}}{2})\}*\Delta{t}
\end{array}
\right.
\end{eqnarray}$

これらのように、値を作ることで、これまでと同じように値を更新することができます。

$\begin{eqnarray}
\left\{
\begin{array}{l}
\theta_{i+1}=\theta+\frac{1}{6}(k_{\theta1}+2k_{\theta2}+2k_{\theta3}+k_{\theta4})\\
\phi_{i+1}=\phi+\frac{1}{6}(k_{\phi1}+2k_{\phi2}+2k_{\phi3}+k_{\phi4})\\
\vdots
\end{array}
\right.
\end{eqnarray}$

この数式をpythonでプログラミングすることで、回るコマが歳差運動しながら倒れないシミュレーションを行うことができます。

プログラミングコードはこちらのページをご参照ください

まとめ

このページでは、微分方程式をシミュレーションする際によく用いられる「ルンゲクッタの方法」について、数値計算式が具体的に分かるように解説しました。

連立で高次（２次）の微分方程式に対して、どのようにすればルンゲクッタを利用できるか分かりにくいかもしれませんが、実は驚くほど簡単にプログラミングに落とし込むことができます。

ただし、微分方程式をいきなりルンゲクッタの方法でプログラミングすると間違いを起こしやすいので、オイラーの方法ルンゲクッタの方法でのシミュレーションをそれぞれエクセルで行い挙動を把握した後、ルンゲクッタの方法をプログラミングすることをお勧めします。