实验五常微方程数值解.docVIP

实验五 常微分方程数值解 欧拉法 算法说明 对于xi，i=0,1,2，…，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为：yi+1=yi+h*f(xi,yi)；i=0,1,2，…，n。 程序中主要符号说明 a为x的下界， b为x的上界， h为步长， n为循环次数（即为x的数值点数减一）， x0、y0为循环初始值， x1、y1为输出值。 计算流程框图 开始 开始 读入数据x0，y0，h，n FOR I=1 TO N X0+H=X1,Y0+H*F(X0,Y0)=Y1 输出X1，Y1 X1=X0 Y1=Y0 NEXT 结束 程序清单 syms a b h n x0 y0 x1 y1; a=0; b=0.1; n=5; x0=0; y0=1; h=(b-a)/n; for j=1:n x1=x0+h; y1=y0-h*0.9*y0/(1+2*x0); x1 y1 x0=x1; y0=y1; end 算例及输入数据说明 算例：求解初值问题 当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)的值 程序运行结果及结果分析 运行结果： x1 = 0.0200 y1 = 0.9820 x1 = 0.0400 y1 = 0.9650 x1 = 0.0600 y1 = 0.9489 x1 = 0.0800 y1 = 0.9337 x1 = 0.1000 y1 = 0.9192 结果分析：欧拉法计算简单，但计算效率并不高，计算精度很低，局部截断误差较大。 改进欧拉法 算法说明 对于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。 程序中主要符号说明 a为x的下界， b为x的上界， h为步长， n为循环次数（即为x的数值点数减一）， x0、y0为循环初始值， yp、yc为运算中间值， x1、y1为输出值。 计算流程框图 开始 开始 读入数据x0，y0，h，n FOR I=1 TO N X0+H=X1 Y0+H*F(X0,Y0)=YP Y0+H*F(X0,Y0)=YC (YP+YC)/2=Y1 输出X1，Y1 X1=X0 Y1=Y0 NEXT 结束 程序清单 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:n x1=x0+h; yp=y0-h*0.9*y0/(1+2*x0); yc=y0-h*0.9*yp/(1+2*x1); y1=(yp+yc)/2; x1 y1 x0=x1; y0=y1; end 算例及输入数据说明 算例：求解初值问题 当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)的值决定 程序运行结果及结果分析 x1 = 0.0200 y1 = 0.9825 x1 = 0.0400 y1 = 0.9660 x1 = 0.0600 y1 = 0.9503 x1 = 0.0800 y1 = 0.9354 x1 = 0.1000 y1 = 0.9212 结果分析： 计算过程比欧拉法较复杂，但改进欧拉法先用欧拉法求出预报值，再利用公式求出校正值，局部截断误差比欧拉法低了一阶，较大程度地提高了计算精度。 龙格库塔法 算法说明 对于xi，i=0,1,2，…，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为：yi+1＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6 　　K1=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。 程序中主要符号说明 h为步长， n为循环次数（即为x的数值点数减一）， x0、y0为循环初始值， k