みなさん，こんにちは
おかしょです．
積分をする方法にはさまざまなものがあります．その方法の一つにルンゲクッタと言うものがあります．
この記事を読むと以下のようなことがわかる・できるようになります．
- ルンゲクッタとは
- ルンゲクッタの計算方法
- MATLABのプログラムの書き方
この記事を読む前に
積分手法はルンゲクッタのほかにもさまざまなものがあります．ルンゲクッタは少し複雑な計算が必要なので，以下の記事で解説している台形法の方が理解はしやすいと思います．まだ読んでいない方は以下のリンクから読んでみてください．
ルンゲクッタとは
ルンゲクッタと一言で言っても，いくつか種類があります．その中でも必要十分な精度を得ることができる4次のルンゲクッタについて解説します．ルンゲクッタを簡単に説明すると，4種類の傾きを求め，それらの傾きから得られるデータの平均値を積分されたデータとして採用する手法のことです．
このようにルンゲクッタを言葉で説明してもわかりにくいので，図を使って解説していきます．
関数\((f(t,y))\)を時間\((t)\)でルンゲクッタを用いて積分することを考えます．
Step 1
まず一つ目の傾きは，被積分関数\((f(t,y))\)をそのまま利用します．これによって求められるデータは積分間隔\(dt\)を使って以下のように表されます．
\begin{eqnarray}
k_1 &=& f(t_i ,y_i) \times dt \tag{1}
\end{eqnarray}
上式を図で表すと以下のようになります．
上の図において水色の線は積分した結果を表しているとします．つまり求めたい解を表しています．時刻\(t_i\)におけるデータを\(y_i\)としています．このデータ\((t_i, y_i)\)を用いてルンゲクッタによって時刻\(t_{i+1}\)のデータ\(y_{i+1}\)を求めたいです．このstep 1では時刻\(t_i\)における傾きを被積分関数によって求めています．したがって，傾きに積分間隔\(dt\)をかけて求められるデータの変化量\(k_1\)を使って時刻\(t_{i+1}\)のデータを求めると上の図の赤い点のようになります．
この赤い点は明らかに真値（水色の線）とは異なる値を示しています．この\(k_1\)のみで求める積分手法を長方形則，もしくはオイラー法と言いこちらの記事で詳しく解説しています．
Step 2
Step 1の傾きだけでは精度が良くないので次の傾きを求めていきます．ここでは先ほど求めたデータの変化量\(k_1\)を用います．
Step 1では計算で用いるデータを\((t_i, y_i)\)としましたが，このデータを\((t_i+\frac{dt}{2}, y_i+\frac{k_1}{2})\)に移動します．つまり，積分間隔\(dt\)とデータの変化量\(k_1\)の半分ずつ移動することになります．
Step 2ではこの時のデータを被積分関数\((f(t,y))\)に代入して傾きを求めます．
\begin{eqnarray}
k_2 &=& f(t_i+\frac{dt}{2} ,y_i+\frac{k_1}{2}) \times dt \tag{2}
\end{eqnarray}
これを図で表すと以下のようになります．
緑色の線は移動した後のデータを用いて求められた傾きを表しています．この傾きを用いて求められるデータの変化量が\(k_2\)で，時刻\(t_{i+1}\)のデータが紫色の点で表されています．
Step 3
先ほどのstep 2と同様に，まずはデータを移動します．Step 3ではstep 2で求められたデータの変化量\(k_2\)を用いるので，以下のようになります．
\begin{eqnarray}
k_3 &=& f(t_i+\frac{dt}{2} ,y_i+\frac{k_2}{2}) \times dt \tag{3}
\end{eqnarray}
これも図で表すと以下のようになります．
オレンジ色の線は移動した後のデータを用いて求められた傾きを表しています．この傾きを用いて求められるデータの変化量が\(k_3\)で，時刻\(t_{i+1}\)のデータが青色の点で表されています．
Step 4
最後は使用するデータを\((t_i, y_i)\)から\((t_i+dt, y_i+k_3)\)に移動します．このデータを用いて求められる傾きを使ってデータの変化量\(k_4\)は以下のように求められます．
\begin{eqnarray}
k_4 &=& f(t_i+dt ,y_i+k_4) \times dt \tag{4}
\end{eqnarray}
これを図で表すと以下のようになる．
赤色の線は移動した後のデータを用いて求められた傾きを表しています．この傾きを用いて求められるデータの変化量が\(k_4\)で，時刻\(t_{i+1}\)のデータが水色の点で表されています．
Step 5
ルンゲクッタ法では最後にこれまでに求めたデータの変化量の平均を求めて，その平均値をデータの変化量として利用します．つまりデータの変化量\(k\)は以下のようになります．
\begin{eqnarray}
k &=& \frac{k_1+2(k_2+k_3)+k_4}{6} \tag{5}
\end{eqnarray}
ここで，上式は通常の平均ではないことに注意してください．
上の図を見るとわかるようにstep 1やstep 4で求められたデータよりもstep 2やstep 3で求められたデータの方が真値を表す水色の線に近いことがわかります．したがって，求めるデータの変化量\(k\)はstep 2とstep 3の時に求めた\(k_2\)と\(k_3\)を重要視すべきです．
そのため，式(5)ではstep 2とstep 3で求められたデータの変化量を2倍して平均値を求めています．このような平均を重み付き平均と呼びます．
この変化量\(k\)を用いて求められたデータこそが時刻\(t_{i+1}\)におけるデータ\(y_{i+1}\)になります．
\begin{eqnarray}
y &=& y_i+k\\
&=& y_i + \frac{k_1+2(k_2+k_3)+k_4}{6} \tag{6}
\end{eqnarray}
ルンゲクッタのプログラムの書き方
次式で表される関数を積分することを考えます．
\begin{eqnarray}
\dot{y} &=& \sin y \tag{1}
\end{eqnarray}
この関数を普通に計算で積分すると以下のようになります．
\begin{eqnarray}
\frac{dy}{dt} &=& \sin y\\
\int \frac{1}{\sin y} dy &=& \int dt\\
\int \frac{\sin y}{\sin^2 y} dy &=& \int dt\\
\int \frac{\sin y}{1-\cos^2 y} dy &=& \int dt\\
\int \frac{\sin y}{(1+\cos y)(1-\cos y)} dy &=& \int dt\\
\frac{1}{2} \int \left(\frac{\sin y}{1-\cos y}+\frac{\sin y}{1+\cos y}\right) dy &=& \int dt\\
\frac{1}{2} \left[ \ln (1-\cos y)-\ln (1+\cos y)\right]_{y_0}^y &=& t\\
\frac{1}{2} \left\{ \ln \left(\frac{1-\cos y}{1+\cos y}\right)-\ln \left( \frac{1-\cos y_0}{1+\cos y_0}\right)\right\} &=& t\\
\ln \left\{\frac{(1-\cos y)(1+\cos y_0)}{(1+\cos y)(1-\cos y_0)}\right\}^{\frac{1}{2}} &=& \ln e^t\\
\frac{(1-\cos y)(1+\cos y_0)}{(1+\cos y)(1-\cos y_0)} &=& e^{2t} \\
1-\cos y &=& (1+\cos y) \frac{1-\cos y_0}{1+\cos y_0}e^{2t}\\
&=& \frac{1-\cos y_0}{1+\cos y_0}e^{2t}+\cos y \frac{1-\cos y_0}{1+\cos y_0}e^{2t}\\
\cos y \left(1+\frac{1-\cos y_0}{1+\cos y_0}e^{2t}\right) &=& 1- \frac{1-\cos y_0}{1+\cos y_0}e^{2t}\\
\cos y &=& \frac{1- \frac{1-\cos y_0}{1+\cos y_0}e^{2t}}{1+\frac{1-\cos y_0}{1+\cos y_0}e^{2t}}\\
y &=& \cos^{-1} \left(\frac{1+\cos y_0-(1-\cos y_0) e^{2t}}{1+\cos y_0+(1-\cos y_0) e^{2t}}\right) \tag{2}
\end{eqnarray}
初期値\(y_0\)を0.05として図を描くと以下のようになります．
これをルンゲクッタで積分します．そのためのプログラムは以下のようになります．
% シミュレーション時間設定
time = 10;
dt = 0.01;
% データ格納庫
y = zeros(time/dt,1);
t = zeros(time/dt,1);
% 初期値設定
y(1) = 0.05;
for i=1:time/dt
% step1
state = y(i);
k1 = sin(state)*dt;
% step2
state = y(i)+1/2*k1;
k2 = sin(state)*dt;
% step3
state = y(i)+1/2*k2;
k3 = sin(state)*dt;
% step4
state = y(i)+k3;
k4 = sin(state)*dt;
% step5
state = (y(i)+(k1+2*(k2+k3)+k4)/6);
y(i+1) = state;
t(i+1) = t(i)+dt;
end
figure(1)
plot(t,y)
grid onこのプログラムをMATLABで実行すると以下のような図が表示されます．
先ほどの図と比較して，誤差をグラフ化すると以下のようになります．
このように誤差が非常に小さいことがわかります．
ルンゲクッタの精度を良くする方法
ルンゲクッタは非常に精度の良い積分手法ですが，1点だけ注意する必要があります．それは積分間隔です．当たり前のことですが，この積分間隔が広すぎると計算の精度は悪くなります．
たとえばロボットの運動を数値シミュレーションするために微分方程式をルンゲクッタで解く場合，積分間隔を広くとった場合と狭くとった場合で応答が全く異なることがあります．これはシステムの持つ固有振動数が深く関わっています．つまり共振現象をルンゲクッタで表現できるかということです．システムの固有振動数がわかっている場合は，それに合わせた積分間隔を設定すれば問題ありません．固有振動数がよくわからない場合は，何も考えずにとにかく小さい数字を設定しておけば大丈夫です．
まとめ
ここではルンゲクッタによる積分の方法を解説しました．この積分則を用いることで精密な数値シミュレーションを行うことができ，多くの方がこの方法で研究をしています．
続けて読む
このブログでは，この積分手法を用いて数値シミュレーションのやり方などを解説しています．以下の記事では数値シミュレーションを行う手順を紹介しているので参考にしてください．
Twitterでは記事の更新情報や活動の進捗などをつぶやいているので気が向いたらフォローしてください．
YouTubeでも制御工学や電子工作などの解説をしているので，ぜひのぞいてみてください．
それでは最後まで読んでいただきありがとうございました．
コメント