常微分方程作业欧拉法与改进欧拉法.docx
《常微分方程作业欧拉法与改进欧拉法.docx》由会员分享，可在线阅读，更多相关《常微分方程作业欧拉法与改进欧拉法.docx（9页珍藏版）》请在冰点文库上搜索。
常微分方程作业欧拉法与改进欧拉法
常微分方程作业欧拉法与改进欧拉法
P7731.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：
代码：
%改进欧拉法
functionEuler（t0,y0,inv,h）
n=round（inv
（2）-inv
（1））/h;
t
（1）=t0;
y
（1）=y0;
fori=1:
n
y1（i+1）=y（i）+h*fun（t（i）,y（i））;
t（i+1）=t（i）+h;
y（i+1）=y（i）+1/2*h*（fun（t（i）,y（i））+fun（t（i+1）,y1（i+1）））
end
plot（t,y,'*r'）
functiony=fun（t,y）;
y=y+1;
调用：
Euler（0,3,[0,2],0.5）
得到解析解：
holdon;
y=dsolve（'Dy=y+1','（y（0）=3）','t'）;
ezplot（y,[0,2]）
图像：
代码：
functionEuler1（t0,y0,inv,h）
n=round（inv
（2）-inv
（1））/h;
t
（1）=t0;
y
（1）=y0;
fori=1:
n
y1（i+1）=y（i）+h*fun（t（i）,y（i））;
t（i+1）=t（i）+h;
y（i+1）=y（i）+1/2*h*（fun（t（i）,y（i））+fun（t（i+1）,y1（i+1）））
end
plot（t,y,'*r'）
functiony=fun（t,y）;
y=y^2-4*t;
调用：
Euler1（0,0.5,[0,2],0.2）
图像：
代码：
functionEuler2（t0,y0,inv,h）
n=round（inv
（2）-inv
（1））/h;
t
（1）=t0;
y
（1）=y0;
fori=1:
n
y1（i+1）=y（i）+h*fun（t（i）,y（i））;
t（i+1）=t（i）+h;
y（i+1）=y（i）+1/2*h*（fun（t（i）,y（i））+fun（t（i+1）,y1（i+1）））
end
plot（t,y,'*r'）
functiony=fun（t,y）;
y=（3-y）*（y+1）;
调用：
Euler2（0,4,[0,5],1）
得到解析解：
holdon;
y=dsolve（'Dy=（3-y）*（y+1）','y（0）=4','t'）;
ezplot（y）
图像：
代码：
functionEuler2（t0,y0,inv,h）
n=round（inv
（2）-inv
（1））/h;
t
（1）=t0;
y
（1）=y0;
fori=1:
n
y1（i+1）=y（i）+h*fun（t（i）,y（i））;
t（i+1）=t（i）+h;
y（i+1）=y（i）+1/2*h*（fun（t（i）,y（i））+fun（t（i+1）,y1（i+1）））
end
plot（t,y,'*r'）
functiony=fun（t,y）;
y=（3-y）*（y+1）;
调用：
Euler2（0,4,[0,5],0.5）
得到解析解：
holdon;
y=dsolve（'Dy=（3-y）*（y+1）','y（0）=4','t'）;
ezplot（y）
图像：
14.考虑满足初始条件（x（0），y（0））=（1,1）的下列方程组：
选定时间步长
t=0.25，n=5.用改进欧拉方法求两个方程组的近似解；
（1）代码：
functionEuler4（t0,int,n,h）
t=t0;
x
（1）=int
（1）;
y
（1）=int
（2）;
fori=1:
n
x1（i+1）=x（i）+h*xfun（t（i）,x（i）,y（i））;
y1（i+1）=y（i）+h*yfun（t（i）,x（i）,y（i））;
t（i+1）=t（i）+h;
x（i+1）=x（i）+1/2*h*（xfun（t（i）,x（i）,y（i））+xfun（t（i+1）,x1（i+1）,y1（i+1）））;
y（i+1）=y（i）+1/2*h*（yfun（t（i）,x（i）,y（i））+yfun（t（i+1）,x1（i+1）,y1（i+1）））;
end
plot（t,x,'o-r'）
holdon
plot（t,y,'*-g'）
holdon
plot（x,y）
functionx=xfun（t,x,y）;
x=y;
functiony=yfun（t,x,y）;
y=-2*x-3*y;
调用函数：
Euler4（0,[1,1],5,0.25）
图像：
（2）代码：
functionEuler5（t0,int,n,h）
t=t0;
x
（1）=int
（1）;
y
（1）=int
（2）;
fori=1:
n
x1（i+1）=x（i）+h*xfun（t（i）,x（i）,y（i））;
y1（i+1）=y（i）+h*yfun（t（i）,x（i）,y（i））;
t（i+1）=t（i）+h;
x（i+1）=x（i）+1/2*h*（xfun（t（i）,x（i）,y（i））+xfun（t（i+1）,x1（i+1）,y1（i+1）））;
y（i+1）=y（i）+1/2*h*（yfun（t（i）,x（i）,y（i））+yfun（t（i+1）,x1（i+1）,y1（i+1）））;
end
plot（t,x,'o-r'）
holdon
plot（t,y,'*-g'）
holdon
plot（x,y）
functionx=xfun（t,x,y）;
x=y+y^2;
functiony=yfun（t,x,y）;
y=-x+0.2*y-x*y+1.2*y^2;
调用函数：
Euler5（0,[1,1],5,0.25）
图像：
升级会员 