资源描述
I.Euler 法function x,y=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);y(i)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1)A0.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.00000.2000 1.20000.4000 1.37330.6000 1.53150.8000 1.68111.0000 1.8269 2.隐式Euler法function x,y=naeulerb(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1y(n+1)=iter(dyfun,x(n+1),y(n),h);endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1).A0.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)ey1=y;y=y0+h*feval(dyfun,x,y);k=k+1;if kKerror(迭代发散);endend dyfun=inline(y-2*x/y);x,y=naeulerb(dyfun,0,1,1,0.2);x,yans =0 1.00000.2000 1.16410.4000 1.30140.6000 1.41460.8000 1.50191.0000 1.55613. 改进Euler法function x,y=naeuler2(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=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;endx=x;y=y;x1=0:0.2:1;y1=(1+2*x1).A0.5;plot(x,y,x1,y1) dyfun=inline(y-2*x/y);x,y=naeuler2(dyfun,0,1,1,0.2);x,y ans =0 1.00000.2000 1.18670.4000 1.34830.6000 1.49370.8000 1.62791.0000 1.7542
展开阅读全文