如何用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);  
}

我没学过龙格库塔法，我只能帮你该一下，剩下的靠你自己了
#include <stdio.h>
 #include <math.h>
 #define N 100000
 #define h 0.01
double f(double x1,double y1)
 {
    return(10*(y1-x1));
 }
 double p(double x2,double y2,double z2)
 {
    return(28*x2-x2*z2-y2);
 }
double q(double x3,double y3,double z3)
 {
    return(x3*y3-8/3*z3);
 }
 double  runge_kutta(double x0,double y0,double z0)
 {  
    double x1,y1,x2,y2,z2,x3,y3,z3,a,b,c;
    double 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");
        return 0;
    }
 }
void main()
 {  
    double x0=1,y0=1,z0=1;
 
   runge_kutta(x0,y0,z0);  
 }

谢谢 哈