Numerical analysis of van der Pol equation by using Runge-Kutta method

常微分方程式、坂井秀隆、東京大学出版会に収録されているコードを転載した。

Gnuplotを用いて、出力された２次元のデータを可視化した図も収録されている。

#include <stdio.h>
#include <math.h>
#define xinitval 2
#define yinitval -3
#define step 0.01
#define max 1024


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

double fy(double x, double y){
  return -x + 1.5*(1-x*x)*y;
}

int main(void){
  int j;
  double K1, K2, K3, K4, L1, L2, L3, L4;
  double x[max], y[max];
  x[0]=xinitval, y[0]=yinitval;
  for(j=0; j<max; j++){
    K1=fx(x[j], y[j])*step;
    L1=fy(x[j], y[j])*step;
    K2=fx(x[j]+(K1/2), y[j]+(L1/2))*step;
    L2=fy(x[j]+(K1/2), y[j]+(L1/2))*step;
    K3=fx(x[j]+(K2/2), y[j]+(L2/2))*step;
    L3=fy(x[j]+(K2/2), y[j]+(L2/2))*step;
    K4=fx(x[j]+(K3/2), y[j]+(L3/2))*step;
    L4=fy(x[j]+(K4/2), y[j]+(L4/2))*step;
    x[j+1]=x[j]+(K1+2*K2+2*K3+K4)/6;
    y[j+1]=y[j]+(L1+2*L2+2*L3+L4)/6;
  }
  for(j=1; j<=max; j++){
    printf("%lf $lf\n", x[j], y[j]);
  }
  return 0;
}