4次ルンゲ・クッタ法

微分方程式なら任せろーバリバリバリー

ルンゲ=クッタ法ってもともとどういうもの?理論は?刻み幅\(h\)を自動的に制御する方法について知りたい方は、ルンゲ=クッタ法の系統的扱いと刻み幅制御へどうぞ。

4次ルンゲ=クッタ法は微分方程式を数値的に解く手段です。

ルンゲ=クッタ法が良く使われる理由は、ひとえにプログラムの実装のしやすさです。ルンゲクッタ法は非常に良い!というアルゴリズムではありませんが、他の方法よりもシンプルで、プログラムに組み込みやすいのです。

”4次”は問題の解の関数をテイラー展開した場合、4次までは一致するように作られた、という意味です。
例えば微分方程式
\(
\displaystyle \frac{dy}{dx}=g(x,y)
\)

を考えます。
数値計算では初期値(例えばx=0の時、y=1など)を決めて、そこからxをhだけ増やし、微分方程式というルールに従って関数\(y(0+h),y(0+h+h),y(0+h+h+h),\cdots\)を作り上げていきます。
この時、4次ルンゲ=クッタ法で求められる答えの関数というのは
\(
\displaystyle y(x+h)=y(x)+hg(x,y)+\frac{h^2}{2!}\frac{dg}{dx}+\frac{h^3}{3!}\frac{d^2g}{dx^2}+\frac{h^4}{4!}\frac{d^3g}{dx^3}+O(h^5)
\)

という関数になります。
\(
g(x,y)=-xy
\)

(ただし、初期条件\(x=0\)で\(y=1\))
である場合、微分方程式の解析解は
\(
\displaystyle y=e^{-x^2/2}
\)

であるため、4次ルンゲ=クッタ法によって導かれる答えは、
\(
\displaystyle y(x+h)=y(x)\left[1-hx+\frac{h^2}{2}(x^2-1)+\frac{h^3}{6}(-x^3+3x)+\frac{h^4}{24}(x^4-6x^2+3)\right]+O(h^5)
\)

となります。

本題に入りましょう。
4次ルンゲ=クッタ法は6つのステップが必要となります。
初期値を\((x_0,y_0)\)と書くと、(\(y_0=y(x_0)\)です。)

  1. \((x_0,y_0)\)より\(k_a\)を求める
  2. \(k_a\)より\(k_b\)を求める
  3. \(k_b\)より\(k_c\)を求める
  4. \(k_c\)より\(k_d\)を求める
  5. \(k_a,k_b,k_c,k_d\)より\(y(x_0+h)\)を求める
  6. \((x_0+h,y(x_0+h))\)を初期値だと思って手順1に戻る。

という感じです。
3章:連立ルンゲ・クッタ法による微分方程式の解を参考にすると、数値計算での4次ルンゲ=クッタ法の計算スキームは以下のようになります。

