《实验五 常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《实验五 常微分方程数值解.docx（7页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实验五 常微分方程数值解一 欧拉法1 算法说明对于xi，i=0,1,2，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为：yi+1=yi+h*f(xi,yi)；i=0,1,2，n。2 程序中主要符号说明a为x的下界，b为x的上界，h为步长，n为循环次数（即为x的数值点数减一），x0、y0为循环初始值，x1、y1为输出值。3 计算流程框图开始读入数据x0，y0，h，nFOR I=1 TO NX0+H=X1,Y0+H*F(X0,Y0)=Y1输出X1，Y1 X1=X0 Y1=Y0NEXT结束4 程序清单syms a b h n x0 y0 x1 y1;a=0;b=0.1;n=5
2、;x0=0;y0=1;h=(b-a)/n;for j=1:nx1=x0+h; y1=y0-h*0.9*y0/(1+2*x0);x1y1x0=x1;y0=y1;end5 算例及输入数据说明算例：求解初值问题当x=0,0.02,0.04,0.10时的数值解。输入数据说明：a=0; a为x 的下界b=0.1; b为x的上界n=5; n为循环次数，即为x的数值点数减一x0=0; x（0）的值y0=1; y(0)的值6 程序运行结果及结果分析运行结果：x1 = 0.0200y1 = 0.9820x1 = 0.0400y1 = 0.9650x1 = 0.0600y1 = 0.9489x1 = 0.0800
3、y1 = 0.9337x1 = 0.1000y1 =0.9192结果分析：欧拉法计算简单，但计算效率并不高，计算精度很低，局部截断误差较大。二 改进欧拉法1 算法说明对于xi，i=0,1,2，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为：yp=yi+h*f(xi,yi) yc=yi+h*f(xi+1,yp) yi+1=(yp+yc)/2；i=0,1,2，n。2 程序中主要符号说明a为x的下界，b为x的上界，h为步长，n为循环次数（即为x的数值点数减一），x0、y0为循环初始值，yp、yc为运算中间值，x1、y1为输出值。3 计算流程框图开始读入数据x0，y0，h，nF
4、OR I=1 TO NX0+H=X1 Y0+H*F(X0,Y0)=YPY0+H*F(X0,Y0)=YC (YP+YC)/2=Y1输出X1，Y1 X1=X0 Y1=Y0NEXT结束4 程序清单syms a b h n x0 y0 yp yc x1 y1;a=0;b=0.1;n=5;x0=0;y0=1;yp=1;yc=1;h=(b-a)/n;for j=1:nx1=x0+h; yp=y0-h*0.9*y0/(1+2*x0);yc=y0-h*0.9*yp/(1+2*x1);y1=(yp+yc)/2;x1y1x0=x1;y0=y1;end5 算例及输入数据说明算例：求解初值问题当x=0,0.02,0.
5、04,0.10时的数值解。输入数据说明：a=0; a为x的下界b=0.1; b为x的上界n=5; n为循环次数，即为x的数值点数减一x0=0; x(0)的值y0=1; y(0)的值yp=1; 由y(0)的值决定yc=1; 由y(0)的值决定6 程序运行结果及结果分析x1 = 0.0200y1 = 0.9825x1 = 0.0400y1 = 0.9660x1 = 0.0600y1 = 0.9503x1 = 0.0800y1 = 0.9354x1 = 0.1000y1 =0.9212结果分析：计算过程比欧拉法较复杂，但改进欧拉法先用欧拉法求出预报值，再利用公式求出校正值，局部截断误差比欧拉法低了一
6、阶，较大程度地提高了计算精度。三 龙格库塔法1 算法说明对于xi，i=0,1,2，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为：yi+1yih*( K1+ 2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)； i=0,1,2，n。2 程序中主要符号说明h为步长，n为循环次数（即为x的数值点数减一），x0、y0为循环初始值，k1,k2,k3,k4为运算中间值，x1、y1为输出值。3 计算流程框图开始读入数据x0，y0，h，nFOR I=1 T
7、O NX0+H=X1 F(X0,Y0)=K1F(X0+H/2,Y0+K1*H/2)=K2 F(X0+H/2,Y0+K2*H/2)=K3 F(X0+H,Y0+K3*H)=K4 Y0+(K1+2K2+2K3+K4)=Y1输出X1，Y1 X1=X0 Y1=Y0NEXT结束4 程序清单syms h n x0 y0 x1 y1 k1 k2 k3 k4;x0=0;y0=1;h=0.2;n=5;for j=1:nx1=x0+h; k1=x0+y0;k2=x0+h/2+y0+k1*h/2;k3=x0+h/2+y0+k2*h/2;k4=x0+h+y0+k3*h;y1=y0+h*(k1+2*k2+2*k3+k4)/6;x1y1x0=x1;y0=y1;end5 输入数据说明算例：求解初值问题输入数据说明：x0=0; x（0）的值y0=1; y(0)的值h=0.2; 步长为0.2n=5; n为循环次数，即为x的数值点数减一6 程序运行结果及结果分析x1 = 0.2000y1 = 1.2428x1 = 0.4000y1 = 1.5836x1 = 0.6000y1 = 2.0442x1 = 0.8000y1 = 2.6510x1 = 1y1 = 3.4365数据分析：龙格-库塔法是应用广泛的高精度单步算法。此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。