常微分方程数值解法课件_002

第7章 常微分方程数值解法7.1 引言引言7.2 简单的数值方法与基本概念（欧拉法）简单的数值方法与基本概念（欧拉法）7.3-库塔方法库塔方法7.5 方程组和高阶方程方程组和高阶方程7.1 引 言 科学技术中常常需要求解常微分方程的定解问题科学技术中常常需要求解常微分方程的定解问题.这类问题最简单的形式，是本章将着重考察的这类问题最简单的形式，是本章将着重考察的一阶方一阶方程的初值问题程的初值问题 我们知道，只有我们知道，只有f(x,y)适当光滑适当光滑譬如关于譬如关于y满足满足莱布尼兹莱布尼兹(Lipschitz)条件条件理论上就可以保证初值问题的解理论上就可以保证初值问题的解yf(x)存在并且唯一存在并且唯一.虽然求解常微分方程有各种各样的解析方法，但虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法题中归结出来的微分方程主要靠数值解法.所谓所谓数值解法数值解法,就是寻求解就是寻求解y(x)在一系列离散节点在一系列离散节点上的近似值上的近似值 y1,y2,yn,yn+1,.相邻两个节点的间距相邻两个节点的间距hn=xn+1-xn称为称为步长步长.今后如不特别说明，总是假定今后如不特别说明，总是假定 hi=h(i=1,2,)为为常数常数,这时节点为这时节点为xn=x0+nh(i=0,1,2,)(等距节点等距节点).截去最后一项截去最后一项 ,可得可得这就是著名的这就是著名的(显式显式)欧拉欧拉(Euler)公式公式.若初值若初值y0已已知，则依公式知，则依公式(2.1)可逐次逐步算出各点数值解可逐次逐步算出各点数值解.即即7.2 简单的数值方法与基本概念7.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法用向前差商代替导数用向前差商代替导数 例例1 用欧拉公式求解初值问题用欧拉公式求解初值问题 解解 取步长取步长h=0.1，欧拉公式的具体形式为，欧拉公式的具体形式为其中其中xn=nh=0.1n(n=0,1,10),已知已知y0=1,由此式可得由此式可得依次计算下去，依次计算下去，部分计算结果部分计算结果见下表见下表.与准确解与准确解 相比，可看出欧拉公式的计算相比，可看出欧拉公式的计算结果精度很差结果精度很差.xn 欧拉公式数值解欧拉公式数值解yn准确解准确解y(xn)误差误差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 欧拉公式具有明显的几何意义欧拉公式具有明显的几何意义,就是就是用折线近似用折线近似代替方程的解曲线代替方程的解曲线，因而常称公式，因而常称公式(2.1)为为欧拉折线法欧拉折线法.还可以通过几何直观来考察欧拉方法的精度还可以通过几何直观来考察欧拉方法的精度.假假设设yn=y(xn),即顶点即顶点Pn落在积分曲线落在积分曲线y=y(x)上，那么，上，那么，按欧拉方法做出的折线按欧拉方法做出的折线PnPn+1便是便是y=y(x)过点过点Pn的切线的切线.从图形上看从图形上看,这这样定出的顶点样定出的顶点Pn+1显著显著地偏离了原来的积分曲地偏离了原来的积分曲线，可见欧拉方法是线，可见欧拉方法是相相当粗糙当粗糙的的.为了分析计算公式的精度，通常可用泰勒展开为了分析计算公式的精度，通常可用泰勒展开将将y(xn+1)在在xn处展开，则有处展开，则有在在yn=y(xn)的前提下，的前提下，f(xn,yn)=f(xn,y(xn)=y(xn n).于于是可得欧拉法是可得欧拉法(2.1)的的公式误差公式误差为为称为此方法的称为此方法的局部截断误差局部截断误差.截去最后一项截去最后一项 ,可得可得这就是著名的这就是著名的(隐式隐式)后退欧拉后退欧拉(Euler)公式公式.若初值若初值y0已知，则依公式已知，则依公式(2.5)可逐次逐步算出各点数值解可逐次逐步算出各点数值解.即即用向后差商代替导数用向后差商代替导数 后退的欧拉公式与欧拉公式有着本质的区别后退的欧拉公式与欧拉公式有着本质的区别,后后者是关于者是关于yn+1的一个直接计算公式，这类公式称作是的一个直接计算公式，这类公式称作是显式的显式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它实际上，它实际上是关于是关于yn+1的一个函数方程的一个函数方程,这类方程称作是这类方程称作是隐式的隐式的.显式显式与与隐式隐式两类方法各有特点，考了到数值稳两类方法各有特点，考了到数值稳定性等其他因素，人们有时需要选用定性等其他因素，人们有时需要选用隐式隐式方法，但方法，但使用使用显式显式算法远比算法远比隐式隐式方便方便.隐式方程通常用迭代法求解，而迭代过程的实隐式方程通常用迭代法求解，而迭代过程的实质是质是逐步逐步显式化显式化.7.2.2 梯形方法 为得到比为得到比欧拉法精度高的计算公式，将欧拉公式欧拉法精度高的计算公式，将欧拉公式和后退的欧拉公式进行算数平均，可得和后退的欧拉公式进行算数平均，可得称为称为梯形公式梯形公式.梯形公式梯形公式是是隐式单步法隐式单步法，用迭代法求解，同后，用迭代法求解，同后退的欧拉方法一样，仍用欧拉法提供迭代初值。退的欧拉方法一样，仍用欧拉法提供迭代初值。7.2.3 单步法的局部截断误差与阶 初值问题初值问题(1.1),(1.2)的的单步法单步法可用可用一般形式一般形式表示为表示为其中多元函数其中多元函数 与与f(x,y)有关，当有关，当 含有含有yn+1时，方法时，方法是是隐式隐式的，若不含的，若不含yn+1则为则为显式方法显式方法，所以显式单步，所以显式单步法可表示为法可表示为(x,y,h)称为称为增量函数增量函数，例如对欧拉法，例如对欧拉法(2.1)有有它的局部截断误差它的局部截断误差(2.3)已由给出已由给出,对一般显式单步法对一般显式单步法则可如下定义则可如下定义.定义定义1 设设y(x)是初值问题是初值问题(1.1),(1.2)的的准确解准确解,称称为显式单步法为显式单步法(2.10)的的局部截断误差局部截断误差.Tn+1之所以称为之所以称为局部的局部的，是假设在，是假设在xn前各步没有前各步没有误差误差.当当yn=y(xn)时，计算一步，则有时，计算一步，则有 所以，局部截断误差可理解为用方法所以，局部截断误差可理解为用方法(2.10)计算计算一步的误差，也即公式一步的误差，也即公式(2.10)中用准确解中用准确解y(x)代替数代替数值解产生的公式误差值解产生的公式误差.根据定义根据定义,显然显然欧拉法的局部欧拉法的局部截断误差截断误差为为显然显然Tn+1=O(h2).一般情形的定义如下一般情形的定义如下 定义定义2 设设y(x)是初值问题的准确解，若存在最大是初值问题的准确解，若存在最大整数整数p使显式单步法使显式单步法(2.10)的的局部截断误差局部截断误差满足满足则称方法则称方法(2.10)具有具有p阶精度阶精度.若将若将(2.10)展开式写成展开式写成则则 称为称为局部截断误差主项局部截断误差主项.以上定义对隐式单步法以上定义对隐式单步法(2.9)也是适用的也是适用的.例如，例如，对后退欧拉法对后退欧拉法(2.5)其局部截断误差为其局部截断误差为这里这里p=1是是1阶方法阶方法，局部截断误差主项局部截断误差主项为为同样对梯形法同样对梯形法(2.7)有有所以梯形方法所以梯形方法(2.7)是是二阶的二阶的.其其局部截断误差主项局部截断误差主项为为7.2.4 改进的欧拉公式我们看到，梯形方法虽然提高了精度，但其算法我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式复杂，在应用迭代公式(2.9)进行实际计算时，每迭代进行实际计算时，每迭代一次，都要重新计算函数一次，都要重新计算函数f(x,y)的值，而迭代又要反的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测复进行若干次，计算量很大，而且往往难以预测.为为了控制计算量，通常只迭代一两次就转入下一步的计了控制计算量，通常只迭代一两次就转入下一步的计算，这算，这就简化了算法就简化了算法.具体地说，我们先用欧拉公式求得一个初步的近具体地说，我们先用欧拉公式求得一个初步的近似值似值 ，称之为，称之为预测值预测值，此预测值，此预测值 的精度可能很的精度可能很差，再用梯形公式差，再用梯形公式(2.7)将它校正一次，即按将它校正一次，即按(2.8)式迭式迭代一次，这个结果代一次，这个结果称之为称之为校正值校正值.这样建立的这样建立的预测预测校正系统校正系统通常称为通常称为改进的欧拉公式改进的欧拉公式:或表示为下列或表示为下列平均化形式平均化形式(2.13)预测预测校正校正 例例2 用改进的欧拉法用改进的欧拉法解解例例1中的初值问题中的初值问题(2.2).解解 仍取步长仍取步长h=0.1,改进的欧拉公式为改进的欧拉公式为部分计算结果见下表部分计算结果见下表 xnyn 误差误差 xnyn 误差误差00.20.4 1.1.184096 1.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.005818同例同例1中的欧拉法的计算结果比较，明显改善了精度中的欧拉法的计算结果比较，明显改善了精度.xn 欧拉公式欧拉公式 改进欧拉公式改进欧拉公式 yn 误差误差 yn 误差误差 0.0 0.2 0.4 0.6 0.8 1.0 1.1.191818 1.358213 1.508966 1.649783 1.784770 00.0086020.016572 0.0257260.0373310.052719 1 1.184096 1.343360 1.485956 1.616476 1.737869 00.0000880.0017190.0027160.0040240.005818 7.3 龙格库塔方法 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方析这两个公式的特点，从而探索一条构造高精度方法的途径法的途径.上节给出了显式单步法的表达式上节给出了显式单步法的表达式(2.10),其局其局部截断误差为部截断误差为O(hp+1)，对欧拉法，对欧拉法Tn+1=O(h2)，即方法，即方法为为p=1阶阶，若用改进欧拉法，若用改进欧拉法(2.13)，它可表为，它可表为此时增量函数为此时增量函数为一般说来，点数一般说来，点数r 越多，精度越高，上式右端相当于越多，精度越高，上式右端相当于增量函数增量函数(xn,yn,h)，为得到便于计算的显式方法，为得到便于计算的显式方法，可类似于改进欧拉法可类似于改进欧拉法(3.1),(3.2)，将公式表示为将公式表示为其中其中这里这里ci,ai,bij均为常数均为常数.(3.4)和和(3.5)称为称为r级显式级显式龙格龙格-库塔库塔(Runge-Kutta)法法,简称简称R-K方法方法.7.3.2 二阶显式R-K方法 对对r=2的的R-K方法方法，由，由(3.4),(3.5)式可得式可得如下如下计算计算公式公式这里这里 c1,c2,a2,b21 均为待定常数，我们希望适当选取均为待定常数，我们希望适当选取这些系数，使公式阶数这些系数，使公式阶数 p 尽量高尽量高.则由此可以看出在改进的欧拉公式中相当于取则由此可以看出在改进的欧拉公式中相当于取(xn,yn),(xn+1,yn+1)两两点处点处斜率的平均值斜率的平均值，近似代替平，近似代替平均斜率，其精度比欧拉公式提高了均斜率，其精度比欧拉公式提高了.c1=c2=1/2,a2=b21=a=1=1时，即为改进的欧拉公式：时，即为改进的欧拉公式：称为称为变形的欧拉公式变形的欧拉公式,也可以表示为也可以表示为:如取如取a=1，则，则c1=0,c2=1,a2=b21=1/2.得计算公式得计算公式 对对r=2的的R-K公式公式能否使局部误差提高到能否使局部误差提高到O(h4)?为此为此 需把需把k2多展开一项，多展开一项，展开式中的项展开式中的项 是不能通过选择参数削掉的，实际上要使是不能通过选择参数削掉的，实际上要使 h3 的项为的项为零，需增加零，需增加3个方程，要确定个方程，要确定4个参数个参数c1,c2,2,21，这，这是不可能的是不可能的.故故r2的显式的显式R-K方法的阶只能是方法的阶只能是p=2，而不能得到三阶公式而不能得到三阶公式.四阶四阶R-K方法方法的每一步需要计算四次函数值的每一步需要计算四次函数值f，可以证明其局部截断误差为可以证明其局部截断误差为O(h5).7.3.3 标准四阶显式R-K方法 然而值得指出的是，龙格然而值得指出的是，龙格-库塔方法的推导基于库塔方法的推导基于泰勒展开泰勒展开方法，因而它要求所求的方法，因而它要求所求的解具有较好的光解具有较好的光滑性质滑性质.反之，如果解的光滑性差，那么，使用龙反之，如果解的光滑性差，那么，使用龙格格-库塔方法求得的数值解，其精度可能反而不如改库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法进的欧拉方法.实际计算时，我们应当实际计算时，我们应当针对问题的具针对问题的具体特点选择合适的算法体特点选择合适的算法.7.3.4 变步长的龙格-库塔方法 单从每一步看，步长越小，截断误差就越小，单从每一步看，步长越小，截断误差就越小，但随着步长的缩小，在一定求解范围内所要完成的步但随着步长的缩小，在一定求解范围内所要完成的步数就增加了数就增加了.步数的增加不但引起计算量的增大，而步数的增加不但引起计算量的增大，而且可能导致舍入误差的严重积累且可能导致舍入误差的严重积累.因此同积分的数值因此同积分的数值计算一样，微分方程的数值解法也有个计算一样，微分方程的数值解法也有个选择步长的问选择步长的问题题.在选择步长时，需要考虑在选择步长时，需要考虑两个问题两个问题：1.怎样衡量和检验计算结果的精度？怎样衡量和检验计算结果的精度？2.如何依据所获得的精度处理步长？如何依据所获得的精度处理步长？我们考察我们考察四阶四阶R-K公式公式(3.13)，从节点，从节点xn出发，出发，先以先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于公式，由于公式的局部截断误差为的局部截断误差为O(h5)，故，故然后将步长折半，即取为步长然后将步长折半，即取为步长 ,从从xn跨两步到跨两步到xn+1，再求得一个近似值再求得一个近似值 ，每跨一步的局部截断误差，每跨一步的局部截断误差是是 ，因此有，因此有比较比较(3.14)式和式和(3.15)式我们看到，步长折半后，误差式我们看到，步长折半后，误差大约减少到大约减少到1/16，即有，即有由此易得下列事后估计式由此易得下列事后估计式这样，我们可以通过检查步长，折半前后两次计算结这样，我们可以通过检查步长，折半前后两次计算结果的偏差果的偏差来判定所选的步长是否合适，具体地说，将区分以下来判定所选的步长是否合适，具体地说，将区分以下两种情况处理：两种情况处理：这种通过加倍或折半处理步长的方法称为变步长这种通过加倍或折半处理步长的方法称为变步长方法方法.表面上看，为了选择步长，每一步的计算量增表面上看，为了选择步长，每一步的计算量增加了，但总体考虑往往是合算的加了，但总体考虑往往是合算的.1.1.对于给定的精度对于给定的精度,如果如果,我们反复将步长我们反复将步长折半计算折半计算,直至直至为为止止,这时取最终得到的这时取最终得到的 作为作为结果；结果；2.2.如果如果为止，这时再将步长折半计算一次，就得到所为止，这时再将步长折半计算一次，就得到所要的结果要的结果.7.5 方程组和高阶方程7.5.1 一阶方程组一阶方程组 前面我们研究了单个方程前面我们研究了单个方程y=f 的数值解，只的数值解，只要把要把y 和和f 理解为向量，那么，所提供的各种计算公理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形式即可应用到一阶方程组的情形.考察考察一阶方程组一阶方程组的初值问题，初始条件给为的初值问题，初始条件给为若采用向量的记号，记若采用向量的记号，记(向量向量)则上述方程组的初值问题可表示为则上述方程组的初值问题可表示为求解这一初值问题的求解这一初值问题的四阶龙格四阶龙格-库塔公式库塔公式为为(向量向量)式中式中(向量向量)或表示为或表示为(分量分量)其中其中(分量分量)这里这里yin是第是第i个因变量个因变量yi(x)在节点在节点xn=x0+nh的近似值的近似值.为了帮助理解这一公式的计算过程，我们再考察为了帮助理解这一公式的计算过程，我们再考察两个方程的特殊情形两个方程的特殊情形这时这时四阶龙格四阶龙格-库塔公式库塔公式具有形式具有形式其中其中这是一步法，利用节点这是一步法，利用节点xn上的值上的值yn,zn，由，由(6.3)式式顺序计算顺序计算K1,L1,K2,L2,K3,L3,K4,L4，然后代入，然后代入(6.2)式即可求得节点式即可求得节点xn+1上的上的yn+1,zn+1.关于高阶微分方程关于高阶微分方程(或方程组或方程组)的初值问题，原则的初值问题，原则上总可以归结为一阶方程组来求解上总可以归结为一阶方程组来求解.例如，考察下列例如，考察下列m阶微分方程阶微分方程初始条件为初始条件为只要引进新的变量只要引进新的变量即可将即可将m阶方程阶方程(6.4)化为如下的一阶方程组：化为如下的一阶方程组：7.5.2 化高阶方程为一阶方程组化高阶方程为一阶方程组初始条件初始条件(6.5)则相应地化为则相应地化为不难证明初始条件不难证明初始条件(6.4),(6.5)和和(6.6),(6.7)是彼此等是彼此等价的价的.特别地，对于下列二阶方程的问题特别地，对于下列二阶方程的问题引进新变量引进新变量z=y,即可化为下列一阶方程组的初值问题即可化为下列一阶方程组的初值问题:针对这个问题应用针对这个问题应用四阶龙格四阶龙格-库塔公式库塔公式(6.2)，有，有由由(6.3)式可得式可得写在最后写在最后成功的基成功的基础在于好的学在于好的学习习惯The foundation of success lies in good habits45 结束语当你尽了自己的最大努力时，失败也是伟大的，所以不要放弃，坚持就是正确的。When You Do Your Best,Failure Is Great,So DonT Give Up,Stick To The End演讲人：XXXXXX 时 间：XX年XX月XX日