在for循环中，利用四阶龙阁库塔公式y[i+1]=y[i]+h*z[i]+h*h/6*(L1+L2+L3)z[i+1]=z[i]+h/6*(L1+2L2+3L3+L4)其中L1=f(i,y[i],z[i])L2=f(i+h/2,y[i]+h/2*z[i],z[i]+h/2*L1)L3=f(i+h/2,y[i]+h/2*z[i]+h*h/4*L1,z[i]+h/2*L2)L4=f(i+h,y[i]+h*z[i]+h*h/2*L2,z[i]+h*L3)这样就可以得到需要的值，最后输出下即可。
#include<stdlib.h>#include<stdio.h>
/*n表示几等分，n+1表示他输出的个数*/
int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double))
{
double h=(b-a)/n,k1,k2,k3,k4
int i
// x=(double*)malloc((n+1)*sizeof(double))
// y=(double*)malloc((n+1)*sizeof(double))
x[0]=a
y[0]=y0
switch(style)
{
case 2:
for(i=0i<ni++)
{
x[i+1]=x[i]+h
k1=function(x[i],y[i])
k2=function(x[i]+h/2,y[i]+h*k1/2)
y[i+1]=y[i]+h*k2
}
break
case 3:
for(i=0i<ni++)
{
x[i+1]=x[i]+h
k1=function(x[i],y[i])
k2=function(x[i]+h/2,y[i]+h*k1/2)
k3=function(x[i]+h,y[i]-h*k1+2*h*k2)
y[i+1]=y[i]+h*(k1+4*k2+k3)/6
}
break
case 4:
for(i=0i<ni++)
{
x[i+1]=x[i]+h
k1=function(x[i],y[i])
k2=function(x[i]+h/2,y[i]+h*k1/2)
k3=function(x[i]+h/2,y[i]+h*k2/2)
k4=function(x[i]+h,y[i]+h*k3)
y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6
}
break
default:
return 0
}
return 1
}
double function(double x,double y)
{
return y-2*x/y
}
//例子求y'=y-2*x/y(0<x<1)y0=1
/*
int main()
{
double x[6],y[6]
printf("用二阶龙格-库塔方法\n")
RungeKutta(1,0,1,5,x,y,2,function)
for(int i=0i<6i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i])
printf("用三阶龙格-库塔方法\n")
RungeKutta(1,0,1,5,x,y,3,function)
for(i=0i<6i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i])
printf("用四阶龙格-库塔方法\n")
RungeKutta(1,0,1,5,x,y,4,function)
for(i=0i<6i++)
printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i])
return 1
}
1、首先这是一个纯C的程序，所以你应该使用一个纯C的编译器进行编译，VC实现的是C++，所以会出现上面的错误。我使用gcc编译运行的结果是： （windows下，你可以使用turbo c来编译）t= 0.00 y(0)=-1.000000e+00 y(1)=0.000000e+00 y(2)=1.000000e+00
t= 0.01 y(0)=-9.999500e-01 y(1)=9.999833e-03 y(2)=9.900498e-01
t= 0.02 y(0)=-9.998000e-01 y(1)=1.999867e-02 y(2)=9.801987e-01
t= 0.03 y(0)=-9.995500e-01 y(1)=2.999550e-02 y(2)=9.704455e-01
t= 0.04 y(0)=-9.992001e-01 y(1)=3.998933e-02 y(2)=9.607894e-01
t= 0.05 y(0)=-9.987503e-01 y(1)=4.997917e-02 y(2)=9.512294e-01
t= 0.06 y(0)=-9.982005e-01 y(1)=5.996401e-02 y(2)=9.417645e-01
t= 0.07 y(0)=-9.975510e-01 y(1)=6.994285e-02 y(2)=9.323938e-01
t= 0.08 y(0)=-9.968017e-01 y(1)=7.991469e-02 y(2)=9.231163e-01
t= 0.09 y(0)=-9.959527e-01 y(1)=8.987855e-02 y(2)=9.139312e-01
t= 0.10 y(0)=-9.950042e-01 y(1)=9.983342e-02 y(2)=9.048374e-01
2、如果一定要使用c++，下面是我调试好的程序：（我用g++调通的）：
#include "stdio.h"
#include "stdlib.h"
void Func(double * y, double *d)
// double y[],d[]
{
d[0]=y[1] /*y0'=y1*/
d[1]=-y[0] /*y1'=y0*/
d[2]=-y[2] /*y2'=y2*/
return
}
void RKT(int t, double * y,double n,double h,int k,double * z)
// int n /*微分方程组中方程的个数，也是未知函数的个数*/
// int k /*积分的步数(包括起始点这一步)*/
// double t /*积分的起始点t0*/
// double h /*积分的步长*/
// double y[] /*存放n个未知函数在起始点t处的函数值,返回时,其初值在二维数组z的第零列中*/
// double z[] /*二维数组,体积为n x k.返回k个积分点上的n个未知函数值*/
{
// extern void Func() /*声明要求解的微分方程组*/
int i,j,l
double a[4],*b,*d
b=(double *)malloc(n*sizeof(double)) /*分配存储空间*/
if(b == NULL)
{
printf("内存分配失败\n")
exit(1)
}
d= (double *)malloc(n*sizeof(double)) /*分配存储空间*/
if(d == NULL)
{
printf("内存分配失败\n")
exit(1)
}
/*后面应用RK4公式中用到的系数*/
a[0]=h/2.0
a[1]=h/2.0
a[2]=h
a[3]=h
for(i=0i<=n-1i++)
z[i*k]=y[i] /*将初值赋给数组z的相应位置*/
for(l=1l<=k-1l++)
{
Func(y, d)
for (i=0i<=n-1i++)
b[i]=y[i]
for (j=0j<=2j++)
{
for (i=0i<=n-1i++)
{
y[i]=z[i*k+l-1]+a[j]*d[i]
b[i]=b[i]+a[j+1]*d[i]/3.0
}
Func(y,d)
}
for(i=0i<=n-1i++)
y[i]=b[i]+h*d[i]/6.0
for(i=0i<=n-1i++)
z[i*k+l]=y[i]
t=t+h
}
free(b) /*释放存储空间*/
free(d) /*释放存储空间*/
return
}
main()
{
int i,j
double t,h,y[3],z[3][11]
y[0]=-1.0
y[1]=0.0
y[2]=1.0
t=0.0
h=0.01
RKT(t,y,3,h,11,(double *)z)
printf("\n")
for (i=0i<=10i++) /*打印输出结果*/
{
t=i*h
printf("t=%5.2f\t ",t)
for (j=0j<=2j++)
printf("y(%d)=%e ",j,z[j][i])
printf("\n")
}
}
