《计算方法上机作业(共6页).doc》由会员分享，可在线阅读，更多相关《计算方法上机作业(共6页).doc（6页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上计算方法第四次上机报告2.用欧拉方法解初值y=10x(1-y) 0=x 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.分析： 该结果与准确结果相比比较接近，但是有一定的误差。6.用四阶龙格库塔公式解第三题中的初值问题，取步长h=0.2，保留五位有效数字。题目目的分析：该题考查四阶龙格-库塔方法和改进欧拉方法求解精确度问题。程序：改进欧拉法：function Heun
2、(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=2:n+1 k1=f(x(j-1),y(j-1); k2=f(x(j-1)+
3、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;结果： Heun(0,1.0,1,10)x0=0.y0=1.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.
4、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.结果分析：由程序结果得改进欧拉法和四阶龙格-库塔法运行结果上在小数第二位后才开始有所不同，但结果相差不大，精确度都是挺高的12.将二阶方程化为一阶方程组，取h=0.1,用四阶龙格-库塔法求有y（0.2）的近似值，保留5位有效数字。解：目的分析： 该题目是要将二阶方程化为一解方程组，并用四阶龙格-库塔法求近似值。程序如下：公式部分：function z=f1(xx,yy)z=5*
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),k-1,y(k);endend结果图：结果分析：由程序结果得四阶龙格-库塔法所得y(0.2)=0.所给初值问题精确解为,y(0.2)=0.两结果在数值上相差不大，精确度高。专心-专注-专业