解きたい微分方程式を連立1次微分方程式の形で書くと一般的にはこう書けます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy_1}{dx} &= f_1(x,y_1,y_2,\cdots,y_N) \\
\frac{dy_2}{dx} &= f_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
\frac{dy_N}{dx} &= f_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
まず係数\(k_a\)を求めます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{a1} &= hf_1(x,y_1,y_2,\cdots,y_N) \\
k_{a2} &= hf_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
k_{aN} &= hf_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
次に\(k_b\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{b1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
k_{b2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
&\vdots \\
k_{bN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
そして\(k_c\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{c1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
k_{c2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
&\vdots \\
k_{cN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に\(k_d\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{d1} &= hf_1(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
k_{d2} &= hf_2(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
&\vdots \\
k_{dN} &= hf_N(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に求めた\(k_a,k_b,k_c,k_d\)を使って\(x+h\)でのそれぞれの関数の値を導きます。、
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
x &= x+h \\
y_{1} &= y_{1}+(k_{a1}+2k_{b1}+2k_{c1}+k_{d1})/6 \\
y_{2} &= y_{2}+(k_{a2}+2k_{b2}+2k_{c2}+k_{d2})/6 \\
&\vdots \\
y_{N} &= y_{N}+(k_{aN}+2k_{bN}+2k_{cN}+k_{dN})/6 \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
となります。

4次ルンゲ=クッタ法のプログラム


実際に例題を解きましょう。
2つの例題を解きます。

[adsense1]

1階微分方程式


1つ目は
\(
\displaystyle \frac{dy}{dx}=-y\sin{x}
\)

です。
一般解は
\(
y=Ae^{\cos{x}}
\)

であり、初期条件\(x=0,y=2\)として解けば解析解は
\(
y=2e^{\cos{x}-1}
\)

です。これを4次ルンゲ=クッタ法を用いて解くには、以下のプログラムで実現できます。

解いて、gnuplotで
plot “fort.10”
数値計算解(赤点)解析解(緑線)と共に出力すると、こんなグラフが得られます。
4次ルンゲクッタ法による計算1

連立1階微分方程式


もう一つの例題は微分方程式
\(
\displaystyle \frac{d^2y}{dx^2}+2\gamma \frac{dy}{dx}+y=0
\)

を解きます。これは物理の世界ではバネの減衰振動を表す運動方程式で(詳しくは減衰振動へ。)、\(0\lt\gamma\lt1\)の場合、解は
\(
y(x)=Ae^{-\gamma x}\cos(x\sqrt{1-\gamma^2}-\alpha)
\)

となります。
\(\gamma=0.15\), 初期条件を今、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\left.y\right|_{x=0} &= 1 \\
\left.\frac{dy}{dx}\right|_{x=0} &= -0.15 \\
\end{aligned}
\right.
\end{eqnarray}
\)
とすると、解析解は
\(
y(x)=e^{-\gamma x}\cos{(x\sqrt{1-0.15^2})}
\)

となります。
プログラムで計算する際は、まず連立1次微分方程式に焼きなおす必要があります。すなわち、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy}{dx}&=v \\
\frac{dv}{dx}&=-2\gamma v -y \\
\end{aligned}
\right.
\end{eqnarray}
\)
として解けばいいわけです。

その時の変更するべき場所は少しで、mainとrkfdの一部を変えればおしまいです。
変えた場所をコメントで
!+—+
を入れたので確認してください。

gnuplotで数値計算解(赤点)解析解(緑線)のグラフを書くとこうでなります。
4次ルンゲクッタ法による計算2

[adsense2]


「4次ルンゲ・クッタ法」への8件のフィードバック

  1. function関数のところのsは何を表しているのでしょうか。是非教えてください。

    1. sは微分方程式のインデックスを表しています。
      すなわち、
      s=1の時、1本目の微分方程式の右辺、
      s=2の時、2本目の微分方程式の右辺、
      ・・・
      s=s’の時、s’本目の微分方程式の右辺、
      ・・・
      を返す、という意味です。
      二階微分方程式の例をご覧になると良く分かるかもしれません。

      ※sによる微分方程式の区別は綺麗な方法ではないので、実用には向かないと思います。その点はご了承下さい。

  2. 1階微分方程式を同じコードで解かせてみたのですが、誤差が2次精度で出てきました。4次精度のルンゲクッタを使っているのに2次精度までしか出ないのはなぜなのでしょうか。。

  3. すみません、私の勘違いで、4次精度が出ていました。ありがとうございます。

  4. dy/dx がxに露わによらない場合はfの中の(x+h/2)は無視して良いのでしょうか?

    具体的には物理でよくx(t)を扱う場面があるのですが、実際dx/xtが陽にtに依ることは少ないです。しかしx-t図がかける以上、(t + h/2 , x + k1/2)という座標がかけるので、無視するのが気持ち悪くて…

    1. 無視して良いです。
      例えば
      \begin{align}
      \frac{dy}{dx}=f(x,y)=y
      \end{align}
      の場合、
      \begin{align}
      f(x+\frac{h}{2}, y+\frac{k_1}{2})=y+\frac{k_1}{2}
      \end{align}
      となるので無いものとしてみて良いです。

      > しかしx-t図がかける以上、(t + h/2 , x + k1/2)という座標がかけるので、
      という点については媒介変数の考えと同じです。媒介変数は多くの問題で消してしまいますが、媒介変数を横軸にとって関数を書くことはできます。
      我々の世界では時刻tが媒介変数であり、時間の進み方を勝手に変えない、と定義しているからこそ、時刻は意味のある媒介変数となっています。すなわち時刻は意味を持たなければならない、という先入観があるからこそ
      >無視するのが気持ち悪くて…
      につながっているのだと思います。

      時刻が固定されているからこそ、運動方程式では”時間”を単位に含む物理量に、媒介変数に対する変化量を押し込めています。媒介変数(時刻)の変化量は変えず、それにかかる係数に意味を見出します。例えば重力加速度がそうですかね。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です