公式

龙格-库塔（Runge-Kutta）法

龙格库塔法的c++编程

龙格库塔法的c++编程

CODE:

#include<stdlib.h>

#include<stdio.h>

/*n表示几等分，n+1表示他输出的个数*/

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[0]=a;

y[0]=y0;

switch(style)

{

case 2:

for(i=0;i<n;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<n;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<n;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<x<1);y0=1;

/*

int main()

{

double x[6],y[6];

printf("用二阶龙格-库塔方法\");

RungeKutta(1,0,1,5,x,y,2,function);

for(int i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\",i,x[i],i,y[i]);

printf("用三阶龙格-库塔方法\");

RungeKutta(1,0,1,5,x,y,3,function);

for(i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\",i,x[i],i,y[i]);

printf("用四阶龙格-库塔方法\");

RungeKutta(1,0,1,5,x,y,4,function);

for(i=0;i<6;i++)

printf("x[%d]=%f,y[%d]=%f\",i,x[i],i,y[i]);

return 1;

龙格库塔求解微分方程数值解

工程中很多的地方用到龙格库塔求解微分方程的数值解，

龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。

此代码只是演示求一个微分方程

/*y'=y-2x/y,x∈[0,0.6]

/*y(0)=1

的解，要求解其它的微分方程，可以自己定义借口函数，替换程序里面的函数：

float f(float , float)即可；

CODE:

/*

/*编程者：

/*完成日期：XXXX，XX，XX

/*e.g:y'=y-2x/y,x∈[0,0.6]

/*　y(0)=1

/*使用经典四阶龙格-库塔算法进行高精度求解

/* y(i+1)=yi+( K1+ 2*K2 +2*K3+ K4)/6

/* K1=h*f(xi,yi)

/* K2=h*f(xi+h/2,yi+K1/2)

/* K3=h*f(xi+h/2,yi+K2/2)

/* K4=h*f(xi+h,yi+K3)

*/

#include <stdio.h>

#include <math.h>

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[3]; //用于存放计算出的常微分方程数值解

int i=0;

int j;

//以下为函数的接口

printf("intput x0:");

scanf("%f",&x0);

printf("input x1:");

scanf("%f",&x1);

printf("input y[0]:");

scanf("%f",&y[0]);

printf("input h:");

scanf("%f",&h);

// 以下为核心程序

n=((x1-x0)/h);

n=3;

for(j=0;j<n;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\",x0,y[i+1]);

++i;

}

}

定义

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。

公式

对于一阶精度的欧拉公式有：

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)

龙格-库塔（Runge-Kutta）法

到目前为止，我们已经学习了多步法，例如：亚当斯-巴什福思（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++编程