計算方法上機作業(yè)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上計算方法第四次上機報告2.用歐拉方法解初值y=10x(1-y) 0<=x<=1Y(0)=0取步長h=0.1,保留5位有效數(shù)字，并與準(zhǔn)確解相比較分析：該題目考察歐拉方法解初值問題程序如下：function Heun(a,b,y0,n)h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc);endfor k=1:n+1 fprintf('x%d=%fty%d=%fn',

2、k-1,x(k),k-1,y(k);endfunction z=f(xx,yy)z=10*xx*(1-yy);運行結(jié)果：>> Heun(0,1,0,10)x0=0.y0=0.x1=0.y1=0.x2=0.y2=0.x3=0.y3=0.x4=0.y4=0.x5=0.y5=0.x6=0.y6=0.x7=0.y7=0.x8=0.y8=0.x9=0.y9=0.x10=1.y10=0.>>分析： 該結(jié)果與準(zhǔn)確結(jié)果相比比較接近，但是有一定的誤差。6.用四階龍格庫塔公式解第三題中的初值問題，取步長h=0.2，保留五位有效數(shù)字。題目目的分析：該題考查四階龍格-庫塔方法和改進(jìn)歐拉方法求解

3、精確度問題。程序：改進(jìn)歐拉法：function Heun(a,b,y0,n)h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc);endfor k=1:n+1 fprintf('x%d=%fty%d=%fn',k-1,x(k),k-1,y(k);end四階龍格-庫塔方法：function odeRK4(a,b,y0,n)h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);for j=

4、2:n+1 k1=f(x(j-1),y(j-1); k2=f(x(j-1)+h/2,y(j-1)+h/2*k1); k3=f(x(j-1)+h/2,y(j-1)+h/2*k2); k4=f(x(j-1)+h,y(j-1)+h*k3); y(j)=y(j-1)+h/6*(k1+k4)+h/3*(k2+k3);endfor k=1:n+1 fprintf('x%d=%fty%d=%fn',k-1,x(k),k-1,y(k);endfunction z=f(xx,yy)z=-yy;結(jié)果：>> Heun(0,1.0,1,10)x0=0.y0=1.x1=0.y1=0.x2=0

5、.y2=0.x3=0.y3=0.x4=0.y4=0.x5=0.y5=0.x6=0.y6=0.x7=0.y7=0.x8=0.y8=0.x9=0.y9=0.x10=1.y10=0.>>>> odeRK4(0,1.0,1,5)x0=0.y0=1.x1=0.y1=0.x2=0.y2=0.x3=0.y3=0.x4=0.y4=0.x5=1.y5=0.>>結(jié)果分析：由程序結(jié)果得改進(jìn)歐拉法和四階龍格-庫塔法運行結(jié)果上在小數(shù)第二位后才開始有所不同，但結(jié)果相差不大，精確度都是挺高的12.將二階方程y''-5y+6y=0 (y0=1,y'0=-1)化為一階

6、方程組，取h=0.1,用四階龍格-庫塔法求有y（0.2）的近似值，保留5位有效數(shù)字。解：目的分析： 該題目是要將二階方程化為一解方程組，并用四階龍格-庫塔法求近似值。程序如下：公式部分：function z=f1(xx,yy)z=5*xx-6*yy;四階龍格-庫塔法：function odeRK4_1(a,b,y0,n)h=(b-a)/n;x=a:h:b;y=y0*ones(1,n+1);z=-ones(1,n+1);for j=2:n+1 k1=z(j-1); l1=f1(z(j-1),y(j-1); k2=z(j-1)+h/2*l1; l2=f1(z(j-1)+h/2*l1,y(j-1)+h/2*k1); k3=z(j-1)+h/2*l2; l3=f1(z(j-1)+h*l2/2,y(j-1)+h/2*k2); k4=z(j-1)+h*l3; l4=f1(z(j-1)+h*l3,y(j-1)+h*k3); y(j)=y(j-1)+h*(k1+2*k2+2*k3+k4)/6; z(j)=z(j-1)+h*(l1+2*l2+2*l3+l4)/6;endfor k=1:n+1 fprintf('x%d=%fty%d=%fn',k-1,x(k)