第九章-常微方程数值解法

第9章常微分方程数值解法,8-1,第9章,常微分方程数值解法,第9章常微分方程数值解法,8-2,第9章目录,1欧拉（Euler）方法1.1Euler法及其简单改进1.2改进的Euler法2龙格库塔（Runge-kutta）方法2.1龙格-库塔方法的基本思想2.2二阶龙格-库塔公式2.3高阶R-K公式2.4变步长R-K法3线性多步法4一阶方程组与高阶方程初值问题5收敛性与稳定性,第9章常微分方程数值解法,8-3,第8章序,许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。,所谓初值问题，是函数及其必要的导数在积分的起始点为已知的一类问题，一般形式为：,我们先介绍简单的一阶问题：,第9章常微分方程数值解法,8-4,第9章序,由常微分方程的理论可知：上述问题的解唯一存在。常微分方程求解求什么？应求一满足初值问题（81）的解函数y=y(x)，如对下列微分方程：,高等数学中，微分方程求解，如对一阶微分方程：y=f(x,y)是求解解函数y=y(x)，使满足上述方程。但能够求出准确的解析函数y(x)的微分方程是很少的，高数中研究微分方程的求解，是分门别类讨论，对不同类型的微分方程，求解方法不一样，因此，要求解微分方程，首先必须认清类型。,第9章常微分方程数值解法,8-5,微分方程数值解,而常微分方程初值问题的数值解法，是要寻求解函数y(x)在一系列点y(xi)（离散点）:,上y(xi)的近似值yi（i=1,2,n），并且还可由这些（xi,yi)（i=1,2,n）构造插值函数作为近似函数。上述离散点相邻两点间的距离hi=xi-1-xi称为步长，若hi都相等为一定数h,则称为定步长，否则为变步长。,由于在实际问题和科学研究中遇到的微分方程往往很复杂，绝大多数很难，甚至不可能求出解析函数y(x)，因此只能考虑求其数值解。,本章重点讨论如下一阶微分方程：,在此基础上介绍一阶微分方程组与高阶微分方程的数值解法。,第9章常微分方程数值解法,8-6,1欧拉（Euler）法,以Euler法及其改进方法为例,说明常微分方程初值问题数值解法的一般概念，Euler法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其它方法的特点。,Euler法用于求解一阶微分方程初值问题：,第9章常微分方程数值解法,8-7,1.1Euler法及其简单改进,Euler公式为：,由x0出发x1,x2,xN，而利用此式可算出对应的y1,y2,yN，式（8-2）称为差分方程（序列yn满足的方程）。,下面是Euler公式的推导：,一、从几何意义出发：y=f(x,y)的解函数y=y(x)在xoy平面上是一族解曲线，而初值问题则是其中一条积分曲线假定y=y(x)的曲线如图8-1从给定的点P0(x0,y0)出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率y=f(x0,y0)，因此可得切线：（点斜式）,第9章常微分方程数值解法,8-8,Euler公式的推导(续1）,几何意义：用折线近似解曲线y=y(x)，折线不会偏离太远，因为每项以f(x,y)（斜率）修正。,切线与x=x1交于P1(x1,y1)，在x0,x1上以切线,近似曲线，,第9章常微分方程数值解法,8-9,Euler公式的推导（续2）,二、利用Taylor级数：将y(x)在xn处展开：,第9章常微分方程数值解法,8-10,Euler公式的推导（续4）,第9章常微分方程数值解法,8-11,Euler公式的推导（续5）,四、利用数值积分公式：在xn,xn+1上对y(x)=f(x,y(x)积分,对右端积分项采用不同的数值积分公式，便可得到各种不同的求解dE初值问题的计算公式。如，以矩形面积代替曲边梯形面积,1）以左矩形面积代替曲边梯形面积如图8-2，亦即以,图8-2,第9章常微分方程数值解法,8-12,3）以梯形公式（面积）代替曲边形如图8-4则有,式（8-5）称为求dE初值问题的梯形公式，梯形公式看作是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数得到的，而Euler公式则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的解微分方程的数值公式，梯形公式比Euler公式更准确一些，误差更小。,Euler公式的推导（续6）,2）以右矩形面积代替曲边梯形(后退的梯形公式）：如图8-3,第9章常微分方程数值解法,8-13,Euler公式注释,注1：Euler公式为显式，后退的Euler公式，梯形公式为隐式；,注2：Euler公式，梯形，后退的Euler公式为单步法，计算yn+1只用yn，而中点法公式为多步法（还可如上二所述，构造多步法）即必须已知yn-1，yn才能计算yn+1，（求y0,y1不能用此公式。y0,y1称为多步法的开始值，y0给定，而y1必须由其它公式算出，然后才能用中点法）；,注3：前面已有Euler法的局部截断误差：,后退Euler法的局部截断误差：,误差阶：如果局部截断误差,则称方法为P阶的。,第9章常微分方程数值解法,8-14,显然，步长h越小，阶数P越高，局部截断误差越小，当然计算精度越高；,注4：梯形法是几阶？梯形法精度比Euler法高，阶数肯定比Euler法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积分公式计算,因此由积分中梯形公式的误差知此,时的局部截断误差为：,梯形法为2阶方法！,Euler法，后退Euler法为1阶方法，而中点法为2阶，,Euler公式注释（续）,第9章常微分方程数值解法,8-15,关于Euler法的整体截断误差注释,注5：关于Euler法的整体截断误差：,实际计算时，yn是y(xn)的近似值，因此，计算过程中除每步所产生的局部截断误差外，还有因前面的计算不准确而引起的误差。在不考虑舍入误差的情况下，称y(xn+1)与yn+1之差为整体截断误差，记为：,下面讨论Euler方法的整体截断误差。,为简便起见，假定函数f(x,y)充分光滑，问题（8-1）解y(x)在a,b上二阶连续可微，于是由式（8-6），局部截断误差有界，即存在M0，使得对任意xa,b，都有|y(x)|M，从而有：,（紧接下屏）,第9章常微分方程数值解法,8-16,结论,对于实际问题来说，由于L，M难以估计，因很难应用，而且上述推导过程中一再放大了误差上限，这样的估计往往也很保守，远远大于实际的误差，但是，从估计式（8-7）却可以得到下面很有用处的结论。,1）当h0时，en0即，,亦即数值解yn，一致收敛于初值问题（8-1）的真y(xn)，,并且，Euler法的整体截断误差的阶为O(h)与h同阶，比局部截断误差低一阶。,2）舍入误差局部截断误差对实际计算结果有影响，并且随h减少而减少或增大。,第9章常微分方程数值解法,8-17,3）计算结果与解法的阶数p，真解的导数y(p+1)有关，p越大，hp+1越小，|y(p+1)()|的上限越大，M也越大，因此为保证精度当然应选阶数p较高的方法。但如果M很大，当f(x,y)是分段连续的函数时，则应采用低阶的方法如用Euler法。,结论（续）,4）计算结果还与开始值的精度有关，为使这种误差的影响不致于超过局部截断误差，对多步法，应采用跟多步法同阶的方法计算开始值。,第9章常微分方程数值解法,8-18,1.2改进的Euler法,梯形公式为二阶方法，但却是隐式格式，即若利用梯形公式求yn+1，就要求解方程（8-5）式，计算量较大，通常在实际计算时，将Euler法与梯形公式合起来使用，即先使用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为预测值，然后用梯形公式去提高精度，将yn+1(0)校正为较准确的值：,由于函数f(x,y)满足Lipschitz条件，容易得出:,第9章常微分方程数值解法,8-19,改进的Euler法(续）,第9章常微分方程数值解法,8-20,预测校正型公式,实际经验表明，式（8-8）的迭代效果主要体现在第一次，由此构成如下的预测校正型公式:,此式称为改进的Euler公式，为上机计算编程方便，常将式（8-9）改写为:,下面分析改进的Euler公式的局部截断误差：,第9章常微分方程数值解法,8-21,改进的Euler公式的局部截断误差分析,假定yn=y(xn)，y(xn+1)的Taylor展式为:,对于改进的Euler公式，由于,这说明改进的Euler法的局部截断误差为O(h3)，比Euler公式高一阶，是二阶方法。,第9章常微分方程数值解法,8-22,改进的Euler公式举例,例1,这些结果在表8-1中，可见计算结果的精度，Euler法与后退Euler法差不多，与准确值相比较Euler法偏小，而后退Euler法偏大；中点法与梯形法精度同为2阶，但梯形法更好一些，这跟它们局部截断误差的符号，阶数和系数的大小是完全一致的。,表见下屏：,第9章常微分方程数值解法,8-23,表格8-1,表8-1y=y,y(0)=1的数值解(h=0.1),第9章常微分方程数值解法,8-24,表格8-2,而表8-2是分别取了不同的h=0.1,h=0.01，h=0.001,h=0.0001，还是利用这些公式，经过若干步的计算（h越小，计算量越大）算到y(1)的近似值，可见：随着h的减小，y(1)的近似值的精度在提高，0.01比0.001差，即0.001比0.01时的y(1)准确。,（紧接下屏）,第9章常微分方程数值解法,8-25,表8-2计算结果说明（续）,但h太小，到h=0.0001时却又变得误差大了，这与前面所说h越小，p阶越高，应该局部截断误差越小，因而计算精度更高矛盾了，为什么会产生这种情况呢？这是由于h太小而引起计算量大因而造成了舍入误差和截断误差的积累，这种情况由于初值问题不同可能会影响更大，偏离更严重，如下面的例2。这种问题实际上是稳定性问题，我们将会讨论方法的稳定性，由此得出对h有一定的要求的稳定性制区域。,第9章常微分方程数值解法,8-26,5收敛性与稳定性,通过前面的讨论可以看到，微分方程数值解法的基本思想是：通过某种离散化手段，将微分方程转化为差分格式求解。从理论上说，这里需要解决两个问题，一是当步长h0时，差分解，这就是差分格式的收敛性问题；二是利用差分格式求解时，初始误差及计算过程中的舍入误差能否得到控制，这就是差分格式的稳定性问题。本节对这两个问题作一简单介绍。,第9章常微分方程数值解法,8-27,5.1收敛性,用某种差分格式解微分方程时，若对求解区间a,b中的任一点x，,当h0时差分解yn一致收敛于微分方程之解y(x)，,对x一致成立，则称该差分格式是收敛的。,显然，只有收敛的差分格式才有实用价值。可以证明，前几节介绍的单步法与多步法对满足解的存在唯一性定理条件的微分方程都是收敛的。,对于收敛的差分格式，从理论上说只要h充分小，差分解yn的精度就可以任意高。,第9章常微分方程数值解法,8-28,5.2稳定性,实际计算中舍入误差几乎不可避免，这类误差在递推计算过程中会不会恶性增长，以致淹没了差分方程的真解，这是微分方程数值解法中的另一个重问题。由于研究的对象不同，关于稳定性的定义很多，这里只介绍其中一种绝对稳定性。由于实际问题中微分方程（8-1）的右端f(x,y)复杂多样，这给研究稳定性带来了困难，因此通常取以下模型方程为基础进行研究。,其中为常数，0表示微分