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

1、P7731利用改进欧拉方法计算下列初值问题，并画出近似解的草图:呻y 1,y(0) 3,012tg代码：液进欧拉法fun ctio n Euler(tO,yO,i nv,h)n=rou nd(i nv(2)-i nv(1)/h;t(1)=tO;y(1)=y0;for i=1: ny1(i+1)=y(i)+h*fu n(t(i),y(i);t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fu n(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,得到解

2、析解：hold on;y=dsolve(Dy=y+1,(y(0)=3),t); ezp lot(y,0,2)图像:dy2学y4t,y(0)0.5,0 t 2, t 0.2;dt1代码：fun ctio n Euler1(t0,y0,i nv,h)n=rou nd(i nv(2)-i nv(1)/h;t(1)=t0;y(1)=y0;for i=1: ny1(i+1)=y(i)+h*fu n( t(i),y(i);t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fu n(t(i),y(i)+ fun (t(i+1),y1(i+1)end plot(t,y,*门function y

3、=fun( t,y); y=yA2-4*t;调用:Euler1(0,0,2,图像:-05-1.5-2-250.2040.6 O.B1.21.41.61.S2dy(3) (3 y)(y 1),y(0) 4,0 t 5, t 1; dt代码：fun ctio n Euler2(t0,y0,i nv,h)n=rou nd(i nv(2)-i nv(1)/h;t(1)=t0;y(1)=y0;for i=1: ny1(i+1)=y(i)+h*fu n( t(i),y(i);t(i+1)=t(i)+h;y(i+1)=y(i)+1/2*h*(fu n( t(i),y(i)+ fun( t(i+1),y1(i

4、+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); ezpl ot(y)图像：(A)lo-dz寸H(0)A-(L+A)*(A-e)HACI-)e_OSPHA-coP-OLI-_g-0寸-0)22n山ML+A)*(Ag)HAXAeunEUo4oun4(_*-A)10_dpu(L)LA-(L)l)un4+(M(m)uru)*LI*2/L+(_)AH(L)A二+(_)芥(L)lgA-(m)un芸+(_)AH(L

5、)LALrLU 一046AH(L)A6芥(L)lw(L)u_-(2)u_)punoHU(LI-uroA-9)2_n山UOQOUnAIP-;寸(。)5)(；)刮寸)14.考虑满足初始条件(x(0),y(0)=(1,1)的下列方程组:t=，n=5.用改进欧拉方法求两个方程组的近似解;(1)代码：fun ctio n Euler4(t0,i nt,n ,h)t=t0;x(1)=i nt(1);y(1)=in t(2);for i=1: nx1(i+1)=x(i)+h*xfu n(t(i),x(i),y(i);y1(i+1)=y(i)+h*yfu n(t(i),x(i),y(i);t(i+1)=t(i

6、)+h;x(i+1)=x(i)+1/2*h*(xfu n( t(i),x(i),y(i)+xfu n(t(i+1),x1(i+1),y1(i+1);y(i+1)=y(i)+1/2*h*(yfu n( t(i),x(i),y(i)+yfu n( 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=xfu n( t,x,y);x=y;function y=yfun( t,x,y);y=-2*x-3*y;调用函数：Euler4(0,1,1,5,dxdx(1):dty,2x20.

7、2y xy 1.2y ;选定时间步长图像:(2)代码：fun ctio n Euler5(t0,i nt,n,h) t=t0;x(1)=i nt(1);y(1)=in t(2);for i=1: nx1(i+1)=x(i)+h*xfu n(t(i),x(i),y(i);y1(i+1)=y(i)+h*yfu n( t(i),x(i),y(i);t(i+1)=t(i)+h;x(i+1)=x(i)+1/2*h*(xfu n( t(i),x(i),y(i)+xfu n(t(i+1),x1(i+1),y1(i+1);y(i+1)=y(i)+1/2*h*(yfu n( t(i),x(i),y(i)+yfu n( t(i+1),x1(i+1),y1(i+1); endplot(t,x,o-r)hold onplot(t,y,*-g)h