龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有： yi+1=yi+hki
其中h为步长，则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同，即局部截断误差为O(h2)。 当用点xi处的斜率近似值k1与右端点xi+1处的斜率k2的算术平均值作为平均斜率k∗的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h(k1+k2)
其中k1=f(xi,yi) ， k2=f(xi+h,yi+hk1) 依次类推，如果在区间[xi,xi+1]内多预估几个店上的斜率值k1,k2,…,km，并用他们的加权平均数作为平均斜率k∗的近似值，显然能够构造出具有很高精度的高阶计数公式。 上述两组公式在形式删过的共同点：都是用f(x,y)在某些点上值得线性组合得出y(xi+1)的近似值yi+1，且增加计算的次数，可以提高截断误差的阶，他们的误差估计可以用f(x,y)在xi处的Taylor展开来估计。
于是可考虑用函数f(x,y)在若干点上的函数值的线性组合老构造金斯公式，构造时要求近似公式在f(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使金斯公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在[xi,xi+1]这一步内计算多个点的斜率值，若够将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法。
一般的龙格-库塔法的形式为
称为P阶龙格-库塔方法。其中ai,bij,cj为待定参数，要求上式yi+1在点(xi,yi)处作Taylor展开，通过相同项的系数确定参数。
Runge--kutta算法C语言实现：
#include "stdio.h"#include "stdlib.h"
void RKT(t,y,n,h,k,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=malloc(n*sizeof(double)); /*分配存储空间*/
if(b == NULL)
{
printf("内存分配失败\n");
exit(1);
}
d=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=0; i<=n-1; i++)
z[i*k]=y[i]; /*将初值赋给数组z的相应位置*/
for(l=1; l<=k-1; l++)
{
Func(y,d);
for (i=0; i<=n-1; i++)
b[i]=y[i];
for (j=0; j<=2; j++)
{
for (i=0; i<=n-1; i++)
{
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=0; i<=n-1; i++)
y[i]=b[i]+h*d[i]/6.0;
for(i=0; i<=n-1; i++)
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,z);
printf("\n");
for (i=0; i<=10; i++) /*打印输出结果*/
{
t=i*h;
printf("t=%5.2f\t ",t);
for (j=0; j<=2; j++)
printf("y(%d)=%e ",j,z[j][i]);
printf("\n");
}
}void Func(y,d)double y[],d[];
{
d[0]=y[1]; /*y0'=y1*/
d[1]=-y[0]; /*y1'=y0*/
d[2]=-y[2]; /*y2'=y2*/
return;
}
In real-time Runge-Kutta algorithm for control system,for satisfying the quickness need of real-time simulation,it is expect to choose larger integration step-size.However,the larger step-size would result in the unsteadiness of numerical value and larger error in numeration. So,based on the existing low-order real-time RK formula,using the stability equation of RK formula,the maximum stability region is found.Then,at the basis of the relation of truncation error and related coefficients,a problem of restricted optimization for min is gotten,and the real-time optimum third-grade second-order RK formula and fourth-grade third-order RK formula are deduced finally.The simulation results show that this algorithm is superior in a certain extent.本回答被网友采纳