9常微分方程的数值解法

第九章常微分方程的数值解法,1,1引言,在许多应用问题中得出的数学模型，常常是一个或一组常微分方程，然而只有极少数的微分方程能够用初等的方法求出它们的解析解，多数微分方程只能用近似的方法求其近似解。近似方法有两类，一类称为近似解析法，如级数解法，逐次逼近法；另一类称之为数值解法。本章主要介绍数值解法，其基本思想是利用近似计算求出微分方程的解在一些离散点处的近似值。,由于高阶的微分方程可以转化为一阶常微分方程组，而一阶常微分方程组又可写成向量形式的单个方程。所以本章主要考虑常微分方程初值问题,2,(9.1),为了使利用数值解法得出的常微分方程的近似解具有实际意义，就必须保证常微分方程的解存在且唯一。,定理9.1设常微分方程(9.1)中的二元函数满足,在区域上连续。,在上D关于y满足Lipschitz条件，即存在常数L，使,(9.2),其中L称为Lipschitz常数。则问题(9.1)在区间a,b上,存在唯一连续可微的解,3,以后我们总是假定f(x,y)在D上关于y满足Lipschitz条件。,数值解法的具体做法是在区间a,b上插入n-1个点，使,记，称hi为从xi到xi+1的步长，在这些节点xi,上通过离散化的方法，求得y(xi)的近似值yi。,数值解法常分为两大类,(1)单步法：计算节点xi+1处的近似值yi+1时只用到了前面,一个点xi处的信息。,(2)多步法：计算节点xi+1处的近似值yi+1时用到了前面,多个点处的信息。,4,定义9.1若对，微分方程的数值解法的计算,公式均准确成立，但至少有一个r+1次的多项式不能准确成立，,则称计算公式是r阶的。,记为准确值，为计算的近似值，称,当公式是r阶时，有,称C为渐进误差常数。,为近似值yi+1的局部截断误差,5,2Euler方法,Euler方法的精度较低，因此在实际计算时，Euler方法没有多大的实用性，但由于Euler方法的推导比较简单，而且能说明一般的数值计算公式构造时的一些技巧及思想，因此我们从Euler法开始介绍常微分方程的数值方法。,一、Euler方法的推导,对微分方程（9.1），将a,bn等分，步长，节点,又设方程（9.1）存在唯一充分连续可微的解y(x)，,现介绍三种推导Euler公式的方法。,6,1Taylor展开,在点xi处将y(xi+1)进行Taylor展开有,省去误差项得到近似公式,局部截断误差,（9.3）,设yi是y(xi)的近似值，则对方程（9.1）有差分方程,7,(9.4),(9.4)称为求解初值问题(9.1)的Euler方法。从(9.3)(9.4)可知，Euler方法是一阶的单步显式法。推导公式(9.4)还可以用其它的方式。,2数值微分方法,其基本思想是用差商代替微商，对充分小的h有,因此有近似公式,8,从而有差分方程（9.4）。,3数值积分法,在上对积分有,对右端的积分式，利用数值积分的左矩形公式有,从而得方程（9.4）,9,二、几何意义,对差分方程（9.4），如图9.1所示，设y=y(x)是（9.1）的解，线段AB为曲线在点A处的切线，差分方程（9.4）的几何意义是利用线段AB来近似代替曲线段AD。局部截断误差,为线段BD，而且BD的值与h的大小和的值有关。,如图9.2所示，Euler方法的几何意义是用折线段,来近似代替方程（9.1）的解,10,三、改进的Euler法,1中点方法,在区间上对积分有,对右端的积分项使用数值积分的中矩形公式有,由此有差分方程,11,（9.5）,局部截断误差,差分方程（9.5）称为中点公式，中点公式（9.5）是二步二阶显式方法。,2梯形方法,在上对积分有,对右端积分项使用数值积分公式中的梯形公式，有,12,由此有差分方程,（9.6）,局部截断误差为,差分方程（9.6）称为梯形公式，它是单步二阶方法。由于（9.6）式的两端均含有yi+1的信息，它是关于未知量yi+1的一个隐函数方程，称之为隐式方法。,13,对梯形公式常采用迭代法来求解yi+1的值，其迭代格式为,（9.7）,定理9.2若f(x,y)在区域,上关于y满足Lipschitz条件（9.2），且，则迭代,公式（9.7）生成的序列收敛于,14,例9.1设有微分方程,分别用Euler法和Euler预测校正法作数值计算，并比较其计算结果（取h=0.1）,这种方法称为Euler预测校正法，计算公式为,称（9.8）为预测式，（9.9）为校正式。,15,解（1）微分方程的解析解为,（2）Euler方法的计算公式为,（3）Euler预测校正法的计算公式为,计算结果如表9.1所示。,16,表9.1,17,3Runge-Kutta方法,Runge-Kutta方法的基本思想是利用函数f(x,y)在某些点上函数值的线性组合来计算y(xi+1)处的近似值yi+1。构造近似公式时要求近似公式在(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前几项重合，从而使近似公式达到所需的阶数，这样既避免了计算函数f(x,y)的偏导数，又提高了精度。,一、Runge-Kutta方法的构造,Runge-Kutta方法的一般形式,(9.10),18,（9.1）的显式单步法，一般每步需要计算p次f(x,y)的值，,其中是待定参数。由（9.10）得,因此对给定的，则（9.10）式为求解微分方程,则Runge-Kutta方法的局部截断误差为,故称之为p级Runge-Kutta方法。假设y=y(xi)是准确的，,（9.11）,若则称方法是m阶的，现以p=2为例来推导,Runge-Kutta方法，当p=2时有,19,（9.12）,将（9.12）中的K2在(xi,yi)处进行Taylor展开,对y(xi+1)在xi处进行Taylor展开有,20,为使局部截断误差Ti的阶数尽量地高，应选择适当的参数,使（9.13）式中右端的h和h2的系数恒为0。,即取,（9.14）,以为自由参数有,（9.15）,此时局部截断误差为,21,对(9.15)式，选取不同的值可得到相应的二级二阶,Runge-Kutta公式，取，得中点公式,取，得Heun公式,取，得Euler预测校正公式,22,二、四阶经典Runge-Kutta公式,（9.17）,（9.17）也称为标准四阶Runge-Kutta方法，它常用来作线性多步法的启动值计算，（9.17）式也可写成,23,例9.2用经典四阶Runge-Kutta方法求解,24,解：取步长h=0.2，计算公式为,25,计算结果如下表所示。,与例9.1的Euler方法计算结果相比较，尽管Runge-Kutta方法的步长放大了一倍，而其数值解的精度还是比Euler法高。,26,Runge-Kutta法简练易于编制程序，而且是单步法。其计算也具有数值稳定的优点。但每一次的计算量都比较大，而且公式的局部截断误差也难以求得。所以从计算量来看，一般来说低精度问题常采用低阶Runge-Kutta公式，高精度的问题采用高阶的Runge-Kutta公式，但四阶以上的Runge-Kutta公式很少被使用。,三、步长的选取,步长的选取在数值解法中非常重要，步长过大，每步计算产生的局部截断误差也较大。若步长取得较小，虽然每步计算的截断误差较小，但在一定的求解范围内需要完成的计算步骤就较多。这不仅增加了计算量，而且还会造成舍入误差的累积。应用中的一种有效措施是在计算的过程中自动调整步长，即变步长技巧。现利用Richardson外推法来构造变步长的技巧。,27,设计算公式是p阶的，从yi出发先取步长为h，经过一步,计算得出的数值解记为，局部截断误差记为,然后将步长折半，取步长为，经过两步计算，从yi出发,计算得出的数值解记为，其局部截断误差为,其中的渐近误差常数与在,上的值有关，但我们可以近似认为,28,故有,以乘（9.19）并减去（9.18）式有,（9.20）,显然以（9.20）的左端值近似yi+1比用,作为yi+1的近似值，其精度要高得多。,如取p=4，有,29,这种技巧与Romberg积分法的思想是一致的。从（9.18）（9.19）又可得到近似值,从（9.21）可知，从的值来选择步长h的大小，若误差,精度为，其方法为,（1）当时，反复加倍步长计算，直到，再以,上一次步长计算所得值作为yi+1。,（2）当时，反复折半步长计算，直到，再以,最后一次计算所得值作为yi+1。,30,从表面上看，判别的工作量是增加了，但当方程的解y(x),变化较大的情况下，总的工作量还是减少了。,外推法也可用来进行误差估计，设方法是p阶的，则有,因此有估计式,（9.22）式的右端常用作的误差估计。,31,4线性多步法,从中点方法(9.5)是二步二阶方法，Euler法是(9.4)是单步一阶方法可知，多步法的局部截断误差比单步法要高。在常微分方程的数值解法中，常用多步法来求解。,一、线性多步法的一般形式,求解微分方程（9.1）的线性多步法，即差分方程,32,恒假设不全为0，称（9.23）式为p+1步法。,计算时需要p+1个值,若，则（9.23）是p+1步的显式方法。,若，则（9.23）是p+1步的隐式方法。,假设是准确值，则,（9.23）式产生的局部截断误差为,33,其中,34,现对（9.31）右端的积分项中的被积函数f(x,y(x)采用插值多项式作逼近，便得出实用中最常用的一种线性多步法Adams方法。,1Adams外推法,由p+1个点,构造f(x,y(x)的插值多项式,35,将（9.32）代入（9.31），记有差分方程,其中,局部截断误差为,36,线性多步法计算公式（9.33）是显格式又称之为Adams-Bashforth公式。,当p=3时有计算公式,（9.34）即常用的Adams四步四阶显式法。,当p=0时，计算公式为,此公式就是Euler公式。,37,2Adams内插法,由p+2个点,作f(x,y(x)的Lagrange插值多项式，并代入（9.31）式右端的,积分项中得差分方程,线性多步法（9.36）是隐格式，又称之为Adams-Moulton公式。,当p=2时，计算公式为,38,局部截断误差为,此即常用的Adams三步四阶隐式方法。,当p=0时，有计算公式,此即Euler梯形公式,实用中常将(9.34)和(9.37)联合起来构成Adams预测校正公式,39,其中称第一式为预测公式，第二式为校正公式，,并用经典的四阶Runge-Kutta作启动值计算。,例9.3利用Adams预测校正公式求解微分方程,解取步长h=0.1，计算公式为,40,计算结果,如表所示。,41,5局部截断误差的估计,一、局部截断误差的估计,当时，线性多步法（9.23）是隐式格式，可写为,（9.40）式常用的求解方法是用迭代公式,42,现取显格式作预测式，隐格式作校正式得到预测校正法为,（9.42）,对（9.42）式，一个常用的估计局部截断误差的方法，称为,Milne估计。假设预测式和校正式都是r阶的，误差常数为,则在y(x)充分可微的条件下有,定理9.4设存在正数M，使则当时，,由迭代公式（9.41）产生的数列收敛于yi+1。,43,（9.43）式减去（9.44）式有,略去高阶无穷小有,由此可得预测式和校正式的局部截断误差估计,44,有了（9.46），（9.47）的近似估计后，将其加在原来的预测式和校正式上，可以改进预测值和校正值的精度，但一般只对预测式这样做，因若对校正式这样做，常会使稳定性变坏，称之为修正预测校正法。估计式（9.47）常用来选择步长，控制局部截断误差。,二、修正的Adams预测校正法,对四阶的Adams方法，由（9.35）和（9.38）式有,45,故有修正Adams预测校正法,预测,修正,校正,局部截断误差估计为,46,6一阶方程组与高阶方程,前面各节介绍的初值问题的数值解法都可以推广到方程组和高阶方程的情形。,一、一阶方程组,设有初值问题,引进向量,47,则（9.48）式可写成向量形式,前面介绍的一切方法，都可以用来求解初值问题（9.49）。,（1）Euler公式,写成分量形式为,写成方程组的形式为,48,（2）经典的Runge-Kutta方法的计算公式,49,写成分量形式为,50,二、高阶方程的情形,对高阶方程，可作变量替换将其转化为一阶方程组，,设有高阶方程,引入变量,则（9.50）可化为,51,再用一阶方程组的数值解法来求解（9.51）。,例9.4写出用Euler预测校正方法求初值问题,的计算公式。,52,解记，则所给初值问题化为,则Euler预测校正法的计算公式为,预测公式,校正公式,53,7收敛性和稳定性,一、收敛性,定义9.3假设yi+1之前函数值是准确的，若用差分方程求出,的解yi+1满足,其中y(xi+1)是xi+1处的准确值，则称差分方程是收敛的。,Euler方法，Runge-Kutta方法是收敛的，以及本章介绍的其它方法也是收敛的。,54,二、稳定性,稳定性即数值稳定性，是指在数值计算过程中误差传播的情况。应用数值方法求解微分方程初值问题时，由于求解过程是按节点逐次递推进行。误差的传播是不可避免的。所以若计算公式不能有效地控制误差的传播，那么误差积累，将使最终的计算结果严重失真。,例9.5分别在h=0.1，h=