如何用C语言编程龙格库塔法求解常微分方程组
#include <stdio.h>
#include <math.h>
#define N 100000
#define h 0.01
float f(float x1,float y1)
{
   return(10*(y1-x1));
}
float p(float x2,float y2,float z2)
{
   return(28*x2-x2*z2-y2);
}
float q(float x3,float y3,float z3)
{
   return(x3*y3-8/3*z3);
}
float  runge_kutta(float x0,float y0,float z0)
{  
   float x1,y1,x2,y2,z2,x3,y3,z3,a,b,c;
   float d1,d2,d3,d4,k1,k2,k3,k4,f1,f2,f3,f4;
   int i;
      y1=y2=y3=y0;
      x2=x1=x3=x0;
      z2=z3=z0;
   FILE *fp;
   fp=fopen("luo_lunzir.txt","w");
   for(i=1;i<=N;i++)
   {   
      d1=f(x1,y1);
      k1=p(x2,y2,z2);
      f1=q(x3,y3,z3);
      
      d2=f(y1+h/2,x1+h*d1/2);
      k2=p(x2+h/2,z2+h/2,y2+h*k1/2);
      f2=q(x3+h/2,y3+h/2,z3+h*f1/2);
      
      d3=f(y1+h/2,x1+h*d2/2);
      k3=p(x2+h/2,z2+h/2,y2+h*k2/2);
      f3=q(x3+h/2,y3+h/2,z3+h*f2/2);
      
      d4=f(y1+h,x1+h*d3);
      k4=p(x2+h,z2+h,y2+h*k3);
      f4=q(x3+h,y3+h,z3+h*f3);
      
      a=x1+h*(d1+2*d2+2*d3+d4)/6;
      b=y2+h*(k1+2*k2+2*k3+k4)/6;
      c=z3+h*(f1+2*f2+2*f3+f4)/6;
      
       y1=y3=y0+i*h;
       z2=z0+i*h;
      x1=x3=x2=x0+i*h;
      
       fprintf(fp," %8f  , %8f,  %8f\n",a,b,c);
       printf("  %8f ,  %8f ,  %8f",a,b,c);
       printf("\n");
   }
   }
main()
{  
   float x0=1,y0=1,z0=1;

   runge_kutta(x0,y0,z0);  
}

谢谢 哈