常微分方程作业欧拉法与改进欧拉法.doc

P77 31.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：代码：%改进欧拉法function Euler(t0,y0,inv,h)n=round(inv(2)-inv(1)/h;t(1)=t0;y(1)=y0;for i=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)endplot(t,y,*r)function y=fun(t,y);y=y+1;调用：Euler(0,3,0,2,0.5)得到解析解：hold on;y=dsolve(Dy=y+1,(y(0)=3),t);ezplot(y,0,2)图像：代码：function Euler1(t0,y0,inv,h)n=round(inv(2)-inv(1)/h;t(1)=t0;y(1)=y0;for i=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)endplot(t,y,*r)function y=fun(t,y);y=y2-4*t;调用：Euler1(0,0.5,0,2,0.2)图像：代码：function Euler2(t0,y0,inv,h)n=round(inv(2)-inv(1)/h;t(1)=t0;y(1)=y0;for i=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)endplot(t,y,*r)function y=fun(t,y);y=(3-y)*(y+1);调用：Euler2(0,4,0,5,1)得到解析解：hold on;y=dsolve(Dy=(3-y)*(y+1),y(0)=4,t);ezplot(y)图像：代码：function Euler2(t0,y0,inv,h)n=round(inv(2)-inv(1)/h;t(1)=t0;y(1)=y0;for i=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)endplot(t,y,*r)function y=fun(t,y);y=(3-y)*(y+1);调用：Euler2(0,4,0,5,0.5)得到解析解：hold on;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） 代码：function Euler4(t0,int,n,h)t=t0;x(1)=int(1);y(1)=int(2);for i=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);endplot(t,x,o-r)hold onplot(t,y,*-g)hold onplot(x,y)function x=xfun(t,x,y);x=y;function y=yfun(t,x,y);y=-2*x-3*y;调用函数：Euler4(0,1,1,5,0.25)图像：（2）代码：function Euler5(t0,int,n,h)t=t0;x(1)=int(1);y(1)=int(2);for i=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);endplot(t,x,o-r)hold onplot(t,y,*-g)hold onplot(