龙格库塔法-原理及程序实现

龙格库塔法

一、基本原理：

“龙格-库塔(Runge-Kutta)方法”是一种在工程上应用广泛的“高精度单步算法”。由于此算法精度高，采取措施对误差进行抑制，

可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即：

yi＋１＝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）

计算公式（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

x0y0k1k2k3 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.2500000.780556-0.761584

-0.645858-0.673410

-0.562189

1.5000000.615459-0.568185

-0.481668-0.500993

-0.420537

1.7500000.492374-0.424257

-0.361915-0.374868

-0.317855

2.0000000.400054-0.320087

-0.275466-0.284067

-0.243598

2.2500000.329940-0.244935

-0.212786-0.218538

-0.189482

2.5000000.275895-0.190295

-0.166841-0.170744

-0.149563

2.7500000.233602-0.150068

-0.132704-0.135399

-0.119703

3.0000000.200020-0.120024

-0.106973-0.108868

-0.097048

3.2500000.172989-0.097256

-0.087300-0.088657

-0.079618

3.5000000.150956-0.079757

-0.072054-0.073042

-0.066030

3.7500000.132790-0.066124

-0.060087-0.060818

-0.055305

4.0000000.117655-0.055371

-0.050580-0.051129

-0.046743

4.2500000.104924-0.046789

-0.042945-0.043363

-0.039833

4.5000000.094123-0.039866

-0.036750-0.037072

-0.034202

4.7500000.084885-0.034226

-0.031675-0.031926

-0.029571

xn=5.000000yn=0.076927