## 龍格-庫塔(Runge-Kutta)方法
龍格-庫塔(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展開，通過相同項的系數確定參數。
當然，經典的龍格-庫塔方法是四階的。也就是在[xi,xi+1]上用四個點處的斜率加權平均作為平均斜率k?的近似值，構成一系列四階龍格-庫塔公式。具有四階精度，即局部截斷誤差是O(h5)。
下面介紹最常用的一種四階龍格-庫塔方法。
設
yi+1=yi+c1K1+c2K2+c3K3+c4K4
這里K1,K2,K3,K4為四個不同點上的函數值，分別設其為
其中c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43均為待定系數。
把K2,K3,K4分別在xi點占城h的冪級數，帶入線性組合式中，將得到的公式與y(xi+1)在xi點上的泰勒展開式比較，使其兩式右端知道h4的系數相等，經過較復雜的解方程過程便可得到關于ai,bij,cj的一組特解。
a2=a3=b21=b32=12
b31=b41=b42=0
a4=b43=1
c1=c4=16
c2=c3=13
帶入之后得到
龍格-庫塔方法的推導基于Taylor展開方法，因而它要求所求的解具有較好的光滑性。如果解的光滑性差，那么，使用四階龍格-庫塔方法求得的數值解，其精度可能反而不如改進的歐拉方法。在實際計算時，應正對問題的具體特點選擇適合的算法。對于光滑性不太好的解，最好采用**低階算法**而將步長h取小。
龍格-庫塔法的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;
}
~~~
ps:如果有時間的話，可能會回過頭來加一分解方程的推到吧…
