2.3.1 Runge-Kutta法のプログラム

test-rk.c

* test-rk.c --- 常微分方程式の初期値問題を Runge-Ketta 法で解く */

#include <stdio.h>
#include <math.h>

int main(void)
{
  /* 開始時刻と終了時刻 */
  double a = 0.0, b = 1.0;
  /* 初期値 */
  double x0;
  /* 変数と関数の宣言 */
  int N,j;
  double t,x,h,f(double,double),k1,k2,k3,k4;
  /* 初期値の設定 */
  x0 = 1.0;
  /* 区間の分割数 N を入力してもらう */
  printf("N="); scanf("%d", &N);
  /* 小区間の幅 */
  h = (b-a) / N;
  /* 開始時刻と初期値のセット */
  t = a;
  x = x0;
  printf("%f %f\n", t, x);
  /* Runge-Kutta 法による計算 */
  for (j = 0; j < N; j++) {
    k1 = h * f(t, x);
    k2 = h * f(t + h / 2, x + k1 / 2);
    k3 = h * f(t + h / 2, x + k2 / 2);
    k4 = h * f(t + h, x + k3);
    x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    t += h;
    printf("%f %f\n", t, x);
  }
  printf("%f %17.15f\n", t, x);
  //printf("       e=%17.15f\n", exp(1.0));
  return 0;
}

/* 微分方程式 x'=f(t,x) の右辺の関数 f の定義 */
double f(double t, double x)
{
  return x;
}

コンパイル＆実行

cc test-rk.c
./a.out

たった $10$ 等分でも、相対誤差が $10^{-6}$ 以下になる。

Runge-Kutta 法の公式は Euler 法よりは大部面倒だが、 それに見合うだけの価値がある。