实验五-常微分方程数值解.doc

实验五 常微分方程数值解 作者： 日期：11 个人收集整理 勿做商业用途实验五 常微分方程数值解一 欧拉法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;x0=0;y0=1;h=（ba)/n；for j=1:nx1=x0+h; y1=y0-h0。9*y0/(1+2x0);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.0800y1 = 0。9337x1 = 0。1000y1 =0.9192结果分析:欧拉法计算简单，但计算效率并不高，计算精度很低，局部截断误差较大.二 改进欧拉法1 算法说明对于xi,i=0,1,2,，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为:yp=yi+hf(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，nFOR 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=(ba)/n;for j=1：nx1=x0+h; yp=y0h0。9*y0/(1+2*x0)；yc=y0h0.9yp/（1+2*x1);y1=（yp+yc)/2;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）的值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结果分析:计算过程比欧拉法较复杂，但改进欧拉法先用欧拉法求出预报值,再利用公式求出校正值，局部截断误差比欧拉法低了一阶，较大程度地提高了计算精度.三 龙格库塔法1 算法说明对于xi，i=0，1，2，n，取步长h为定值时，有h=xi+1-xi,EURLER法的计算公式为：yi+1yih( K1+ 2K2 +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 TO 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+k3h；y1=y0+h(k1+2k2+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数据分析:龙格-库塔法是应用广泛的高精度单步算法。此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。