数值分析 2.doc

Euler法和预估校正法求解初值问题Euler法和预估校正法求解初值问题摘要在数学与计算科学中，Euler法是一种一阶数值方法，通常用于对给定初值的常微分方程（初值问题）的求解。Euler法的基本思想是迭代，就是逐次替代，然后求出所要求的解，并达到一定的精度。Euler法思想是简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大，因此Euler法一般不用于实际计算。为提高精度，需要在Euler法的基础上进行改进，即为预估校正法。预估校正法的精度为二阶，思想是采用区间两端的函数值的平均值作为直线方程的斜率。预估校正法先用Euler法求出预报值，再利用梯形公式求出校正值，局部截断误差比Euler法低了一阶，较大程度地提高了计算精度。并编写MATLAB程序实现两种数值解法，通过作图对比其精度，加深对两种方法的认识。 关键字：Euler法，预估校正法，MATLAB软件EULER METHOD AND FORECAST CORRECTION METHOD FOR SOLVING INITIAL VALUE PROBLEMSABSTRACTIn mathematics and computer science, the Euler method is a numerical method. It is usually used to solve the equations of the given initial value(initial value problems),Eulers basic method is iterative, that is to say, the ideal is successive substitution, then, find out the required solution and achieved a certain accuracy. Euler method simply means take as the starting point of the next step to calculate the tangent of the end point, when numbers increase, errors due to the accumulation of more and more big. So, the Euler method is generally not used for practical calculation. In order to improve the accuracy, we need to be on the basis of Euler method was improved, the forecast correction method. Forecasts for the second order correction method of the precision, using the average value of a function as a linear equation at each end of the range of the slope. Forecast correction method with Euler method first predicted value, using trapezoid formula to find the correction, the local truncation error lower than the Euler method, greatly improve the calculation accuracy. By write MATLAB program to realize two methods, and through comparing the drawing accuracy, deepen understanding of the two methods.Key words: Euler method, forecast correction method, MATLAB目录1 欧拉法.11.1 Euler方法简介.11.1.1 Euler格式.21.1.2欧拉方法的误差估计.22 预估校正法.62.1预估校正法简介.621.1预估校正法.62.2.2 预估校正法的误差估计.63.实例以及结果分析.4 3.1 Euler法与预估校正法的Matlab实例及实现.3.1.1 实例1的求解及Matlab实现.73.1.2 实例2的求解及Matlab实现.3.1.3实例3的求解及Matlab实现.参考文献.10 附录.11 Euler法和预估校正法求解初值问题Euler法1.1 Euler方法一阶常微分方程的初值问题,其一般形式为 （1）我们知道，只要函数f（x,y）适当光滑-譬如关于y满足Lipschitz条件理论上就可以保证初值问题（1）的解存在且唯一。所谓数值解法，就是求问题（1）在某些离散节点a=x0x1xn=b的近似解y0,y1,y2,yn的方法。y0,y1,y2,yn就称为问题（1）的数值解。hn=xn+1-xn成为xn到xn+1的步长，为了方便取为常量h。1.1.1 Euler格式将微分方程离散化，用向前差商代替微分，带入（1）中的微分方程，得： （n=1,2,3）化简得： （n=1,2,3） （2）如果用近似代入上式便可得的近似值，计算式为： （n=1,2,3）这样问题（1）的近似解可通过求解下面的差分初值问题： （3）得到，按（3）式由初值可逐次求出，。Euler方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的解。即由公式（3）算出的近似值。这组公式求问题（1）的数值解就是著名的欧拉公式。1.1.2欧拉方法的误差估计对于欧拉公式（3）我们看到，当n=1,2时公式右端的都是近似的，所以用它计算的会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。假定用（3）式时右端的没有误差，即=那么由此计算出 （4）局部误差指的是，按（4）式计算由到这一步的计算值与精确值之差。为了估计它，由taylor展开得到的精确值是 （5）（4）、（5）两式相减（注意到）得： （6）即截断误差是阶的，而数值算法的精度定义为：若一种算法的局部截断误差为，则称该算法具有p阶精度。显然p越大精度越高。式（6）说明，欧拉方法式一阶方法，因此它的精度不高。预估校正法2.1预估校正法21.1预估校正法用数值方法离散化问题（1），两端积分得对右端积分使用梯形公式可得用，代替，则得计算公式 （7）很明显可以注意到（7）式为隐式形式。预估校正法是先用欧拉公式求的一个近似解，称为预测值，然后用梯形公式进行矫正并求得近似解，即2.2.2 预估校正法的误差估计由（7）式知，所以预估校正法是二阶的，精度较欧拉法要高，实用性更广泛。实例以及结果分析3.1 Euler法与预估校正法的Matlab实例及实现3.1.1 实例1的求解及Matlab实现设f(t,y)=-r*y，r=-10;初值t=0, y0=2。且y=-10y 的真实解为y(t)=2e-10t。设t0.02，r=-10。解：根据欧拉法和预估校正法的算法编程（程序见附录），然后得出在不同步长下，两算法在不同节点处的值，再与题中给出的解函数值对照，其结果如下欧拉法结果为：预估校正法结果：结果分析：由上述图易见，在同等步长下，预估校正法的误差0.0005欧拉法的误差0.05；而使用同一方法时，步长小的收敛性更好；采用预估校正法，当h=0.1时，预测值与真值已十分接近，所以预估校正法的收敛性比欧拉法好.31.2 实例2的求解及Matlab实现在区间0，1上以h=0.1用欧拉法，预估校正法求解微分方程 dy/dx=-y+x+1，初值y（0）=1；其精确解为y=x+exp(-x),且将计算结果与精确解进行比较，对两个算法的收敛性进行分析比较。解：首先，根据欧拉法和预估校正法的算法编程（程序见附录），然后得出在不同步长下，两算法在不同节点处的值，再与题中给出的解函数值对照，其结果如下（此处仅比较h=0.1时的值）：真值A =1.0000 0.9094 0.8351 0.7742 0.7239 0.6823 0.6476 0.6182 0.5931 0.5712 0.5518欧拉y=1.0000 0.9000 0.8190 0.7535 0.7004 0.6572 0.6218 0.5925 0.5680 0.5472 0.5291预估校正y =1.0000 0.9095 0.8354 0.7745 0.7243 0.6827 0.6480 0.6186 0.5935 0.5716 0.5521再对各值进行图上分析：结果分析：由上述图易见，在同等步长下，预估校正法的误差0.0005 Euler(fun,0,1,0.1,10)T = Columns 1 through 10 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178 Column 11 1.00001.7848 Euler(fun,0,1,0.2,5)T = 0 0.2000 0.4000 0.6000 0.8000 1.0000 1.0000 1.2000 1.3733 1.5315 1.6811 1.8269预估校正法：M文件：Euler_modify.mfunction E=Euler_modify(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n); y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n)+feval(fun,x(n+1),z0);endT=x;y Euler_modify(fun,0,1,0.1,10)T = Columns 1 through 9 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 1.00