计算方法-上机作业.docVIP

用下列方法求dy/dx=2/3*x*y^(-2),0=x=1 ; y(0)=1的数值解（取h=0.1），并将计算结果与准确解y=(1+x^2)^(1/3)进行比较 欧拉方法 改进欧拉方法 经典R-K法 解： 欧拉方法： function [x,y]=Euler1(fun,a,b,y0,n) h=(b-a)/n; x=a:h:b; y=0*x; y(1)=y0; for k=1:n y(k+1)=y(k)+h*feval(fun,x(k),y(k)); end 输入： [x,y]=Euler1(f,0,1,1,10) plot(x,y,or) hold on xx=0:0.01:1; yy=(1+xx.^2).^(1/3); plot(xx,yy) 输出： 改进欧拉算法： function [x,y]=imEuler2(fun,a,b,y0,n) h=(b-a)/n; x=a:h:b; y=0*x; y(1)=y0; for k=1:n s1=feval(fun,x(k),y(k)); ys=y(k)+h*s1; s2=feval(fun,x(k+1),ys); y(k+1)=y(k)+0.5*h*(s1+s2); end 输入： [x,y]=imEuler2(f,0,1,1,10) plot(x,y,or) hold on xx=0: