微分方程式
ルンゲクッタ(RK4)法; Runge=Kutta Method

概要

---

ルンゲクッタ(RK4)法とは微分方程式の数値解法の一つ。オイラー法やホイン法にくらべ誤差が少ない。
この手法よりも誤差が少ないさらに高次のルンゲクッタ法や、時間刻み幅を調整することができる適応型ルンゲクッタ法があるが、ルンゲクッタ(RK4)法以上の精度が要求されることは少ない。

微分方程式と初期状態が以下に与えられるとき、 $$ \frac{d v}{d t} = f(t, v), \quad v(t_0) = v_0 $$
ルンゲクッタ(RK4)法では\(t\)の刻み幅を\(h=t_{n+1}-t_n\)としたとき、\(v_n\)の次状態\(v_{n+1}\)を以下として算出する。 $$ \begin{eqnarray} \Delta v_1 &=& h f(t_n, v_n) \\ \Delta v_2 &=& h f(t_n + h / 2, v_n + \Delta v_1 / 2) \\ \Delta v_3 &=& h f(t_n + h / 2, v_n + \Delta v_2 / 2) \\ \Delta v_4 &=& h f(t_n + h, v_n + \Delta v_3) \\ v_{n+1} &=& v_n + \frac{\Delta v_1 + 2 \Delta v_2 + 2 \Delta v_3 + \Delta v_4}{6} \end{eqnarray} $$
ソースコード

---

namespace ODEs {

    /// <summary>ルンゲクッタ(RK4)法</summary>
    public class RungeKuttaMethod : ODEsMethod{

        /// <summary>コンストラクタ</summary>
        /// <param name="v0">初期値</param>
        /// <param name="diff_func">微分方程式</param>
        /// <param name="diff_t">時間刻み幅</param>
        public RungeKuttaMethod(double v0, Func<double, double, double> diff_func, double diff_t) : base(v0, diff_func, diff_t) { }

        /// <summary>次状態に更新</summary>
        public override void Reflash() {
            double dv1 = DiffT * diff_func(T, V);
            double dv2 = DiffT * diff_func(T + DiffT / 2, V + dv1 / 2);
            double dv3 = DiffT * diff_func(T + DiffT / 2, V + dv2 / 2);
            double dv4 = DiffT * diff_func(T + DiffT, V + dv3);

            V += (dv1 + 2 * (dv2 + dv3) + dv4) / 6;
            T += DiffT;
        }
    }

    /// <summary>ルンゲクッタ(RK4)法 連立方程式</summary>
    public class SRungeKuttaMethod : SODEsMethod{

        /// <summary>コンストラクタ</summary>
        /// <param name="v0">初期値</param>
        /// <param name="diff_func">微分方程式</param>
        /// <param name="diff_t">時間刻み幅</param>
        public SRungeKuttaMethod(Vector v0, Func<double, Vector, Vector> diff_func, double diff_t) : base(v0, diff_func, diff_t) { }

        /// <summary>次状態に更新</summary>
        public override void Reflash() {
            Vector dv1 = DiffT * diff_func(T, V);
            Vector dv2 = DiffT * diff_func(T + DiffT / 2, V + dv1 / 2);
            Vector dv3 = DiffT * diff_func(T + DiffT / 2, V + dv2 / 2);
            Vector dv4 = DiffT * diff_func(T + DiffT, V + dv3);

            V += (dv1 + 2 * (dv2 + dv3) + dv4) / 6;
            T += DiffT;
        }
    }
}

関連項目

---

微分方程式数値解法 基本クラス
微分方程式数値解法の比較