第六章-常微分方程--更正版

下午5时44分41秒,1,下午5时44分41秒,1,常微分方程数值解法,第六章,1、引言2、欧拉法、梯形法和改进欧拉法3、龙格库塔法4、多步法5、Adam法6、Gear法7、精度和误差分析8、数值稳定性9、刚性系统,内容提要,1、应具有熟练的数学推导能力（重在推导、运算而不在于证明）（1）建立数学模型（2）简化数学模型2、应熟悉常用的数值计算方法（1）线性方程组的数值解法（2）非线性方程组的数值解法（3）曲线的拟合与插值（4）微分方程数值解法（5）常用的优化方法（6）代数特征值问题的数值解法,工科学生对数学的要求,（1）稀疏矩阵技术（2）基于稀疏矩阵技术的高阶线性方程组的数值解法（3）非线性方程组的数值解法（4）微分代数方程组的数值解法（5）与最小二乘问题为核心的优化方法（6）代数特征值问题的数值解法,电气学生应特别重视的数值计算方法,对于一个常微分方程：,通常会有无穷个解。如：,因此，我们要加入一个限定条件。通常会在端点出给出，如下面的初值问题：,本课程我们仅仅学习常微分方程的数值解法。所研究的常微分方程的形式为：其中，x(t)是随时间而变化的状态变量，依赖于初值x0，而这种微分方程的求解问题称为常微分方程的初值问题。对于非线性微分方程组，较难通过解析方法求出它的解，通常只能通过数值计算的方法来求解。,常微分方程数值解问题的提出,对式(1)进行数值求解的过程，就是根据t0时刻的初始值x0，依次计算t1时刻x(t1)的近似值x1，t2时刻x(t2)的近似值x2。其中相邻时间的间隔被称为步长，通常在整个积分区域t(t0,tN)，步长hn+1=tn+1-tn都被取定值。基本的算法就是从tn时已知的xn、xn-1和f(xn,tn)、f(xn-1,tn-1)推出tn+1时的值。,数值解的基本做法,有一在太空中的火箭绕一天体高速飞行，假定飞行可以认为在一个平面内，求其运动的规律。如图1所示。,图1,微分方程数值解的实例（1）,根据牛顿万有引力定律，可列出四个微分方程：,微分方程数值解的实例（2）,图2,如果已知某一时刻t四个变量x，y，vx，vy的值，如图所示。,微分方程数值解的实例（3）,在计算机中数值计算火箭运动的程序流程图如图3所示。,图3,任何实用的数值算法都必须满足以下的标准：1.数值计算的精确度2.数值计算的稳定性3.数值计算的效率数值计算的精确度是指每一步数值计算的误差都是有界的。其中整体误差x(tn)-xn数值计算的稳定性是指每一步数值计算产生的误差不至于影响到以后的计算。数值计算的效率则与计算量和步长大小有关。,微分方程数值算法的选择准则,第2节欧拉法、梯形法和改进欧拉法,用表示式（1）的精确解，将在t=tn点泰勒展开，并计算级数在t=tn+1时的值，可得下式：,函数的泰勒级数展开（1）,15,公式称为的n阶泰勒公式.,公式称为n阶泰勒公式的拉格朗日余项.,泰勒中值定理:,阶的导数,时,有,其中,则当,如果时间步长h=tn+1-tn,则由式（1）可知，所以如果高次项非常小，则可用由式（2）等式右边计算出来的xn+1来x(tn+1)作为的近似值。,函数的泰勒级数展开（2）,通常，泰勒级数法可以表示为式中其中整数p称为阶。对于较大的p，用泰勒级数法可以非常精确，但计算效率却不高。,函数的泰勒级数展开（3）,当p=1时，泰勒级数法变为：式(4)称为欧拉法。实际上，引言中的例子采用的就是欧拉法。,欧拉法,欧拉法的几何意义,根据数值积分的左矩形公式，有因此，有,欧拉法的数值积分推导,试用欧拉法计算下列初值问题。取步长h0.5，并与精确值比较。,欧拉法算例（1）,解：由欧拉法得：有：,如果计算xn+1时，所取的斜率不是tn点上的导数f(xn,tn)，而是tn+1点上的导数f(xn+1,tn+1)，就得到后退欧拉法以后会说明，后退欧拉法比欧拉法具有好得多的数值稳定性。,后退欧拉法,根据数值积分的右矩形公式，有因此，有,后退欧拉法的数值积分推导,试用后退欧拉法计算下列初值问题取步长h0.5，并与精确值比较。,后退欧拉法算例,梯形法,如果计算xn+1时，所取的斜率不是tn点上的导数f(xn,tn)，而是tn点上的导数f(xn,tn)和tn+1点上的导数f(xn+1,tn+1)的平均值，就得到梯形法,梯形法,梯形法的数值积分推导,根据数值积分的梯形公式，有因此,试用梯形法计算下列初值问题。取步长h0.5，并与精确值比较。,梯形法算例,解：由梯形法得,在梯形公式的右端中包含有未知的xn+1，这类数值方法称为隐式方法，一般情况下不能直接求解上述方程，而需要采用迭代的方法来求解。,改进欧拉公式,一种简单的做法是先用欧拉法计算出xn+1的近似值，然后将这个近似值再代入到梯形公式中，即采用如下的格式：这种格式就称为改进的欧拉公式，也叫预报校正格式。,改进欧拉公式（2）,试用改进的欧拉公式计算下列初值问题。取步长h0.5，并与精确值比较。,改进的欧拉公式算例（1）,解：,用Euler方法求解问题取h=0.1解设f(x,y)=-y+x+1，x0=0,y0=1,xn=x0+nh=0.1n（n=0,1,10）Euler格式为由y0=1出发，按上面公式的计算结果如表所示,xnyn01.0000000.11.0000000.21.0100000.31.0290000.41.0561000.51.0904900.61.1314410.71.1782970.81.2304670.91.2874201.01.348678,1欧拉方法/*EulersMethod*/,欧拉公式：,亦称为欧拉折线法/*Eulerspolygonalarcmethod*/,欧拉公式的改进：,隐式欧拉法/*implicitEulermethod*/,一般先用显式计算一个初值，再迭代求解。,隐式欧拉法的局部截断误差：,即隐式欧拉公式具有1阶精度。,梯形公式/*trapezoidformula*/,显、隐式两种算法的平均,注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。,改进欧拉法/*modifiedEulersmethod*/,注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,预报公式,校正公式,1EulersMethod,简单,精度低,稳定性最好,精度低,计算量大,精度提高,计算量大,精度提高,显式,多一个初值,可能影响精度,第3节龙格库塔法,龙格库塔法的基本思想（1）,考察差商由微分中值定理，可得因此微分方程的数值解为,龙格库塔法的基本思想（2）,这里的称为区间（tn,tn+1)上的平均斜率，记作k*,即因此只要对平均斜率k*提供一种算法，由式（6）便相应地得到一种微分方程数值计算公式。用这种观点很容易看出欧拉公式、后退欧拉公式、改进欧拉公式的特点。,改进欧拉公式利用了(tn,tn+1)两个点的斜率值取算术平均作为平均斜率k*的近似值其中k2是通过已知信息xn利用欧拉公式预报的。,龙格库塔法的基本思想（3）,2龙格-库塔法/*Runge-KuttaMethod*/,建立高精度的单步递推格式。,单步递推法的基本思想是从(ti,xi)点出发，以某一斜率沿直线达到(ti+1,xi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,首先推广改进欧拉公式。随意考察区间（tn,tn+1)内一点我们希望用tn,tn+p两个点的斜率值k1和k2加权平均得到平均斜率k*，即令这里的为待定常数。,二阶龙格库塔公式的推导（1）,二阶龙格库塔公式的推导（2）,仿照改进欧拉算法。设计算法如下：分别将k1和k2Taylor展开，可得,二阶龙格库塔公式的推导（3）,x（tn+1）二阶Taylor展开式为,二阶龙格库塔公式的推导（4）,比较系数可以发现，要使得两者具有同样的局部截断误差，需要满足,称满足这一条件的所有数值计算格式为二阶龙格库塔公式,二阶龙格库塔公式的推导（5）,当取p=1,=1/2时，所得的公式即为改进欧拉公式。当取p=1/2,=1时，所得的公式称为变型的欧拉公式或中点格式，为,2020/5/30,龙格-库塔(R-K)法的基本思想,Euler公式可改写成,则xn+1的表达式与x(tn+1)的Taylor展开式的前两项完全相同,即局部截断误差为O(h2)。,Runge-Kutta方法是一种高精度的单步法,简称R-K法,同理，改进Euler公式可改写成,上述两组公式在形式上共同点:都是用f(x,t)在某些点上值的线性组合得出x(tn+1)的近似值xn+1,且增加计算的次数f(x,t)的次数,可提高截断误差的阶。如欧拉法:每步计算一次f(x,t)的值,为一阶方法。改进欧拉法需计算两次f(x,t)的值，为二阶方法。,局部截断误差为O(h3),于是可考虑用函数f(x,t)在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式在(tn,xn)处的Taylor展开式与解x(t)在tn处的Taylor展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求高阶导数,又提高了计算方法精度的阶数。或者说,在tn,tn+1这一步内多计算几个点的斜率值，然后将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格库塔（Runge-Kutta）法的基本思想。,一般龙格库塔方法的形式为,其中ai,bij,ci为待定参数，要求上式xn+1在点(tn,xn)处作Tailor展开，通过相同项的系数确定参数。,称为P阶龙格库塔方法。,二阶龙格库塔法在tn,tn+1上取两点tn和tn+a2=tn+a2h,以该两点处的斜率值K1和K2的加权平均(或称为线性组合)来求取平均斜率k*的近似值K，即,式中:K1为xi点处的切线斜率值K1=hf(xn,tn)=hx(tn)K2为tn+a2h点处的切线斜率值,比照改进的欧拉法,将tn+a2视为tn+1，即可得,确定系数c1、c2、a2、b21，可得到有2阶精度的算法格式,2020/5/30,58,这里有4个未知数，3个方程。,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,令,对应项的系数相等，得到,2020/5/30,59,二级R-K方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2,b21，可使每步计算两次函数值的二阶R-K方法达到二阶精度。能否在计算函数值次数不变的情况下,通过选择不同的参数值,使得二阶R-K方法的精度再提高呢?答案是否定的！无论四个参数怎样选择,都不能使公式的局部截断误差提高到三阶。这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶。若要获得更高阶得数值方法,就必须增加计算函数值的次数。,三阶龙格库塔法,60,为进一步提高精度，在区间tn,tn+1上除两点tn和tn+a2=tn+a2h,以外，再增加一点tn+a3=tn+a3h，用这三点处的斜率值K1、K2和K3的加权平均得出平均斜率K*的近似值K，这时计算格式具有形式:,同理使局部截断误差达到O(h4)，使对应项的系数相等，得到系数方程组：,参数的选择不唯一，从而构成一类不同的三阶R-K公式，下面给出一种常用的三阶R-K公式，形似simpson公式：,2020/5/30,62,2020/5/30,63,四阶(经典)龙格库塔法,如果需要再提高精度，用类似上述的处理方法，只需在区间xi,xi+1上用四个点处的斜率加权平均作为平均斜率K*的近似值，构成一系列四阶龙格库塔公式。具有四阶精度，即局部截断误差是O(h5)。推导过程与前面类似，由于过程复杂，这里从略，只介绍最常用的一种四阶经典龙格库塔公式。,四阶龙格库塔公式,龙格库塔法常用的是四阶公式：其中K4的计算公式为：,由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,例1.使用高阶R-K方法计算初值问题,解:,(1)使用三阶R-K方法,67,其余结果如下:,(2)如果使用四阶R-K方法,ixik1k2k3yi1.00000.10000.10000.11030.12561.11112.00000.20000.12350.13760.15951.24993.00000.30000.15620.17640.20921.42844.00000.40000.20400.23420.28661.66645.00000.50000.27770.32590.41631.9993,2020/5/30,其余结果如下:,ixik1k2k3k4yi1.00000.10000.10000.11030.11130.12351.11112.00000.20000.12350.13760.13920.15631.25003.00000.30000.15620.17640.17910.20421.42864.00000.40000.20400.23420.23890.27811.66675.00000.50000.27770.32590.33480.40062.0000,试用四阶龙格库塔公式计算下列初值问题。取步长h1，并与精确值比较。,四阶龙格库塔公式算例,解：由四阶龙格库塔公式得：,odefun=(t,x)(1-2*x*t/(1+t*t),t,y=ode45(odefun,0:2,0),试用数值方法计算分别应用欧拉法，后退欧拉法，梯形法和龙格库塔法。解：将上式转换为常系数微分方程组形式令，则,综合算例,容易看出其解析解是我们将用这一精确解与以下数值计算解作比较。,运用欧拉法得或用矩阵形式表示为：,综合算例欧拉法,或用矩阵形式表示：,综合算例后退欧拉法,或用矩阵形式表示为,综合算例梯形法,运用龙格库塔法可得：,综合算例龙格库塔法,其中精确解为cos(t),综合算例结果,从上图中可以看出用梯形法和龙格库塔法计算得到的结果与精确解几乎无法区分，而欧拉法和后退欧拉法由于是一种一阶方法，则不如上述的两种高阶方法精确。由欧拉法得出的数值计算结果的数值要比精确值稍大，并且随着时间的推移不断变大；后退欧拉法的计算结果则比精确值稍小，随着时间的推移数值不断变小。,结果比较,可以看出欧拉法和后退欧拉法的误差在数值上是相等的，但方向相反。,各种计算方法的整体误差,梯形法和龙格库塔法的数值误差经过放大后的图,梯形法和龙格库塔法的数值误差,第4节精确度与误差分析,数值计算方法的精确度受到主要两个因素的影响：计算机的舍入误差和截断误差。计算机只能对有限位的数进行运算，一般数值必须进行舍入，产生舍入误差，这种误差比较难以消除。截断误差是指计算时用近似表达式去替代精确表达式时产生的误差。在数值积分算法执行中，最有效的方法是用最少的计算量得到最准确的结果。通常，越高阶的方法能得到越精确的解，但也需要更多的计算量。为了减少计算量，就需要在允许范围内选取尽可能大的步长，但步长的选择又受到一些因素影响，如数值积分中每一步所产生的误差。,精确度与误差分析（1）,局部截断误差（localtruncationerror,LTE),局部截断误差，可用下式表示：式中x(tn+1)表示tn+1时的精确值，而xn+1表示数值计算的近似值。局部截断误差假定误差是在当前步长产生的，即认为x(tn)=xn。,试找出欧拉法、后退欧拉法和梯形法的局部截断误差表达式。解：欧拉法欧拉算法的表达式为：如果xn=x(tn)，则,局部截断误差例子（1）,且那么,后退欧拉法后退欧拉法的表达式为：同样可计算得,可以看出欧拉法和后退欧拉法的局部截断误差在数值上是相等的，但方向相反。梯形法,计算得：,一阶的欧拉法与后退欧拉法产生的误差与h2有关；二阶的方法，如梯形法产生的误差与h3有关。,Matlab常用ODE,第5节多步法,另一种求取式（1）近似解x(t)的方法是用一k阶的多项式去近似x(t)：式中系数0，1，k为常数。多步法是指xn+1近似地用xn,xn-1和f(xn,tn),f(xn-1,tn-1)来表达，而不像单步法只用到前一步长的数据。多步法的通式为,基本思路（1）,将上述两式联系起来，由于k阶多项式有k+1个系数(0,k),而数值积分方法有2p+3个系数，所以两式的系数需要满足下式：用一组线性空间中的基函数1(t),2(t),k(t)去确定系数:令j(t)=tj,j=0,1,k则多步法方程的形式变为：,基本思路（2）,令p=0,k=1，显然满足式（7），可得：其导数为由于多步法方程为,基本思路（3）,基本思路（4）,将2个基函数代入上述方程，得到如下2个方程：利用两个基函数的已知特性可得：,由此可得a0=1，由于tn+1-tn=h，所以由式有因此b-1和b0可任意选择。如果令a0=1，b-10，b0=1。则可得到欧拉法：如果令a0=1，b-11，b0=0，得后退欧拉法：,基本思路（5）,注意到在后退欧拉法中b-10，xn+1的表达式中有f(xn+1,tn+1)项，不能直接计算xn+1，此方法称为隐式法；而当b-10时，xn+1的表达式中无f(xn+1,tn+1)项，能直接计算xn+1，此方法称为显式法。,基本思路（6）,现令p=0，k=2，此时2p+3=k+1，所以系数可以唯一确定。令2(t)=t2,2(t)=2t,结合式前面的推导，可得如下方程：,基本思路（7）,如取tn=0，则tn+1=h，可得a0=1，b-11/2，b0=1/2，所以有上式所表示的二阶积分方法就是梯形法，它也是一种隐式方法。,基本思路（8）,第6节Adam法,上面曾讨论过多步法的通用表达式为：如果x(t)是一小于或等于k阶的多项式，则必满足以下约束条件：此时数值计算的多步算法可以给出xn+1的精确解。,基本思路（1）,Adam方法是令系数a1=a2=ap=0，由式（8）得a0=1，则多步法的公式变为：其中p=k-1。之后可选择显式法还是隐式法。显式Adam方法，又称“Adams-Bashforth”法，是令b-10，由式（9）得,显式Adam方法（1）,式（10）也可写成矩阵形式：选择所希望的阶数，即可通过式（11）计算系数bi。,显式Adam方法（2）,试导出三阶Adams-Bashforth法解：令k=3，有得到,显式Adam方法例子（1）,显式Adam方法例子（2）,从而，得到三阶Adams-Bashforth法为：注意在算法的执行过程中，内存中要保留xn，xn1，xn2的数据。,当令b-10，p=(k-2)时，Adam方法称为隐式Adam方法，又称“Adams-moulton”法。此时有：由式（9）得：,隐式Adam方法（1）,隐式Adam方法（2）,写成矩阵形式为：,试推导三阶Adams-moulton法解：令k=3，有计算得到,隐式Adam方法例子（1）,隐式Adam方法例子（2）,从而，三阶Adams-moulton法为：注意在算法的执行过程中，内存中要保留xn，xn-1的数据。并且需要采用迭代解法，因为是隐式的。避免迭代的一般性做法是采用预测校正方法，即用显式Adam做预测，用隐式Adam方法做校正。,第7节Gear法,另一著名的多步法为Gear法，与Adam法令除a0以外的ai都为零相反，Gear法除b-1外所有bi都为零，因为b-10，所以所有的Gear法都为隐式法。k阶Gear公式的推导为：令p=k-1和b0=b1=0，得k阶的Gear法：,基本思路（1）,与Adam法类似，式（15）中的k+1个系数可由下式得出：,基本思路（2）,Gear法例子（1）,试推导三阶的Gear法。解：令k=3，由式（16）得：计算得到,Gear法例子（2）,从而求得三阶Gear法为：注意到在算法的执行过程中，内存中要保留xn，xn1，xn2的数据，并且需要采用迭代解法，因为是隐式的。,（1）采用高阶单步法起动；（2）依次采用多步法。,关于多步法的起动过程,第8节数值稳定性,由上节的讨论可知，步长大小的选择会直接影响到局部截断误差，本节则要介绍步长大小的选择对数值计算方法稳定性的影响。数值计算方法稳定性：任一步产生的误差在以后均能逐步衰减，则称这种方法是稳定的。对于稳定的数值计算方法，积分步长的选择将只决定于局部截断误差。,问题的提出,为了分析步长大小的选择对数值计算方法稳定性的影响，引入下面简单的测试方程：容易得出其解为当0时，x(t)随着时间的推移而趋向无穷大。下面将讨论应用欧拉法、后退欧拉法、多步法数值积分计算时的稳定区域。,数值稳定性分析的方法测试方程法,对测试方程应用欧拉法：有可以看出只有当|1+h|1时，才能使当0时，x(t)随着时间的推移而趋向于零。,欧拉法的数值稳定性分析（1）,所以对于0，系统稳定的条件是h落在以（1，0）为圆心的单位圆内，如图所示。可以得出值越大，步长就必须选择得越小。,欧拉法的数值稳定性分析（2）,相似地，对测试方程使用后退欧拉法：有,后退欧拉法的数值稳定性分析（1）,后退欧拉法的数值稳定性分析（2）,因为当1。即对于0，系统稳定的条件是h落在以（1，0）为圆心的单位圆外，如图所示。可以看出使用后退欧拉法时积分步长可以取得任意大，而不至于影响到解的稳定性，这时步长的选择只取决于局部截断误差。,多步法的数值稳定性分析（1）,对测试方程应用多步法：整理上式得,多步法的数值稳定性分析（2）,相应的特征方程为设z1,z2,zp+1为上式的根，则有：,当0时，xn+1应随着n趋向于无穷而趋向于零，只有当|zj|1（j=1,2,p+1）时才能满足。定义绝对稳定性为整体误差随着n的增大而减少。称多步法是绝对稳定的，如果对于一给定的h，方程的根满足|zi|1（i=1,2,k）。下面将推导多步法的稳定区域。,多步法的数值稳定性分析（3）,多步法的数值稳定性分析（4）,令其中而由于z为复数，其可表示为,最终有多步法的绝对稳定区域的边界可由上式的函数在复平面中绘出，其中h()是从0到2。,多步法的数值稳定性分析（5）,绘出三阶的Gear法和三阶Adam法（隐式和显式）的绝对稳定区域。解：Gear法：令p=k-1，b0=b1=0代入上式得：,多步法的数值稳定性例子（1）,将三阶的Gear法的系数代入得：将从0取到2，三阶的Gear法的绝对稳定域即可绘出，如图中阴影区域。,多步法的数值稳定性例子（2）,三阶的Gear法的绝对稳定域,多步法的数值稳定性例子（3）,令p=k-1，a1=a2=ap，得：同样，代入三阶系数有：所以三阶的Adams-moulton法的绝对稳定域如图中阴影区域所示。,Adams-moulton法的稳定性分析（1）,三阶的Adams-moulton法的绝对稳定域,Adams-moulton法的稳定性分析（2）,令p=k-1，b-10，a1=a2=ap=0，同理有三阶的Adams-Bashforth法的绝对稳定域如图中阴影区域所示。,Adams-Bashforth法的稳定性分析（1）,三阶的Adams-Bashforth法的绝对稳定域,Adam