欧拉法matlab程序.doc

1.Euler法function x,y=naeuler(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1).0.5;plot(x,y,x1,y1) dyfun=inline(y-2*x/y);x,y=naeuler(dyfun,0,1,1,0.2);x,yans = 0 1.0000 0.2000 1.2000 0.4000 1.3733 0.6000 1.5315 0.8000 1.68111.0000 1.82692.隐式Euler法function x,y=naeulerb(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h);endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1).0.5;plot(x,y,x1,y1)function y=iter(dyfun,x,y,h)y0=y;e=1e-4;K=1e+4;y=y+h*feval(dyfun,x,y);y1=y+2*e;k=1; while abs(y-y1)e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if kK error(迭代发散); endend dyfun=inline(y-2*x/y);x,y=naeulerb(dyfun,0,1,1,0.2);x,yans = 0 1.0000 0.2000 1.1641 0.4000 1.3014 0.6000 1.4146 0.8000 1.5019 1.0000 1.5561 3.改进Euler法 function x,y=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2; end x=x;y=y; x1=0:0.2:1;y1=(1+2*x1).0.5; plot(x,y,x1,y1) dyfun=inline(y-2*x/y);x,y=naeuler2(dyfun,0,1,1,0.2);x,yans = 0 1.0000 0.2000 1.1867 0.4000 1.3483 0.6000 1.4937 0.8000 1.6279 1.0000 1.7542