龍格庫塔法
工程上應用的高精度單步演演算法
龍格-庫塔(Runge-Kutta)方法是一種在工程上應用廣泛的高精度單步演演算法。由於此演演算法精度高，採取措施對誤差進行抑制，所以其實現原理也較複雜。該演演算法是構建在數學支持的基礎之上的。它是如此常用，以至於經常被稱為“RK4”或者就是“龍格庫塔法”。該方法主要是在已知方程導數和初值信息，利用計算機模擬時應用，省去求解微分方程的複雜過程。
龍格庫塔法的c++編程
CODE:
#include
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=a;
y=y0;
switch(style)
{
case 2:
for(i=0;i
{
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=0;i
{
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=0;i
{
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<1);y0=1;
#include
#include
float f(float den,float p0)//要求解的微分方程的右部的函數 e.g: y-2x/y
{
float rus;
// den=w/(W0+sl);
// rus=k*A*f/Wp*sqrt(RTp)-(k-1)*.....
rus=p0-2*den/p0;
return(rus);
}
//float fden()
//{
//}
void main()
{
float x0; //範圍上限
int x1; //範圍下限
float h;//步長
int n; //計算出的點的個數
float k1,k2,k3,k4;
float y; //用於存放計算出的常微分方程數值解
int i=0;
int j;
//以下為函數的介面
printf("intput x0:");
scanf("%f",&x0);
printf("input x1:");
scanf("%f",&x1);
printf("input y:");
scanf("%f",&y);
printf("input h:");
scanf("%f",&h);
// 以下為核心程序
n=((x1-x0)/h);
n=3;
for(j=0;j
{
k1=h*f(x0,y[i]); //求K1
k2=h*f((x0+h/2)，(y[i]+k1/2)); //求K2
k3=h*f((x0+h/2)，(y[i]+k2/2)); //求K3
k4=h*f((x0+h),(y[i]+k3)); //求K4
y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求yi+1
x0+=float(0.2);
printf("y[%f]=%f\n",x0,y[i+1]);
++i;
對於一階精度的歐拉公式有:
yi+1=yi+h*K1
K1=f(xi,yi)
當用點xi處的斜率近似值K1與右端點xi+1處的斜率K2的算術平均值作為平均斜率K*的近似值，那麼就會得到二階精度的改進歐拉公式:
yi+1=yi+h*(K1+ K2)/2
K1=f(xi,yi)
K2=f(xi+h,yi+h*K1)
依次類推，如果在區間[xi,xi+1]內多預估幾個點上的斜率值K1、K2、……Km，並用他們的加權平均數作為平均斜率K*的近似值，顯然能構造出具有很高精度的高階計算公式。經數學推導、求解，可以得出四階龍格-庫塔公式，也就是在工程中應用廣泛的經典龍格-庫塔演演算法:
yi+1=yi+h*(K1+ 2*K2 +2*K3+ K4)/6
K1=f(xi,yi)
K2=f(xi+h/2,yi+h*K1/2)
K3=f(xi+h/2,yi+h*K2/2)
K4=f(xi+h,yi+h*K3)
（1）的局部截斷誤差是。
龍格-庫塔法具有精度高，收斂，穩定（在一定條件下），計算過程中可以改變步長，不需要計算高階導數等優點，但仍需計算 在一些點上的值，如四階龍格-庫塔法每計算一步需要計算四次 的值，這給實際計算帶來一定的複雜性，因此，多用來計算“表頭”。
二、小程序
#include
#include
#define f(x,y) (-1*(x)*(y)*(y))
void main(void)
{
double a,b,x0,y0,k1,k2,k3,k4,h;
int n,i;
printf("input a,b,x0,y0,n:");
scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n);
printf("x0\ty0\tk1\tk2\tk3\tk4\n");
for(h=(b-a)/n,i=0;i!=n;i++)
{
k1=f(x0,y0);
k2=f(x0+h/2,y0+k1*h/2);
k3=f(x0+h/2,y0+k2*h/2);
k4=f(x0+h,y0+h*k3);
printf("%lf\t%lf\t",x0,y0);
printf("%lf\t%lf\t",k1,k2);
printf("%lf\t%lf\n",k3,k4);
y0+=h*(k1+2*k2+2*k3+k4)/6;
x0+=h;
}
printf("xn=%lf\tyn=%lf\n",x0,y0);
}
運行結果：
input a,b,x0,y0,n:0 5 0 2 20
x0 y0 k1 k2 k3 k4
0.000000 2.000000 -0.000000 -0.500000 -0.469238
-0.886131
0.250000 1.882308 -0.885771 -1.176945 -1.129082
-1.280060
0.500000 1.599896 -1.279834 -1.295851 -1.292250
-1.222728
0.750000 1.279948 -1.228700 -1.110102 -1.139515
-0.990162
1.000000 1.000027 -1.000054 -0.861368 -0.895837
-0.752852
1.250000 0.780556 -0.761584 -0.645858 -0.673410
-0.562189
1.500000 0.615459 -0.568185 -0.481668 -0.500993
-0.420537
1.750000 0.492374 -0.424257 -0.361915 -0.374868
-0.317855
2.000000 0.400054 -0.320087 -0.275466 -0.284067
-0.243598
2.250000 0.329940 -0.244935 -0.212786 -0.218538
-0.189482
2.500000 0.275895 -0.190295 -0.166841 -0.170744
-0.149563
2.750000 0.233602 -0.150068 -0.132704 -0.135399
-0.119703
3.000000 0.200020 -0.120024 -0.106973 -0.108868
-0.097048
3.250000 0.172989 -0.097256 -0.087300 -0.088657
-0.079618
3.500000 0.150956 -0.079757 -0.072054 -0.073042
-0.066030
3.750000 0.132790 -0.066124 -0.060087 -0.060818
-0.055305
4.000000 0.117655 -0.055371 -0.050580 -0.051129
-0.046743
4.250000 0.104924 -0.046789 -0.042945 -0.043363
-0.039833
4.500000 0.094123 -0.039866 -0.036750 -0.037072
-0.034202
4.750000 0.084885 -0.034226 -0.031675 -0.031926
-0.029571
xn=5.000000 yn=0.076927
到目前為止，我們已經學習了多步法，例如：亞當斯-巴什福思（Adams
-Bashorth）法，亞當斯-莫爾頓（Adams-Monlton）法，都是 常微分方程的積分
方法。它們需要在每一次 迭代時重新計算一遍等式右邊的結果（非線性隱含問
題忽略計算多個 f （ω）值的可能性）
龍格-庫塔（Runge-Kutta）法是一種不同的處理，作為多級方法為人們所知。
它要求對於一個簡單的校正計算多個 f 的值。
下面，我們列出了 3 種最流行的龍格-庫塔（Runge-Kutta）法：
改進的 歐拉方法（精度：p=2）：
V a = V n + Δtf (V n,tn)
2
Δt）二階格式
V n+1 = V n +Δtf (V a,tn +
2
Hevn’s 方法（p=2）：
這是另一種二階格式：
V a = V n +Δtf (V n,tn)
V n = V n +
+1 Δt[ f (V n,tn) + f (V a,tn +Δt)]
2
注意：f (Vn,tn）在運算中應該只被計算一次。
四次龍格-庫塔（Runge-Kutta）法（p=4）：
這是一個 4 階格式。這次我們寫的形式有點不同：
a = Δtf (V n,tn)
b = Δtf (V n + 1 a,tn + 1
2 2 Δt)
c = Δtf (V n + 1 b,tn + Δt)
1
2 2
d = Δtf (V n + c,tn +Δt)
V n =V n +
+1 1 (a +2b +2c +d）。
6
龍格庫塔法的 c++編程
基本信息
- 外文名
- Runge-Kutta methods
- 特點
- 精度高
- 提出時間
- 1900年左右
- 應用學科
- 數學
- 類別
- 高精度單步演算法