龙格库塔_在线百科全书查询
龙格库塔
龙格库塔法的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++编程
