4.4 C言語によるRunge-Kutta法のプログラム例

1階の微分方程式 (2次元) に帰着されたので、 前節の解説に従えば良いわけだが、 1つくらいプログラム例をあげておく。

rk2ex2.c

/*
 * rk2ex2.c (Runge-Kutta method for free fall)
 */

#include <stdio.h>

double g = 9.8;

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;
  // 初期値 (100m の高さから初速 0 で降りる)
  x0 = 100.0; y0 = 0.0;
  // 分割数, 最終時刻 -> 時間刻み
  printf("# N, Tmax: "); scanf("%d%lf", &N, &Tmax);
  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 - g;
}

コンパイル&実行

% cc -O3 -o rk2ex2 rk2ex2.c
% ./rk2ex2 > rk2ex2.data
500 4.55
%

図 15: 高さ 100m からの自由落下 (横軸時刻, 縦軸高さ)