资源描述:
《实验五常微分方程数值解》由会员上传分享,免费在线阅读,更多相关内容在应用文档-天天文库。
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,nFORI=1TONX0+H=>X1,Y0+H*F(X0,Y0)=Y1输出X1,Y1X1=>X0Y1=>Y0NEXT结束4.程序清单symsabhnx0y
2、0x1y1;a=0;b=0.1;n=5;x0=0;y0=1;h=(b-a)/n;forj=1:nx1=x0+h;y1=y0-h*0.9*y0/(1+2*x0);x1y1x0=x1;y0=y1;end1.算例及输入数据说明算例:求解初值问题当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)的值2.程序运行结果及结果分析运行结果:x1=0.0200y1=0.9820x1=0.0400y1=0.9650x1
3、=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+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
4、为循环初始值,yp、yc为运算中间值,x1、y1为输出值。3.计算流程框图开始读入数据x0,y0,h,nFORI=1TONX0+H=>X1Y0+H*F(X0,Y0)=>YPY0+H*F(X0,Y0)=>YC(YP+YC)/2=>Y1输出X1,Y1X1=>X0Y1=>Y0NEXT结束1.程序清单symsabhnx0y0ypycx1y1;a=0;b=0.1;n=5;x0=0;y0=1;yp=1;yc=1;h=(b-a)/n;forj=1:nx1=x0+h;yp=y0-h*0.9*y0/(1+2*x0);yc=y0-h*0.9*yp/(1+2*x1
5、);y1=(yp+yc)/2;x1y1x0=x1;y0=y1;end1.算例及输入数据说明算例:求解初值问题当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)的值决定2.程序运行结果及结果分析x1=0.0200y1=0.9825x1=0.0400y1=0.9660x1=0.0600y1=0.9503x1=0.0800y1=0.9354x1=0.100
6、0y1=0.9212结果分析:计算过程比欧拉法较复杂,但改进欧拉法先用欧拉法求出预报值,再利用公式求出校正值,局部截断误差比欧拉法低了一阶,较大程度地提高了计算精度。一.龙格库塔法1.算法说明对于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。2.程序
7、中主要符号说明h为步长,n为循环次数(即为x的数值点数减一),x0、y0为循环初始值,k1,k2,k3,k4为运算中间值,x1、y1为输出值。3.计算流程框图开始读入数据x0,y0,h,nFORI=1TONX0+H=>X1F(X0,Y0)=>K1F(X0+H/2,Y0+K1*H/2)=>K2F(X0+H/2,Y0+K2*H/2)=>K3F(X0+H,Y0+K3*H)=>K4Y0+(K1+2K2+2K3+K4)=>Y1输出X1,Y1X1=>X0Y1=>Y0NEXT结束1.程序清单symshnx0y0x1y1k1k2k3k4;x0=0;y0=1;
8、h=0.2;n=5;forj=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