对于最简单的常微分方程形式,y'=F(x,y)与y(x0)=y0,我们假定F(x,y)适当光滑,在一定的区间内进行求值。这里讲的是差分的方法,通过寻求一系列离散点求解近似解y1,y2,y3……下面的所有,以函数Y'=Y-2*X/Y,Y(0)=1为例。
1、Euler方法
最原始的Euler格式是用向前差商的方法推导的,Yn+1=Yn+h*F(Xn,Yn)。基于这个道理,我们如果将前面算过的两个值都利用起来,那么我们就可以得到两步Euler格式,Yn+1=Yn-1+2*h*F(Xn,Yn)。同样的,如果用向后差商的方法,我们就可以得到隐式Euler格式,Yn+1=Yn+h*F(Xn,Yn+1),这里要用到Yn+1的值。如果我们将两者进行综合,就可以得到Euler格式的预报校正系统,改进Euler格式,就是将前者的结果代入后者中进行校正检验。此外,我们还可以对这小段区间[Xn,Yn]进行二分法,Yn+1=Yn+h*F(Xn+1/2,Yn+h*F(Xn,Yn)/2),从而得到变形Euler格式。
既然有了方法,那就编程吧。
double Euler1(double x,double y,double h)//改进的欧拉格式
{
double y1=Euler3(x,y,h);
return y+(f(x,y)+f(x+h,y1))*h/2;
}
double Euler2(double x,double y,double h)//变形的欧拉格式
{
double y1=y+f(x,y)*h/2;
return y+h*f(x+h/2,y1);
}
double Euler3(double x,double y,double h)//单步欧拉格式
{
return y+h*f(x,y);
}
double Euler4(double x,double y1,double y2,double h)//两步欧拉格式
{
return y1+2*h*