3.3.2 Runge-Kutta法のCプログラム

Euler法のプログラムをたたき台にして Runge-Kutta 法のプログラムを作る。 2次元であるから、それほど複雑ではない (しかしエラーなしでやるのは結構難しい)。

rk2ex1.c

/*
 * rk2ex1.c (Runge-Kutta method for van der Pol equation)
 */

#include <stdio.h>

double mu = 1.0;

int main(void)
{
  int i, N;
  double t, x, y, dt, k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y;
  double fx(double, double, double), fy(double, double, double), x0, y0;
  double Tmax;
  // 初期値
  x0 = 0.1; y0 = 0.1;
  // 最終時刻
  Tmax = 50.0;
  // 時間刻み
  printf("# N: "); scanf("%d", &N);
  dt = Tmax / N;
  // 初期値
  t = 0.0;
  x = x0;
  y = y0;
  printf("# t x y\n");
  printf("%f %f %f\n", t, x, y);
  // Runge-Kutta 法
  for (i = 0; i < N; i++) {
    k1x = dt * fx(x, y, t);
    k1y = dt * fy(x, y, t);
    k2x = dt * fx(x + k1x / 2, y + k1y / 2, t + dt / 2);
    k2y = dt * fy(x + k1x / 2, y + k1y / 2, t + dt / 2);
    k3x = dt * fx(x + k2x / 2, y + k2y / 2, t + dt / 2);
    k3y = dt * fy(x + k2x / 2, y + k2y / 2, t + dt / 2);
    k4x = dt * fx(x + k3x, y + k3y, t + dt);
    k4y = dt * fy(x + k3x, y + k3y, t + dt);
    x += (k1x + 2 * k2x + 2 * k3x + k4x) / 6;
    y += (k1y + 2 * k2y + 2 * k3y + k4y) / 6;
    t = (i + 1) * dt;
    printf("%f %f %f\n", t, x, y);
  }
  return 0;
}

double fx(double x, double y, double t)
{
  return y;
}

double fy(double x, double y, double t)
{
  return - x + mu * (1.0 - x * x) * y;
}

使い方 (コンパイル、実行、可視化) は、 Euler 法のプログラム euler2ex1.c と同じである。

rk2ex1.data という名前のファイルにデータを記録したとして

test4a.gp

splot "rk2ex1.data" with l
set term png
set output "rk2ex1_txy.png"
replot

test4b.gp

plot "rk2ex1.data" with l, "rk2ex1.data" using 1:3 with l
set term png
set output "rk2ex1_tx_ty.png"
replot

test4c.gp

plot "rk2ex1.data" using 2:3 with l
set term png
set output "rk2ex1_xy.png"
replot