对于常微分方程的数值解法,四阶Runge—Kutta是一个常用的方法,其精度相对较高,实现起来简单,因而有广泛应用。
下面我们简要讨论一下。
设一阶常微分方程:
u'=f(t,u) a<t<b
u(t(0))=u(0)
Runge-Kutta非线性高阶单步法,p阶R-K法的整体阶段误差为O(h^p)
R-K四阶算法为:
u(i+1)=u(i)+h*(k1+3*k2+3*k3+k4)/8
k1=f(t(i),u(i))
k2=f(t(i+h/3),u(i+h*k1/3))
k3=f(t(i+h/3),u(i+h*k2/3))
k4=f(t(i+h),u(i+h*k3)) */
#include <iostream>
#include <cmath>
using namespace std;
class RK
{
private:
double k1,k2,k3,k4;
double h,b,u,a;
public:
void seth(double l=0){h=l;} //设步长
void setf(double xa=0,double xb=0,double y=0) //设初值和范围(xa,xb)
{
b=xb;
a=xa;
u=y;
}
double f(double t,double u) //函数值,修改它以适应各自需要
{
//函数设定
double f=u-2*t/u;
return f;
}
/**//*---------------------------*/
void dork() //R-K 主函数
{
for(int count=0;count<(b-a)/h;count++)
{
k1=f(a+count*h,u);
k2=f(a+count*h+h/3,u+h*k1/3);
k3=f(a+count*h+2*h/3,u-h*k1/3+h*k2);
k4=f(a+count*h+h,u+h*k1-h*k2+h*k3);
u=u+h*(k1+3*k2+3*k3+k4)/8;
cout<<u<<endl;
}
}
};
void main()
{
RK my;
my.seth(0.1);
my.setf(0,1,1);
my.dork();
}
该程序对数据直接进行显示,如要画图,可以加入几行,输出数据,然后进行画图。