龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有：
yi+1=yi+hkiyi+1=yi+hki 其中 hh为步长，则 yi+1yi+1的表达式与 y(xi+1)y(xi+1)的Taylor展开式的前两项完全相同，即 局部截断误差为 O(h2)O(h2)。 当用点 xixi处的斜率近似值 k1k1与右端点 xi+1xi+1处的斜率 k2k2的算术平均值作为平均斜率 k∗k∗的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h(k1+k2)yi+1=yi+h(k1+k2) 其中 k1=f(xi,yi)k1=f(xi,yi) ， k2=f(xi+h,yi+hk1)k2=f(xi+h,yi+hk1) 依次类推，如果在区间 [xi,xi+1][xi,xi+1]内多预估几个点上的斜率值 k1,k2,…,kmk1,k2,…,km，并用他们的加权平均数作为平均斜率 k∗k∗的近似值，显然能够构造出具有很高精度的高阶计数公式。 上述两组公式在形式删过的共同点：都是用 f(x,y)f(x,y)在某些点上值得线性组合得出 y(xi+1)y(xi+1)的近似值 yi+1yi+1，且增加计算的次数，可以提高截断误差的阶，他们的误差估计可以用 f(x,y)f(x,y)在 xixi处的Taylor展开来估计。于是可考虑用函数f(x,y)f(x,y)在若干点上的函数值的线性组合老构造金斯公式，构造时要求近似公式在f(xi,yi)f(xi,yi)处的Taylor展开式与解y(x)y(x)在xixi处的Taylor展开式的前面几项重合，从而使金斯公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在[xi,xi+1][xi,xi+1]这一步内计算多个点的斜率值，若够将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法。 一般的龙格-库塔法的形式为
称为P阶龙格-库塔方法。 其中 ai,bij,cjai,bij,cj为待定参数，要求上式 yi+1yi+1在点 (xi,yi)(xi,yi)处作Taylor展开，通过相同项的系数确定参数。当然，经典的龙格-库塔方法是四阶的。也就是在[xi,xi+1][xi,xi+1]上用四个点处的斜率加权平均作为平均斜率k∗k∗的近似值，构成一系列四阶龙格-库塔公式。具有四阶精度，即局部截断误差是O(h5)O(h5)。 下面介绍最常用的一种四阶龙格-库塔方法。 设
yi+1=yi+c1K1+c2K2+c3K3+c4K4yi+1=yi+c1K1+c2K2+c3K3+c4K4 这里 K1,K2,K3,K4K1,K2,K3,K4为四个不同点上的函数值，分别设其为 其中 c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43均为待定系数。 把 K2,K3,K4K2,K3,K4分别在 xixi点占城h的幂级数，带入线性组合式中，将得到的公式与 y(xi+1)y(xi+1)在 xixi点上的泰勒展开式比较，使其两式右端知道 h4h4的系数相等，经过较复杂的解方程过程便可得到关于 ai,bij,cjai,bij,cj的一组特解。 a2=a3=b21=b32=12a2=a3=b21=b32=12 b31=b41=b42=0b31=b41=b42=0 a4=b43=1a4=b43=1 c1=c4=16c1=c4=16 c2=c3=13c2=c3=13 带入之后得到龙格-库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应正对问题的具体特点选择适合的算法。对于光滑性不太好的解，最好采用低阶算法而将步长hh取小。
龙格-库塔法的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; } 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889ps:如果有时间的话，可能会回过头来加一分解方程的推到吧…
