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

.,1,第九章常微分方程的数值解法,1、引言2、初值问题的数值解法-单步法3、龙格-库塔方法4、收敛性与稳定性5、初值问题的数值解法多步法6、方程组和刚性方程7、习题和总结,主要内容,.,2,主要内容,研究的问题数值解法的意义,1、引言,.,3,现实世界中大多数事物,内部联系非常复杂,找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系,此种关系的数学表达就为,微分方程,1.什么是微分方程?,其状态随着时间、地点、条件的不同而不同,.,4,如何利用数值方法求解微分方程（组）的问题。,2.数值求解微分方程的意义,如何建立数学模型已在建模课程中得到讨论，各类微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，,本章专门讨论,.,5,寻找解析解的过程称为求解微分方程组。,一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解(或古典解),称为真解或解。,3.什么是微分方程(组)的解析解?,3.什么是微分方程(组)的解析解?,.,6,4.什么是微分方程的数值解?,虽然求解微分方程有许多解析方法,但解析方法只能够求解一些特殊类型的方程,从实际意义上来讲,我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值.,寻找数值解的过程称为数值求解微分方程。,把这样一组近似解称为微分方程在该范围内的,数值解,.,7,在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数,求解数值解,很多微分方程根本求不到问题的解析解！,重要手段。,.,8,常微分方程的数值解法常用来求近似解,根据提供的算法,通过计算机,便捷地实现,5.常微分方程数值解法的特点,数值解法得到的近似解(含误差)是一个离散的函数表.,.,9,本章主要讨论一阶常微方程的初值问题,6.基本知识,其中f(x,y)是已知函数，(1.2)是定解条件也称为初值条件。,各种数值解法,.,10,则称f(x,y)对y满足李普希兹条件，L称为Lipschitz常数.,常微分方程的理论指出：,当f(x,y)定义在区域G=(axb,y),若存在正的常数L使：,就可保证方程解的存在唯一性,(Lipschitz)条件,.,11,我们以下的讨论，都在满足上述条件下进行.,若f(x,y)在区域G连续，关于y,一阶常微分方程的初值问题的解存在且唯一.,满足李普希兹条件,一阶常微分方程组常表述为：,方程组,初值条件,.,12,写成向量形式为,高阶常微分方程定解问题如二阶定解问题：,.,13,也就解决了高阶方程的定解问题.,这些解法都可以写成向量形式,用于一阶常微分方程组的初值问题.,.,14,简单的数值方法与基本概念,1.简单欧拉法（Euler）2后退的欧拉法3梯形法4改进Euler法,2、初值问题的数值解法单步法,.,15,1.简单的欧拉(Euler)方法,考虑模型：,在精度要求不高时,通过欧拉方法的讨论,弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.,欧拉方法,最简单而直观实用方法,.,16,2.欧拉方法的导出,把区间a,b,分为n个小区间,步长为,要计算出解函数y(x)在一系列节点,节点,处的近似值,N等分,.,17,对微分方程(1.1)两端从,进行积分,.,18,右端积分用左矩形数值求积公式,得,.,19,亦称为欧拉折线法/*Eulerspolygonalarcmethod*/,每步计算,只用到,或用向前差商近似导数,依上述公式逐次计算可得：,例题,.,20,3.欧拉公式有明显的几何意义,依此类推得到一折线,故也称Euler为单步法。,公式右端只含有已知项,所以又称为显格式的单步法。,.,21,也称欧拉折线法.,就是用这条折线近似地代替曲线,欧拉方法,.,22,从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。,4.欧拉法的局部截断误差：,在假设第i步计算是精确的前提下，考虑截断误差,称为局部截断误差,/*localtruncationerror*/。,.,23,Ri的主项/*leadingterm*/,欧拉法的局部截断误差：,欧拉法具有1阶精度。,.,24,如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.,Taylor展开式,一元函数的Taylor展开式为:,若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。,Ri的主项/*leadingterm*/,.,25,5.欧拉公式的改进：,隐式欧拉法/*implicitEulermethod*/,由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。,一般先用显式计算一个初值，再迭代求解。,隐式欧拉法的局部截断误差：,即隐式欧拉公式具有1阶精度。,.,26,6.梯形公式/*trapezoidformula*/,显、隐式两种算法的平均,注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。,中点欧拉公式/*midpointformula*/,假设，则可以导出即中点公式具有2阶精度。,需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。,.,27,简单,精度低,稳定性最好,精度低,计算量大,精度提高,计算量大,精度提高,显式,多一个初值,可能影响精度,.,28,改进欧拉法/*modifiedEulersmethod*/,注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,.,29,3龙格-库塔法/*Runge-KuttaMethod*/,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,建立高精度的单步递推格式。,.,30,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,Step2:将K2代入第1式，得到,.,31,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有个未知数，个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,注意到，就是改进的欧拉法。,Q:为获得更高的精度，应该如何进一步推广？,.,32,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,.,33,由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。,深入研究龙格-库塔法请看此处！,.,34,4收敛性与稳定性/*ConvergencyandStability*/,1.收敛性/*Convergency*/,例：就初值问题考察欧拉显式格式的收敛性。,解：该问题的精确解为,欧拉公式为,对任意固定的x=xi=ih，有,.,35,2.稳定性/*Stability*/,例：考察初值问题在区间0,0.5上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。,1.00002.00004.00008.00001.60001013.2000101,1.00002.50001016.25001021.56251023.90631039.7656104,1.00002.50006.25001.56261013.90631019.7656101,1.00004.97871022.47881031.23411046.14421063.0590107,Whatiswrong?!,.,36,一般分析时为简单起见，只考虑试验方程/*testequation*/,常数，可以是复数,.,37,例：考察隐式欧拉法,可见绝对稳定区域为：,注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。,.,38,例：隐式龙格-库塔法,而显式14阶方法的绝对稳定区域为,其中2阶方法的绝对稳定区域为,无条件稳定,.,39,例1用欧拉方法，隐式欧拉方法和欧拉中点公式计算,的近似解，取步长h=0.1，并与精确值比较,解:欧拉方法的算式为：,欧拉隐式方法在本题可解出方程，不必迭代，公式为：,欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式,.,40,本题精确解为y=e-x,计算结果见表9-1,例2用欧拉公式和梯形公式的预估校正法计算：,的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:,解:本例中欧拉公式为：,.,41,梯形公式只校正一次的格式为,结果列入表9.2,.,42,1.Runge-Kutta法的一般形式2.2阶Runge-Kutta方法3.经典Runge-Kutta方法4R-K-Fehhlberg方法5.隐式R-K方法6.变步长方法,龙格库塔法深入研究,.,43,1.Runge-Kutta法的一般形式,Runge-Kutta法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为：,.,44,式(9.11)称L级p阶Runge-Kutta方法(简称R-K法)。当L1就是欧拉法，当L2时为改进的欧拉法。,其中,它的局部截断误差是,(9.11),2.2级2阶Runge-Kutta方法,令L=2，则,.,45,其局部截断误差是：,这是3个未知量的两个方程，其中有一个自由参数，方程组有无穷多解，从而得到一族2级2阶R-K方法。,.,46,3.经典Runge-Kutta方法,我们可以构造出一族3级3阶，一族4级4阶和一族5级4阶等R-K方法。最常用的4级4阶是如下的经典R-K方法：,.,47,4R-K-Fehhlberg方法,R-K-Fehhlberg方法是在R-K方法的基础上引进误差和步长控制的办法。即利用5阶R-K法,估计4阶R-K的局部误差，其中,.,48,注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！,.,49,5.隐式R-K方法,类似于显式R-K公式，稍加改变，就得到隐式R-K方法。,它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。,.,50,通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-K法有:,1级2阶中点公式:,2级2阶梯形公式：,2级4阶R-K公式：,.,51,6.变步长方法,在单步法中每一积分步步长实际上是相互独立的，步长的选择具有了灵活性。根据合理地选择每一积分步的步长，既保证精度的要求，又可以减少计算量，从而减小舍入误差。其方便的控制手段是基于误差的事后估计式。,对于给定的精度,如果,反复将步长折半进行计算，直至为止，这时再将步长折半一次，就得到所要的结果。,这种通过加倍或折半处理步长的计算方法称为变步长方法。,注：推荐使用精度好计算量低的变步长方法。用四阶显式R-K方法做变步长方法是实践中较好的方法！,.,52,例分别用改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：,.,53,节点改进欧拉法四阶龙格库塔法准确解01110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051,表,（注：已指出过准确解,）,.,54,单步法的相容性、收敛性和稳定性,1.单步法的相容性2.单步法的收敛性3.单步法的稳定性4.相容性和收敛性的关系5.相容性和方法阶的关系6.稳定性和收敛性7.绝对稳定性和绝对稳定域,.,55,单步法的相容性,定义一：对于（9.1.1）常微分方程初值问题单步法的形式可以变表示为（9.2.19）其中h为步长若对求解区间中任一固定的x，当时皆成立,则称由（9.2.19）确定的单步法与微分方程初值问题是相容的注意到上式左边极限为由（9.1.1）知它应等于从而由相容性定义得称相容性条件。,.,56,单步法的收敛性,定义二：设y(x)是（9.1.1）的解，是单步法（9.2.19）产生的数值解，对于每一个固定的,，当即。若成立，则称该方法是收敛的。,.,57,单步法的稳定性,定义三：若一个数值方法在基点处的值有的扰动，在此后各基点（mn）处的值产生的偏差均不超过，则称该方法是稳定的。单步法的稳定性有以下定理,定理二：若单步法的增量函数对变量y满足Lipschtiz条件，则单步方法是稳定的。,.,58,相容性和收敛性的关系,定理一：若单步法的增量函数对变量y满足Lipschtiz条件，即存在与h,x无关的常数L，对区域D=任意两点（x，y1）,(x,y2)成立，则单步法收敛的充分必要条件是相容性条件成立。（读者自证）,.,59,相容性和方法阶的关系,若单步法是p阶方法则成立若单步法满足相容性条件，得所以=0也就是说单步法的阶数一定要是正数。由于我们考虑的单步法皆为正整数，p至少为一。因此我们考虑的单步法都满足相容性条件。,.,60,稳定性和收敛性的关系,若单步法的增量函数满足定理二的条件即单步法是稳定的则单步法收敛的充分必要条件是相容性条件成立。,.,61,绝对稳定性和绝对稳定域,稳定性问题是一个比较复杂的问题。为了简化讨论一般仅对试验方程进行考察。这里假定Re0和的允许范围，称为该方法的绝对稳定域。,.,62,绝对稳定性和绝对稳定域2,将Euler方法应用到试验方程得误差方程是要求误差不增长则,.,63,则Euler方法的绝对稳定域是以为半径的圆的内部。为了保证稳定性步长有所控制。假如当时h应满足，当时，h应满足等等。同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。,.,64,1.Adams方法2.米尔尼方法、汉明方法及辛普森方法3.预测校正方法4.多步法的相容性、稳定性和收敛性,5初值问题的数值解法多步法,.,65,考虑型如的k步法，称为阿当姆斯（Adams）方法,1.Adams方法,拿其中一种为例，推导其公式。对上面所说的法一般形式若取，有方程组同样解之，得到3步4阶隐式Adams公式和它的余项。,.,66,其它请读者自证。我们仅将结论列于下表。Adams显式公式,.,67,Adams隐式公式,注：在这些方法中，4阶的Adams预测校正方法具有方法简洁、使用方便、易排序、高精度等优点。尤其当函数f比较复杂时更能显出它的优越性。,.,68,同理得到5个待定参数方程组。解之得，。构成著名的Miline4步4阶显式公式和它的余项。,，,2.米尔尼方法、汉明方法及辛普森方法,.,69,同理得到5个参数方程组。求解后就构成著名的3步4阶隐式Hamming公式和它的余项。若取，也得到5个参数方程组。求解后就构成Simpson隐式公式和其余项。,米尔尼方法、汉明方法及辛普森方法,.,70,不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。,3.预测校正方法,.,71,对于线性多步公式，一般地，不难验证，如果预测公式是阶或阶精度，校正公式为阶精度，则用预测公式提供初值，校正公式迭代一次的效果也能达到阶精度，再迭代下去，效果就不明显了。预测-校正技术即保证了计算精度，又使隐式计算显式化，克服了隐式公式要反复迭代的困难。至于实际计算步长的选取，我们对预测-校正公式使用外推原理，得到误差估计式，用来调整计算步长，使达到控制误差要求。,对于线性多步方法常用的预测-校正公式有Miline-Hamming方法和4阶Adams显隐式预测-校正公式,.,72,将Miline公式和Hamming公式结合，构成预测-校正公式Miline公式和Hamming公式的局部截断误差分别为,修正Miline-Hamming公式,.,73,利用外推原理，即加上后消去局部截断误差的主项。使说明经过外推后的算法其精度提高一阶。忽略误差项，上式可改写为,.,74,由于和是在计算过程中获得的数据，称为Miline公式和Hamming公式的事后误差估计式。我们可以用它们来调节计算步长的大小，即选择一个合适的的步长，使，其中的是要求达到的计算精度。,.,75,又可得到Miline公式和Hamming公式的修正公式，它们分别是从而构成如下的修正Hamming预测-校正公式，简称修正Hamming公式：,.,76,预测：修正：校正：修正：在应用这套公式时，先由同阶单步法提供初值，。计算时可取。,.,77,将4步4阶显式Adams公式作为预测公式和3步4阶式Adams公式作为校正公式，构成Adams预测-校正公式。它们的局部截断误差分别是,Adams预测-校正公式,.,78,利用外推原理，将上两式作线性组合消去局部截断误差的主项。使计算精度至少提高一阶，同时得到两个修正公式，将它们和上两式构成了如下修正Adams公式：预测：修正：校正:,.,79,修正：同理，在计算时，调节计算步长使,由同阶单步法提供初值,，。计算时，可取。理论分析得出它们的绝对稳定区域，修正Hamming法：;4阶Adams预测-校正法：其中,我们也可以在教学演示上看出修正的预测校正格式的误差非常的小。,.,80,多步法的相容性条件k步法的一般形式为其中由对单步法的讨论可知,当时,方法阶数至少为1。即对y1,y=x应精确成立。令y1，所以令y=x得,4.多步法的相容性、稳定性和收敛性,.,81,所以我们称为线性多步法的相容性条件。,.,82,k步法需要k个出发值，而初值问题只提供了一个初值，在使用k步法时尚需要其它方法作补充k-1个出发值，今假定它们是，当这k-1个值都应收敛于共同的极限y（a）即在讨论多步法收敛性时我们总假定（9.3.12）成立定义四：,多步法的收敛问题,.,83,的解收敛于初值问题的解y(x)。这里xa+nh.定义五：称多项式为k步法的特征多项式。如果特征多项式的零点皆不大于1，且等于1的零点是单重的，称根条件成立。称多项式为第二特征多项式。显然根条件可以表示定理三：k步法收敛的充要条件为：（1）相容性条件成立。（2）特征多项式的零点皆不大于1，且等于1的零点是单重的。（称2为）特征根条件。,.,84,多步法的稳定性,关于线性多步法成立以下定理：定理四：若函数f(x,y)对变量y满足Lipschtiz条件在与h,x无关的常数L，对区域D=任意两点（x，y1）,(x,y2)成立k步法的相容性、收敛性、稳定性有以下关系对于常微分方程右端函数f(x,y)，在相容性条件成立情况下，k步法的收敛性和稳定性是等价的。事实上在相容性条件成立时，收敛性和稳定性的充要条件都是特征根条件成立。,.,85,多步法的绝对稳定性和绝对稳定域,将线性多步法的公式应用到试验方程进行考察。这里假定Re0,即试验方程本身是稳定的。记得是常系数差分方程，其特征方程为记它的k个特征根为并设它们是互异的。显然根与有关，不妨记为注意到当时这时由特征方程得由线性多步法的相容性条件得是一个根。不妨设,差分方程的解为其中系数由线性多步法的出发值确定。,.,86,另一方面，y(0)=1的试验方程的精确解为，设多步法截断误差为，由此可得，我们称为主根，其它根都为增根。定义五：线性多步法的绝对稳定区域对给定的，如果特征方程的特征根皆按模小于1，则线性多步法关于u是绝对稳定的。使得成立的构成绝对稳定区域。注：从误差角度来看绝对稳定区域的方法是一个理想的方法。这样，绝对稳定区域越大，方法适用性越广，因而越优越。,.,87,6方程组和刚性方程,9.1一阶方程组9.2化高阶方程为一阶方程组9.3刚性方程组,.,88,9.1一阶方程组,前面我们研究了单个方程y=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形。考虑一阶方程组的出值问题若采用向量的记号，记,.,89,则常微分方程组的出值问题可以表示为前几节我们主要讨论了常微分方程的出值问题的数值解法。只要将y和f改写为向量，那么前面提供的各种计算公式即可适用于一阶常微分方程组的出值问题。Runge-Kutta方法对于方程组四级四阶显示Runge-Kutta公式,.,90,.,91,若写成分量形式就是i=1,2,m为了帮助理解这一公式的计算过程，我们再考虑两个方程的特殊情形：,.,92,这时四阶龙格库塔公式具有形式,.,93,其中,.,94,这是一步法，利用节点上的值，由上式顺序计算，然后代入即可求得节点上的。,.,95,9.2化高阶方程为一阶方程组,关于高阶微分方程的初值问题，原则上总可以归结为一阶方程组来求解，例如，考察下列m阶微分方程初始条件为只要引进新的变量即可将m阶方程化为如下的一阶方程组,.,96,初始条件则相应的化为,.,97,9.3刚性方程组,在求解方程组时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性问题。若线性系统，其中，中A的特征值满足条件0(j=1,N),.,98,且则称系统为刚性方程，称s为刚性比,.,99,7习题与总结,1数值例题2数值练习3总结,.,100,1、数值例题,我们已经学习了很多数值算法，他们的效果到底如何呢？下面我们来分析一道例题，看看那些方法，就这个问题，最能接近真实值求初值问题的数值解，取h=0.1并同精确解比较（1）用欧拉法来计算这个初值问题（2）用各阶的RungeKutta方法来计算这个初值问题（3）用四阶的Adams定步长，Adams定步长预测校正方法来计算这个初值问题。然后将数值结果列成表格：,.,101,.,1