经典四阶龙格库塔法

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法,经常被称为“RK4”或者就是“龙格库塔法”。

初值问题表述如下。

对于该问题的RK4由如下方程给出:

其中:

RK4法是四阶方法,也就是说每步的误差是h5,而总积累误差为h4阶。

RK4计算程序如下:

//visualsan@yahoo.cn
#include <iostream>
using namespace std;

#include <vector>
#include <math.h>
using namespace std;
class RK
{
public:
class DiffFunc
{
public:
double operator()(double x,double y)
{
// y'=cos(x)
return cos(x);
// y'=x*y-1
// return x*y-1;
}
} m_df;

RK(double xend,double x0,double y0,double h=1e-3)
{
m_max_x = xend;
m_x0 = x0;
m_y0 = y0;
m_h = h;
m_half_h=h/2.0;
}

void Solve()
{
double yn=m_y0,xstart=m_x0;
while( xstart<m_max_x)
{
double y=yn+K(xstart,yn)*m_h;//y(n+1)=y(n)+h*y'
m_x.push_back(xstart);
m_y.push_back(yn);
cout<<xstart<<""<<yn<<endl;
yn=y;
xstart+=m_h;
}
}
//求解后的点
std::vector<double> m_x,m_y;
//步长h
double m_h;
//初始点x0,y0
double m_x0,m_y0;
//x范围
double m_max_x;
private:
double m_half_h;
double m_ptx,m_pty;
double K1(double x,double y)
{
double v=m_df(x,y);
m_ptx=x;
m_pty=y;
return v;
}
double K2(double _k1)
{
double v=m_df(m_ptx+m_half_h,m_pty+m_half_h*_k1);
return v;
}
double K3(double _k2)
{

double v=m_df(m_ptx+m_half_h,m_pty+m_half_h*_k2);
return v;
}
double K4(double _k3)
{
double v=m_df(m_ptx+m_h,m_pty+m_h*_k3);
return v;
}
double K(double x,double y)
{
double _k1=K1(x,y),_k2=K2(_k1),_k3=K3(_k2),_k4=K4(_k3);
return (_k1+2.0*_k2+2.0*_k3+_k4)/6.0;
}

};

main()
{
RK s1(1,0,0);
s1.Solve();
}

计算实例:

y'=cos(x) y(0)=0, -> solution y=sin(x)

y'=x*cos(x) y(0)=0

y'=1.0/sqrt(x*x+y*y),y(0)=0.1


---------------------------------------------------------------------------------

----------------------------------------------------visualsan@yahoo.cn--------



优质内容筛选与推荐>>
1、JS深度比较两个对象是否相等
2、python代码风格指南:pep8 中文翻译
3、sort() 用法
4、Python简介
5、嵌入式领域TEST


长按二维码向我转账

受苹果公司新规定影响,微信 iOS 版的赞赏功能被关闭,可通过二维码转账支持公众号。

    阅读
    好看
    已推荐到看一看
    你的朋友可以在“发现”-“看一看”看到你认为好看的文章。
    已取消,“好看”想法已同步删除
    已推荐到看一看 和朋友分享想法
    最多200字,当前共 发送

    已发送

    朋友将在看一看看到

    确定
    分享你的想法...
    取消

    分享想法到看一看

    确定
    最多200字,当前共

    发送中

    网络异常,请稍后重试

    微信扫一扫
    关注该公众号