Lorenz Equation

/*
Lorenz Model

runge-kutta method 3

*/

#include 
#include 

#define h 0.01
#define r 28.0
#define b 2.66666667

double FUNC1();
double FUNC2();
double FUNC3();

main()
{
  int j; /*j; step number*/
  double ny1=0.0,y1=10.0;
  double ny2=0.0,y2=12.0;
  double ny3=0.0,y3=25.0;
  double nt=0.0,t=0.0;
  double k11=0.0,k12=0.0,k13=0.0,
    k21=0.0,k22=0.0,k23=0.0,
    k31=0.0,k32=0.0,k33=0.0,
    k41=0.0,k42=0.0,k43=0.0;
   
  for(j=0; j<=10000;j++)
    {
      k11=h*FUNC1(t,y1,y2,y3);
      k12=h*FUNC2(t,y1,y2,y3);
      k13=h*FUNC3(t,y1,y2,y3);

      k21=h*FUNC1(t+0.5*h,y1+0.5*k11,y2+0.5*k12,y3+0.5*k13);
      k22=h*FUNC2(t+0.5*h,y1+0.5*k11,y2+0.5*k12,y3+0.5*k13);
      k23=h*FUNC3(t+0.5*h,y1+0.5*k11,y2+0.5*k12,y3+0.5*k13);


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

      k41=h*FUNC1(t+h,y1+k31,y2+k32,y3+k33);
      k42=h*FUNC2(t+h,y1+k31,y2+k32,y3+k33);
      k43=h*FUNC3(t+h,y1+k31,y2+k32,y3+k33);

      ny1=y1+0.166666667*(k11+2.0*k21+2.0*k31+k41);
      ny2=y2+0.166666667*(k12+2.0*k22+2.0*k32+k42);
      ny3=y3+0.166666667*(k13+2.0*k23+2.0*k33+k43);

      nt=t+h;
     
      /*printf("%d %f \n",j,y1);*/
      /*printf("%f %f \n",y2,y3);*/
      printf("%f %f %f \n",y1,y2,y3);

      y1=ny1;
      y2=ny2;
      y3=ny3;

      t=nt;
    }
  return 0;
}

/*FUNC 1 */
double FUNC1(double t, double y1, double y2, double y3)
{
  return ( 10.0*(y2-y1) );
}


/*FUNC 2 */
double FUNC2(double t, double y1, double y2, double y3)
{
  return ( -y1*y3+r*y1-y2 );
}

/*FUNC 3 */
double FUNC3(double t, double y1, double y2, double y3)
{
  return ( y1*y2-b*y3 );
}
出力結果
Figure 1: r=28, b=8/3. x0=10, y0=12. z0=25
Figure 2: Y, Z

Back to C Language