程序求助

F_Inzaghi

用C语言编写的程序,结果和matlab差很多画相图的程序
#include <stdio.h>
#include <math.h>           /* definitions for math functions */

/*--- function prototypes: ------------------------------------------*/
double max(double,double);
double min(double,double);
/*-------------------------------------------------------------------*/

#define NDIM 7
#define NPARAM 1




void
function(double *v, double time, double *param, double *deriv){
double x1, x2, x3, x4, x5, x6, x7, R;

  x1 = v[0]; x2 = v[1]; x3 = v[2]; x4 =v[3]; x5 =v[4]; x6 =v[5]; x7=v[6];
  R= param[0];
  
        deriv[0] = -16*x1+R*(-x2*x3+3*x3*x4+8*x6*x7)-4;
        deriv[1] = -10*x2+R*(2.8*x1*x3-1.2*x3*x5+7.2*x4*x6);
        deriv[2] = -2*x3-R*(6*x1*x2+10*x1*x4+2*x2*x5+14*x4*x7);
        deriv[3] = -26*x4+R*(-(14/13)*x1*x3-(4/13)*x2*x6+(38/13)*x3*x7);
        deriv[4] = -8*x5+R*(4*x1*x6+2*x2*x3);
        deriv[5] = -8*x6-R*(4*x1*x5+12*x1*x7+8*x2*x4);
        deriv[6] = -40*x7-R*(0.8*x1*x6+1.2*x3*x4);
}

void
runge_kutta (double *x,double *t, double *tau, double *param, double *xout){
double                 F1[NDIM], F2[NDIM], F3[NDIM], F4[NDIM], xtemp[NDIM];
double                 t_half, half_tau, t_full;
int                    i;

        half_tau = 0.5 * *tau;
        function (x, *t, param, F1);
        t_half = *t + half_tau;
        for (i = 0;i < NDIM; i++) xtemp = x + half_tau*F1;
        function (xtemp, t_half, param, F2);
        for (i = 0;i < NDIM; i++) xtemp = x + half_tau*F2;
        function (xtemp, t_half, param, F3);
        t_full = *t + *tau;
        for (i = 0;i < NDIM; i++) xtemp = x + *tau * F3;
        function (xtemp, t_full, param, F4);
        for (i = 0;i < NDIM; i++)
    xout = x + *tau/6 *(F1 + F4 + 2*(F2 + F3));
}


int
main(int argc, char **argv) {
double   x[NDIM]={-0.25,0.0,-0.1,0.0,0.0,0.0,0.0}, param[NDIM]={44.5},xout[NDIM];
double   time=0.0,tau=0.0001;
int      i,j;
FILE     *fout;


      fout=fopen("phase.dat","w");
      for (i=0;i<500000;i++){
          runge_kutta(x,&time,&tau,param,xout);
          for (j=0;j<NDIM;j++){
               x[j]=xout[j];
               }
          time=time+tau;
          runge_kutta(x,&time,&tau,param,xout);
          i=i+1;
          if (i==10000){
               time=0.0;
               tau=0.001;
               }
          if (i>150000){
              fprintf(fout,"%f %f %f %f\n",time,x[0],x[1],x[2]);
              //printf("%f %f %f %f\n",time,x[0],x[1],x[2]);
              }
          }
      fclose(fout);
}
对了，变步长的也用了，结果也不对,不知道到底是怎么回事.请高人指点
E-mail:F_inzaghi@126.com