微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值,这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数,这就是欧拉算法实现的依据。欧拉(Euler)算法是数值求解中毁闭最基本、最简单的方法,但其求解精度较低,一般不在工程中单独进行运算。所谓数值求解,就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。对于常微分方程:
dy/dx=f(x,y),x∈[a,b]
y(a)=y0
可以将区间[a,b]分成n段,那么方程在第xi点有y'(xi)=f(xi,y(xi)),再用向前差商近似代替导数则为:(y(xi+1)-y(xi))/h= f(xi,y(xi)),在这里,h是步长,即相邻两个结点间的距离。因此可以根据xi点和yi点的数值计算出yi+1来:
yi+1= yi+h*f(xi ,yi),i=0,1,2,L
这就是欧拉格式,若初值yi+1是已知的,则可依据上式逐步算出数值解y1,y2,L。
为简化分析,人们常在yi为准确即yi=y(xi)的前提下估计误差y(xi+1)-yi+1,这种误差称为局部截断误差。
如果一种数值方法的局部截断误差为O(h^p+1),则称它的精度是p阶的,或称之为p阶方法。欧拉格式的局部截断误差为O(h^2),由此可知欧拉格式仅为一阶方法。
欧拉公式:
y(xi+1)=yi+h*f(xi,yi)
且xi=x0+i*h (i=0,1,2,…,n-1)
局部截断误差是O(h^2)
改进的欧拉算法
先用纤辩裂欧拉法求得一个初步的近似值,称为预报值,然后用它替代梯形法右端的yi+1再直接计算fi+1,得到校正值yi+1,这样建立的预报-校正系统称为改进的欧拉格式:
预报值 y~i+1=yi+1 + h*f(xi,yi)
校正值 yi+1 =yi+(h/2)*[f(xi,yi)+f(xi+1,y~i+1)]
它有下列平均化形式:
yp=yi+h*f(xi,yi)
且 yc=yi+h*f(xi+1,yp)
且 yi+1=(xp+yc)/2
它的局部截断误差为O(h^3),可见,改进欧拉格式较欧拉格式提高了精度,其截断误差比欧拉格式提高了一阶。
注:欧拉法用差商 [y(xi+1)-y(xi)]/h 近似代替y(xi)的导数,局部截断误差较大;改进欧拉法先用欧拉法求出预报值,再利用梯形公式求出校正值,局部截断误差比欧拉法低了一阶,较大程度地提高了计算精度。
改进欧拉算法
#include<iostream.h>
#define N 20
void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n)
{
int i
float yp,yc,x=x0,y=y0,h=(xn-x0)/n
cout<<"x[0]="<<x<<'t'<<"y[0]"<<灶磨y<<endl
for(i=1i<=ni++)
{
yp=y+h*f1(x,y)
x=x0+i*h
yc=y+h*f1(x,yp)
y=(yp+yc)/2.0
cout<<"x["<<i<<"]="<<x<<"y["<<i<<"]="<<y<<endl
}
}
void main()
{
float xn=5.0,x0=0.0,y0=2.0
float f1(float ,float)
ModEuler(f1,x0,y0,xn,N)
}
float f1(float x,float y)
{
return -x*y*y
}
用Euler折线法求解微分方程,可以按下列步骤求解。
1、确定初值,y0=1
2、确定枯族时间范围,t【0,4】渣败做
3、确定时间步长,h=(tf-t0)/n
4、自定义函数,即 f=func(t,y)
f=@(t,y)y-2*t/y
5、自定义Euler折线法如衡函数,即 [t,y]=euler(func,t0,y0,tf,h)
6、求解主程序
y0=1
t0=0tf=4
n=20
h=(tf-t0)/n
[t,y]=euler(@func,t0,y0,tf,h)
7、执行结果
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)