Fitzhugh-Nagumo Equation

/*
runge-kutta method

*/

#include 
#include 

#define h 0.01 /*h*/

#define a 0.25
#define gamma 3.0
#define epsilon 0.001


double FUNC1();
double FUNC2();

main()
{
  int j; /*j; step number*/
  double ny1=0.0,y1=0.6;
  double ny2=0.0,y2=0.005;
  double nt=0.0,t=0.0;
  double k11=0.0,k12=0.0,k21=0.0,k22=0.0,
    k31=0.0,k32=0.0,k41=0.0,k42=0.0;
   
  for(j=0; j<=100000;j++)
    {
      k11=h*FUNC1(t,y1,y2);
      k12=h*FUNC2(t,y1,y2);
      k21=h*FUNC1(t+0.5*h,y1+0.5*k11,y2+0.5*k12);
      k22=h*FUNC2(t+0.5*h,y1+0.5*k11,y2+0.5*k12);

      k31=h*FUNC1(t+0.5*h,y1+0.5*k21,y2+0.5*k22);
      k32=h*FUNC2(t+0.5*h,y1+0.5*k21,y2+0.5*k22);

      k41=h*FUNC1(t+h,y1+k31,y2+k32);
      k42=h*FUNC2(t+h,y1+k31,y2+k32);
     
      ny1=y1+0.166666667*(k11+2.0*k21+2.0*k31+k41);
      ny2=y2+0.166666667*(k12+2.0*k22+2.0*k32+k42);

      nt=t+h;
     
      printf("%d %f \n",j,ny1);
      /*printf("%f %f \n",ny1,ny2);*/

      y1=ny1;
      y2=ny2;
      t=nt;
    }
  return 0;
}

/*function 1 */
double FUNC1(double t, double y1, double y2)
{
  return ( y1*(1.0-y1)*(y1-a)-y2 );
}

/*function 2 */
double FUNC2(double t, double y1, double y2)
{
  return ( epsilon*(y1-gamma*y2) );
}
出力結果
相図
進行パルス波

Back to C Language