2.3.2 Runge-Kutta法のプログラム例

Euler法のプログラムを Runge-Kutta法を用いるように書き換えるのは簡単である。

rk1ex1.c

/*
 * rk1ex1.c (Runge-Kutta method for Malthusian model)
 */

#include <stdio.h>

double k = 1.0;

int main(void)
{
  int i, N;
  double t, x, dt, k1, k2, k3, k4;
  double f(double, double), x0;
  double Tmax;
  // 初期値
  x0 = 1.0;
  // 最終時刻
  Tmax = 1.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  printf("# t x\n");
  printf("%f %f\n", t, x);
  // Runge-Kutta 法
  for (i = 0; i < N; i++) {
    k1 = dt * f(x, t);
    k2 = dt * f(x + k1/2, t + dt / 2);
    k3 = dt * f(x + k2/2, t + dt / 2);
    k4 = dt * f(x + k3, t + dt);
    x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    t = (i + 1) * dt;
    printf("%f %20.14f\n", t, x);
  }
  return 0;
}

double f(double x, double t)
{
  return k * x;
}

cc でコンパイル&実行

% cc rk1ex1.c
% ./a.out
# N: 100
# t x
0.000000 1.000000
中略
1.000000 2.718282
%

解曲線の図を描くための手順は Euler 法のときと同様 (データファイルの名前を変えるくらい) なので省略する。

問5      2.718282 は小数点以下6位までしか表示していないが、 小数点以下15位まで表示するように C プログラムを書き直して実行せよ (printf() の書式の選び方の問題, 「浮動小数点数の入出力と四則演算」)。 そうすると 2.718281828234403 となる。 この2.718281828234403 の誤差はいくらか。 $N$ を変えることで誤差がどのように変わるか調べよ (両側対数目盛でプロットすることを勧める, 「対数グラフを描く」)。