资源描述:
《Euler法和改进的Euler法实验报告》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、用Euler法和改进的Euler法求u’=-5u(0≤t≤1),u(0)=1的数值解,步长h=0.1,0.05,并比较两个算法的精度。解:1)当步长h=0.1时编写程序如下所示clfclearclc%直接求解微分方程y=dsolve('Dy=-5*y','y(0)=1','t')%Euler法h=0.1;t=0:h:1;n=length(t);u=zeros(1,n);u(1)=1;zbu(1,1)=t(1);zbu(2,1)=u(1);fori=2:nf=-5*u(i-1);u(i)=u(i-1)+h*f;zbu(1,i)=t(i
2、);zbu(2,i)=u(i);endzbu%改进的Euler法v=zeros(1,n);v0=zeros(1,n);v(1)=1;zbv(1,1)=t(1);zbv(2,1)=v(1);fori=2:nf=-5*v(i-1);v0(i)=v(i-1)+h*f;v(i)=v(i-1)+h/2*(f-5*v0(i));zbv(1,i)=t(i);zbv(2,i)=v(i);endzbvplot(t,u,'r*','markersize',10)holdon,plot(t,v,'r.','markersize',20)holdon,ez
3、plot(y,[0,1])holdon,title('Euler法和改进的Euler法比较(h=0.1)),gridonlegend('Euler法','¸改进的Euler法','解析解')%解真值h=0.1;t=0:h:1;n=length(t);fori=1:ny(i)=1/exp(5*t(i));%通过第一部分程序直接解得的解析解zby(1,i)=t(i);zby(2,i)=y(i);endzby我们可以得到计算后的结果图像如图一所示图1Euler法和改进的Euler法比较(h=0.1)同时,我们得到Euler法,改进的Eul
4、er法和解析解的在各点处数值分别如下所示:t坐标0.00.10.20.30.40.50.60.70.80.91.0欧拉1.00000.50000.25000.12500.06250.03130.01560.00780.00390.00200.0010改进欧拉1.00000.62500.39060.24410.15260.09540.05960.03730.02330.01460.0091真值1.00000.60650.36790.22310.13530.08210.04980.03020.01830.01110.0067表1Eule
5、r法和改进的Euler法在各点数值比较(h=0.1)为了比较Euler法和改进的Euler法的算法精度,在这里我们利用相对误差的概念进行评判。对于Euler法和改进的Euler法的每个的估计值有:相对误差=估计值-真值真值从而我们可以通过计算得到如下的相对误差表:t坐标0.00.10.20.30.40.50.60.70.80.91.0欧拉00.17560.32040.43980.53820.61930.68620.74130.78670.82420.8551改进欧拉00.03050.06180.09420.12750.16180.1
6、9720.23360.27120.30990.3498表2Euler法和改进的Euler法在各点相对误差比较(h=0.1)为了评定算法精度,我们对每种算法的在所有点处的相对误差求平均,可以得到Euler法的平均相对误差为0.5443,改进的Euler法的平均相对误差为0.1670。由此我们可以得出改进的欧拉法的算法进度更高。1)当步长h=0.05时程序编写如下clfclearclc%直接求解微分方程y=dsolve('Dy=-5*y','y(0)=1','t')%Euler法h=0.01;t=0:h:1;n=length(t);u=
7、zeros(1,n);u(1)=1;zbu(1,1)=t(1);zbu(2,1)=u(1);fori=2:nf=-5*u(i-1);u(i)=u(i-1)+h*f;zbu(1,i)=t(i);zbu(2,i)=u(i);endzbu%改进的Euler法v=zeros(1,n);v0=zeros(1,n);v(1)=1;zbv(1,1)=t(1);zbv(2,1)=v(1);fori=2:nf=-5*v(i-1);v0(i)=v(i-1)+h*f;v(i)=v(i-1)+h/2*(f-5*v0(i));zbv(1,i)=t(i);zb
8、v(2,i)=v(i);endzbvplot(t,u,'r*','markersize',10)holdon,plot(t,v,'r.','markersize',20)holdon,ezplot(y,[0,1])holdon,tit