3.3.0.1 tanshin-rk.c

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

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

double omega2; // ω^2

int main(void)
{
  /* 開始時刻と終了時刻 */
  double a = 0.0, b = 1.0;
  /* 変数と関数の宣言 */
  int N, j;
  double t,x,y,h,x0,y0,fx(double, double, double), fy(double, double, double);
  /* Runge-Kutta法に必要な変数 */
  double k1x,k1y,k2x,k2y,k3x,k3y,k4x,k4y;
  /* 初期値 */
  x0 = 1.0; y0 = 0.0;
  /* 角振動数ω */
  double omega = 1.0;
  omega2 = omega * omega;
  /* 区間の分割数 N を入力してもらう */
  printf("N="); scanf("%d", &N);
  /* 小区間の幅 */
  h = (b-a) / N;
  /* 開始時刻と初期値のセット */
  t=a;
  x=x0; y=y0;
  printf("t=%f, (x,y)=(%f,%f)\n", t, x, y);
  /* Euler 法による計算 */
  for (j = 0; j < N; j++) {
    k1x = h * fx(t, x, y);
    k1y = h * fy(t, x, y);
    k2x = h * fx(t + h / 2, x + k1x / 2, y + k1y / 2);
    k2y = h * fy(t + h / 2, x + k1x / 2, y + k1y / 2);
    k3x = h * fx(t + h / 2, x + k2x / 2, y + k2y / 2);
    k3y = h * fy(t + h / 2, x + k2x / 2, y + k2y / 2);
    k4x = h * fx(t + h, x + k3x, y + k3y);
    k4y = h * fy(t + h, x + k3x, y + k3y);
    x += (k1x + 2 * k2x + 2 * k3x + k4x) / 6;
    y += (k1y + 2 * k2y + 2 * k3y + k4y) / 6;
    t += h;
    //printf("t=%f, (x,y)=(%f,%f)\n", t, x, y);
    printf("t=%f, (x,y)=(%f,%f) (%f,%f)\n", t, x, y,
           x0 * cos(omega *t) + y0 / omega * sin(omega * t),
           - x0 * omega * sin(omega * t) + y0 * cos(omega * t));
  }
  return 0;
}

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

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

結果を可視化すると面白い。 単振動の解は周期関数であるから、曲線 $(x(t),y(t))$ は閉曲線になるはずであるが、 Euler法は比較的早くズレが生じる。