《九章常微分方程的数值解法》由会员分享，可在线阅读，更多相关《九章常微分方程的数值解法（113页珍藏版）》请在金锄头文库上搜索。
1、 第九章第九章 常微分方程的数值解法常微分方程的数值解法 1、引言引言 2、初值问题的数值解法初值问题的数值解法-单步法单步法 3、龙格龙格-库塔方法库塔方法 4、收敛性与稳定性收敛性与稳定性 5、初值问题的数值解法初值问题的数值解法多步法多步法 6、方程组和刚性方程方程组和刚性方程 7、习题和习题和总结总结主主 要要 内内 容容主主 要要 内内 容容研究的问题研究的问题数值解法的意义数值解法的意义1、 引引 言言现实世界中大多数事物内部联系非常复杂找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系此种关系的数学表达就为微分方程微分方程微分方程微分方程1.1.什么
2、是微分方程什么是微分方程 ? ?其状态随着其状态随着时间、地点、条件时间、地点、条件的不同而不同的不同而不同如何利用数值方法求解微分方程（组）的问题。如何利用数值方法求解微分方程（组）的问题。2.2.数值求解微分方程的意义数值求解微分方程的意义如何建立数学模型已在建模课程中得到讨论，各类微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，本章专门本章专门讨论讨论寻找解析解的过程称为求解微分方程组。寻找解析解的过程称为求解微分方程组。 一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解解析解( (或古典解或古典解),),
3、称为真解或解。称为真解或解。3.什么是微分方什么是微分方程程(组组)的解析解的解析解?3.什么是微分方程什么是微分方程(组组)的解析解的解析解?4.4.什么是微分方程的数值解什么是微分方程的数值解? ?虽然求解微分方程有许多解析方法,但解析方法只能够求解一些特殊类型的方程,从实际意义上来讲我们更关心的是某些 特定的自变量在某一个定义范围内的一系列离散点一系列离散点上的近似值.寻找数值解的过程称为数值求解微分方程。寻找数值解的过程称为数值求解微分方程。把这样一组近似解称为把这样一组近似解称为微分方程在该范围内的微分方程在该范围内的数值解数值解在大量的实际方程中出现的函数起码的连续性都无法保证，更
4、何况要求阶的导数求解数值解很多微分方程很多微分方程根本求不到根本求不到问题的解析解！问题的解析解！重要手段。常微分方程的数值解法常用来求近似解根据提供的算法通过计算机通过计算机便捷地实现便捷地实现5.常微分方程数值解法的特点常微分方程数值解法的特点数值解法得到的近似数值解法得到的近似解解( (含误差含误差) )是一个是一个离散的函数表离散的函数表. .本章主要讨论一阶常微方程的初值问题6.基本知识基本知识其中f (x,y)是已知函数，(1.2)是定解条件也称为初值条件。初值条件。各各种种数数值值解解法法则称 f (x,y) 对y 满足李普希兹条件，L 称为Lipschitz常数.常微分方程的理
5、论指出：常微分方程的理论指出：当 f (x,y) 定义在区域 G=(axb,y)若存在正的常数 L 使：就可保证方程解的存在唯一性就可保证方程解的存在唯一性(Lipschitz)条件条件我们以下的讨论，都在满足上述条件下进行.若 f (x,y) 在区域 G连续，关于y一阶常微分方程的初值问题的解存在且唯一问题的解存在且唯一.满足李普希兹满足李普希兹满足李普希兹满足李普希兹条件条件条件条件一阶常微分方程组常表述为：方程组方程组方程组方程组初值条件初值条件初值条件初值条件写成向量形式为高阶常微分方程定解问题如二阶定解问题：也就解决了高阶方程的定解问题也就解决了高阶方程的定解问题. .这些解法都可以
6、写成向量形式用于一阶常微分方程组的初值问题用于一阶常微分方程组的初值问题. .简单的数值方法与基本概念简单的数值方法与基本概念 1. 简单简单欧拉法（欧拉法（Euler） 2后退的欧拉法后退的欧拉法 3梯形法梯形法 4改进改进EulerEuler法法 2、初值问题的数值解法、初值问题的数值解法单步法单步法1. 简单的欧拉简单的欧拉(Euler)方法方法考虑模型：在精度要求不高时通过欧拉方法的讨论弄清常微方程初值弄清常微方程初值弄清常微方程初值弄清常微方程初值问题数值解法的一问题数值解法的一问题数值解法的一问题数值解法的一些基本概念和构造些基本概念和构造些基本概念和构造些基本概念和构造方法的思路
7、方法的思路方法的思路方法的思路. .欧欧拉拉方方法法最简单而直观最简单而直观最简单而直观最简单而直观实用方法实用方法实用方法实用方法2. 欧拉方法的导出欧拉方法的导出把区间a,b分为n个小区间步长为要计算出解函数 y(x) 在一系列节点节点处的近似值N N等分等分对微分方程(1.1)两端从进行积分右端积分用右端积分用左矩形数值左矩形数值求积公式求积公式得x0x1亦亦称为称为欧拉折线法欧拉折线法 /* Eulers polygonal arc method*/ 每步计算只用到或用向前差商近似导数依上述公式逐次计算可得：例题3.欧拉公式有明显的几何意义依此类推得到一折线故也称Euler为单步法。公
8、式右端只含有已知项所以又称为显格式的单步法。也称欧拉折线法.就是用这条折线近似地代替曲线欧拉方法欧拉方法从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。非常粗糙。4.欧拉法的局部截断误差：在假设第 i 步计算是精确的前提下，考虑截断误差称为局部截断误差/* local truncation error */。若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。Ri 的的主项主项/* leading term */ 欧拉法的局部截断误差：欧拉法具有 1 阶精度。在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri =
9、 y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p p阶方法阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.Taylor展开式,一元函数的Taylor展开式为:若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该算法有则称该算法有p 阶精度。阶精度。Ri 的的主项主项/* leading term */5. 欧拉公式的改进欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数x0x
10、1)(,()(1101xyxfhyxy+ + )1,., 0(),(111 = =+ += =+ + + +niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故同时出现在等式的两边，不能直接得到，故称为称为隐式隐式 /* implicit */ 欧拉公式，而前者称为欧拉公式，而前者称为显式显式 /* explicit */ 欧拉公式。欧拉公式。一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即即隐式欧拉公式具有隐式欧拉公式具有 1 阶精度。阶精度。6.梯形公式梯形公式 /
11、* trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。阶精度，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式隐式公式，计算时不得不用到公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。 中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数x0x2x1假设假设 ，则可以导出，则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。需要需要2个初值
12、个初值 y0和和 y1来启动递推来启动递推过程，这样的算法称为过程，这样的算法称为双步法双步法 /* double-step method */，而前面的三种算法都是而前面的三种算法都是单步法单步法 /* single-step method */。方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多多一个初值一个初值, 可能影响精度可能影响精度 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1:
13、 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+ += =+ +Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+ +iy),(),(2111+ + + + + += =iiiiiiyxfyxfhyy注：注：此法亦称为此法亦称为预测预测-校正法校正法 /* predictor-corrector method */。可以证明该算法具有可以证明该算法具有 2 阶精度，同时可以看到它是个阶精度，同时可以看到它是个单单步步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将。后面将
14、看到，它的看到，它的稳定性高稳定性高于显式欧拉法。于显式欧拉法。3 龙格龙格 - 库塔法库塔法 /* Runge-Kutta Method */ 考察改进的欧拉法，可以将其改写为：考察改进的欧拉法，可以将其改写为：斜率斜率一定取一定取K1 K2 的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗吗？单步递推法的单步递推法的基本思想基本思想是从是从 ( xi , yi ) 点出发，以点出发，以某一斜某一斜率率沿直线达到沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所点。欧拉法及其各种变形所能达到的最高精度为能达到的最高精度为2阶阶。建立高精度的单步递推格式。建立高精度
15、的单步递推格式。首先希望能确定系数首先希望能确定系数 1、 2、p，使得到的算法格式有使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1: 将将 K2 在在 ( xi , yi ) 点作点作 Taylor 展开展开将将改进欧拉法推广为：改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+ + += = =+ + += =+ + Step 2: 将将 K2 代入第代入第1式，得到式，得到Step 3: 将将 yi+1 与与 y( xi+1 ) 在在 xi 点点的的泰勒泰勒展开作比较展开作比较要求要求 ，则必
16、须有：，则必须有：这里有这里有 个未知个未知数，数， 个方程。个方程。32存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙格阶龙格 - 库库塔格式塔格式。注意到，注意到， 就是改进的欧拉法。就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？其中其中 i ( i = 1, , m )， i ( i = 2, , m ) 和和 ij ( i = 2, , m; j = 1, , i 1 ) 均为待均为待定系数，确定这些系数定系数，确定这些系数的步骤与前面相似。的步骤与前面相似。 ).,(.),(),(),(.
17、1122112321313312122122111 + + + + + + += =+ + + += =+ + += = =+ + + + += =mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法 /* Classical Runge-Kutta Method */ ：注：注： 龙格龙格-库塔法库塔法的主要运算在于计算的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高
18、精度阶数的关系：度阶数的关系：753可可达到的最高精度达到的最高精度642每步须算每步须算Ki 的的个数个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。深入研究龙格深入研究龙格-库塔法请看库塔法请看此处此处！4 收敛性与稳定性收敛性与稳定性 /* Convergency and Stability */ 1.收敛性收敛性 /* Convergency */ 若若某某算算法法对对于于任任意意固固定定的
19、的 x = xi = x0 + i h，当当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是则称该算法是收敛收敛的。的。 例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉欧拉公式为公式为对对任意固定的任意固定的 x = xi = i h ，有有 2.稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30
20、.40.5精确解精确解改进欧拉法改进欧拉法 欧欧拉拉隐式隐式欧拉欧拉显式显式 节点节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?!若若某某算算法法在在计计算算过过
21、程程中中任任一一步步产产生生的的误误差差在在以以后后的的计计算算中中都都逐逐步步衰衰减减，则则称称该该算算法法是是绝绝对对稳稳定定的的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */常数，可以常数，可以是复数是复数当当步长取为步长取为 h 时，将某算法应用于上式，并假设只在初值时，将某算法应用于上式，并假设只在初值产生误差产生误差 ，则若此误差以后逐步衰减，就称该，则若此误差以后逐步衰减，就称该算法相对于算法相对于 绝对稳定绝对稳定， 的全体构成的全体构成绝对稳定区域绝对稳定区域。我
22、们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定区域比的绝对稳定区域比 B 的的大大。h h= =h例：例：考察显式欧拉法考察显式欧拉法由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：0-1-2ReImg例：例：考察隐式欧拉法考察隐式欧拉法可见绝对稳定区域为：可见绝对稳定区域为：210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的好。例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中
23、其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImgk=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定例例1 用欧拉方法用欧拉方法，隐式欧拉方法隐式欧拉方法和和欧拉中点公式欧拉中点公式计算计算的近似解，取步长的近似解，取步长h=0.1，并与精确值比较并与精确值比较解解:欧拉方法的算式为：欧拉方法的算式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉中点公式是两步法，第一步欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式，以后用公式本题精确解为本题精确解为y=e-x,计算结果见表计算结果见
24、表9-1例例2 用用欧拉公式欧拉公式和和梯形公式梯形公式的预估校正法计算：的预估校正法计算：的数值解，取的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较梯形公式只迭代一次，并与精确值比较.方程的解析解为方程的解析解为:解解: 本例中欧拉公式为：本例中欧拉公式为：梯形公式只校正一次的格式为梯形公式只校正一次的格式为结果列入表结果列入表9.21. 1. Runge-Kutta Runge-Kutta 法法的一般形式的一般形式 2. 22. 2阶阶Runge-Kutta Runge-Kutta 方法方法 3. 3. 经典经典Runge-Kutta Runge-Kutta 方法方法 4 4R
25、-K-Fehhlberg R-K-Fehhlberg 方法方法 5. 5. 隐式隐式R-KR-K方法方法 6. 6. 变步长方法变步长方法 龙格库塔法深入研究龙格库塔法深入研究1.1.Runge-Runge-Kutta Kutta 法法的一般形式的一般形式Runge-Kutta 法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为 ：式(9.11) 称L级级p阶阶Runge-Kutta方法方法(简称R-K法法)。当L1就是欧拉法，当L2时为改进的欧拉法。 其中它的局部截断误差是(9.11)2. 22. 2级级2 2阶阶Runge-KuttaRunge-Kutta方法
26、方法令令 L=2L=2，则则 3. 3. 经典经典Runge-KuttaRunge-Kutta方法方法我们可以构造出一族3级3阶，一族4级4阶和一族5级4阶等R-K方法。最常用的4级4阶是如下的经典经典R-KR-K方法方法：4 4R-K-Fehhlberg R-K-Fehhlberg 方法方法R-K-R-K-FehhlbergFehhlberg方法方法是在R-K方法的基础上引进误差和步长控制的办法。即利用5阶R-K法 估计4阶R-K的局部误差，其中 注：注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶
27、R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！ 5. 5. 隐式隐式R-KR-K方法方法类似于显式R-K公式，稍加改变，就得到隐式隐式R-K方法方法。它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-KR-K法有法有:1级级2阶中点公式阶
28、中点公式 :2级级2阶梯形公式：阶梯形公式： 2级级4阶阶R-K公式：公式： 6.6.变步长方法变步长方法在单步法中每一积分步步长实际上是相互独立的，步长的选择具有了灵活性。根据合理地选择每一积分步的步长，既保证精度的要求，又可以减少计算量，从而减小舍入误差。其方便的控制手段是基于误差的事后估计式。对于给定的精度 ,如果 ,反复将步长折半进行计算，直至 为止,这时取最终得到的作为结果；如果 为止，这时再将步长折半一次，就得到所要的结果。这种通过加倍或折半处理步长的计算方法称为变步长方法。 注注：推荐使用精度好计算量低的变步长方法。用四阶显式R-K方法做变步长方法是实践中较好的方法！ 例 分别用
29、改进的欧拉格式和四阶龙格库塔格式解初值问题（取步长h=0.2）：节点 改进欧拉法 四阶龙格库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.1832160.4 1.348312 1.341667 1.3416410.6 1.493704 1.483281 1.4832400.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051表（注：已指出过准确解）单步法的相容性、收敛性和稳定性单步法的相容性、收敛性和稳定性 1.1.单步法的相容性单步法的相容性 2.2.单步法的收敛性单步法的收敛性 3.3.单步法的稳
30、定性单步法的稳定性 4.4.相容性和收敛性的关系相容性和收敛性的关系 5.5.相容性和方法阶的关系相容性和方法阶的关系 6.6.稳定性和收敛性稳定性和收敛性 7.7.绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域 单步法的相容性单步法的相容性定义一定义一：对于（：对于（9.1.1）常微分方程初值问题）常微分方程初值问题单步法的形式可以变表示为单步法的形式可以变表示为 （9.2.19）其中其中 h为步长为步长若对求解区间中任一固定的若对求解区间中任一固定的x，当当 时皆成立时皆成立 则称由（9.2.19）确定的单步法与微分方程初值问题是相容相容的注注意到上式左边极限为由（由（9.1.19.1.1）
31、知它应等于从而由相容性定义得称相容性条件相容性条件。 单步法的收敛性单步法的收敛性定义二：定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19）产生的数值解，对于每一个固定的 , ， 当 即 。若成立 ， 则称该方法是收敛的。单步法的稳定性单步法的稳定性定义三定义三： 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （mn）处的值产生的偏差均不超 过 ，则称该方法是稳定稳定稳定稳定的。 单步法的稳定性有以下定理定理二定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是稳定稳定稳定稳定的。相容性和收敛性的关系相容性和收敛性的关系 定理一：定
32、理一：定理一：定理一： 若单步法的增量函数对变量y满足Lipschtiz 条件，即存在与 h , x 无关的常数 L，对区域D= 任意两点（x，y1）,(x,y2)成立，则单步法收敛的充分必要条件是相容性条件成立。（读者自证）相容性和方法阶的关系相容性和方法阶的关系 若单步法若单步法是是p p阶方法则成立阶方法则成立 若单步法满足相容性条件，得若单步法满足相容性条件，得 所以所以 =0=0也就是说单步法的阶数一定要是正数。由于我们考虑也就是说单步法的阶数一定要是正数。由于我们考虑的单步法皆为正整数，的单步法皆为正整数，p p至少为一。因此我们考虑的至少为一。因此我们考虑的单步法都满足相容性条件
33、。单步法都满足相容性条件。 稳定性和收敛性的关系稳定性和收敛性的关系若单步法的增量函数满足定理二的条件即单若单步法的增量函数满足定理二的条件即单步法是稳定的则单步法收敛的充分必要条件步法是稳定的则单步法收敛的充分必要条件是是相容性条件成立相容性条件成立。 绝对稳定性和绝对稳定绝对稳定性和绝对稳定域域 稳定性问题是一个比较复杂的问题。为了简化讨论一稳定性问题是一个比较复杂的问题。为了简化讨论一般仅对试验方程般仅对试验方程 进行考察。这里假定进行考察。这里假定Re0,Re0h0和的允许范围，称为该方法的和的允许范围，称为该方法的 绝对稳定域绝对稳定域。 绝对稳定性和绝对稳定域绝对稳定性和绝对稳定域
34、2 2将Euler方法应用到试验方程得误差方程是要求误差不增长则则Euler 方法的绝对稳定域是以 为半径的圆的内部。为了保证稳定性步长有所控制。假如当 时h应满足 ，当 时， h 应满足 等等。同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。 1.1.AdamsAdams方法方法 2.2.米尔尼方法、汉明方法及辛普森方法米尔尼方法、汉明方法及辛普森方法 3.3.预测校正方法预测校正方法 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性 5 初值问题的数值解法初值问题的数值解法多步法多步法 考虑型如
35、考虑型如 的的k步法，步法， 称为称为阿当姆斯阿当姆斯（Adams）方法方法1.1.AdamsAdams方法方法拿其中一种为例，推导其公式。对上面所说的法一般形拿其中一种为例，推导其公式。对上面所说的法一般形式式若取若取 ，有，有方程组方程组同样解之，得到同样解之，得到3步步4阶隐式阶隐式Adams公式和它的余项。公式和它的余项。其它请读者自证。我们仅将结论列于下表。其它请读者自证。我们仅将结论列于下表。 Adams显式公式显式公式kp公式11223344Adams隐式公式隐式公式kp公式12233445注：注：在这些方法中，在这些方法中，4阶的阶的Adams预测校正方法具有方法预测校正方法具
36、有方法简洁、使用方便、易排序、高精度等优点。尤其当函数简洁、使用方便、易排序、高精度等优点。尤其当函数f比较复杂时更能显出它的优越性。比较复杂时更能显出它的优越性。同同理理得得到到5个个待待定定参参数数方方程程组组。解解之之得得 ， ， ， ， 。构成著名构成著名的的Miline 4Miline 4步步4 4阶显式公式和它的余项。阶显式公式和它的余项。 ， 2.2.米尔尼方法、汉明方法及米尔尼方法、汉明方法及辛普森方法辛普森方法 同同理理得得到到5 5个个参参数数方方程程组组。求求解解后后就就构构成成著著名名的的3 3步步4 4阶隐式阶隐式HammingHamming公式和它的余项。公式和它的
37、余项。 若若取取，也也得得到到5个个参参数数方方程程组组。求求解解后后就就构构成成Simpson隐式公式和其余项。隐式公式和其余项。 米尔尼方法、汉明方法及米尔尼方法、汉明方法及辛普森方法辛普森方法 不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。 3. 预测校正方法预测校正方法 对于线性多步公式，一般地，不难验
38、证，如果 预测公式是阶或阶精度，校正公式为阶精度， 则用预测公式提供初值，校正公式迭代一次的 效果也能达到阶精度，再迭代下去，效果就不 明显了。预测-校正技术即保证了计算精度，又 使隐式计算显式化，克服了隐式公式要反复迭 代的困难。至于实际计算步长的选取，我们对 预测-校正公式使用外推原理，得到误差估计式 ，用来调整计算步长，使达到控制误差要求。对于线性多步方法常用的预测-校正公式有Miline-Hamming方法和4阶阶Adams显隐式显隐式预测-校正公式 将Miline Miline 公式公式和HammingHamming公式公式结合，构成预测-校正公式 Miline公式公式和Hammin
39、g公式公式的局部截断误差分别为修正修正Miline-Hamming公式公式利用外推原理，即加上后消去局部截断误差的主项。使 说明经过外推后的算法其精度提高一阶。忽略误差项，上式可改写为 由于 和是 在计算过程中获得的数据，称为Miline Miline 公式公式和HammingHamming公式公式的事后误差估计式。我们可以用它们来调节计算步长的大小，即选择一个合适的的步长，使 ，其中的是要求达到的计算精度。 又可得到 MilineMiline公公式式和HammingHamming公公式式的修正公式，它们分别是 从从而而构构成成如如下下的的修修正正HammingHamming预预测测- -校校
40、正正公公式式，简称修正简称修正HammingHamming公式：公式： 预测：修正：校正： 修正： 在应用这套公式时，先由同阶单步法提供初值 ， ， ， 。计算 时可取。 将将4步步4阶阶显显式式Adams公公式式作作为为预预测测公公式式和和3步步4阶阶式式Adams公公式式作作为为校校正正公公式式，构构成成Adams预预测测-校正公式。校正公式。 它们的局部截断误差分别是它们的局部截断误差分别是 Adams预测预测-校正公式校正公式利用外推原理，将上两式作线性组合消去局部截断误差的主项。使计算精度至少提高一阶，同时得到两个修正公式，将它们和上两式构成了如下修正修正Adams公式公式：预测：
41、修正： 校正: 修正： 同理，在计算时，调节计算步长 使 , 由同阶单步法提供初值 , ， ， 。计算 时，可取 。 理论分析得出它们的绝对稳定区域，修正Hamming法： ; 4阶Adams预测预测-校正法校正法： 其中 我们也可以在教学演示上看出修正的预测校正格式的误差非常的小。 多步法的相容性条件多步法的相容性条件 k步法的一般形式为 其中 由对单步法的讨论可知,当 时,方法阶数至少为1。即对 y1,y=x 应精确成立。令 y 1，所以令y=x得 4.4.多步法的相容性、稳定性和收敛性多步法的相容性、稳定性和收敛性所以所以我们称为我们称为线性多步法的相容性条件。线性多步法的相容性条件。
42、k步法需要k个出发值，而初值问题只提供了一个初值，在使用k步法时尚需要其它方法作补充k-1个出发值，今假定它们 是 ，当 这k-1个值都 应收敛于共同的极限y（a）即 在讨论多步法收敛性时我们总假定（9.3.12）成立 定义四定义四： 多步法的收敛问题多步法的收敛问题 的解 收敛于初值问题的解 y(x)。这里 xa+nh. 定义五定义五：称多项式 为k步法的特征多项式。如果特征多项式的零点皆不大于1，且等 于1的零点是单重的，称根条件成立。称多项式 为第二特征多项式。 显然根条件可以表示 定理定理三三：k步法收敛的充要条件为： （1）相容性条件成立。 （2）特征多项式的零点皆不大于1，且等于1
43、的零点是单重的。 （称2为）特征根条件。 多步法的稳定性 关于线性多步法成立以下定理：关于线性多步法成立以下定理： 定定理理四四：若若函函数数f(x,y)f(x,y)对对变变量量y满满足足Lipschtiz 条条件件在在与与h,x无关的常数无关的常数L，对区域对区域D= D= 任意两点任意两点（x x，y1y1）,(x,y2),(x,y2)成立成立 k步法的相容性、收敛性、稳定性有以下关系步法的相容性、收敛性、稳定性有以下关系 对对于于常常微微分分方方程程右右端端函函数数f(x,y)，在在相相容容性性条条件件成成立立情情况况下，下，k步法的收敛性和稳定性是等价的。步法的收敛性和稳定性是等价的。
44、 事事实实上上在在相相容容性性条条件件成成立立时时，收收敛敛性性和和稳稳定定性性的的充充要要条条件都是特征根条件成立。件都是特征根条件成立。多步法的绝对稳定性和绝对稳定域多步法的绝对稳定性和绝对稳定域 将线性多步法的公式应用到试验方程 进行考察。这里假定Re 0,即试验方程本身是稳定的。 记 得 是常系数差分方程，其特征方程为 记它的 k个特征根为 并设它们是互异的。显然根与 有关，不妨记为 注意到当 时 这时由特征方程得 由线性多步法的相容性条件得 是一个根。不妨设, 差分方程的解为 其中系数由线性多步法的出发值确定。另一方面，y(0)=1的试验方程的精确解为 ， 设多步法截断误差为 ，由此
45、可得 ，我们称 为主根，其它根都为增根。 定义五定义五：线性多步法的绝对稳定区域对给定的 ，如果特征方程的特征根 皆按模小于1，则线性多步法关于u是绝对稳定的。使得 成立的 构成绝对稳定区域。注：从误差角度来看绝对稳定区域的方法是一个理想的方法。这样，绝对稳定区域越大，方法适用性越广，因而越优越。 6 6 方程组和刚性方程方程组和刚性方程 9.1 一阶方程组 9.2 化高阶方程为一阶方程组 9.3 刚性方程组9.1 一阶方程组 前面我们研究了单个方程y=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形。 考虑一阶方程组的出值问题 若采用向量的记号，记
46、则常微分方程组的出值问题可以表示为 前几节我们主要讨论了常微分方程的出值问题的数值解法。只要将 y 和f 改写为向量，那么前面提供的各种计算公式即可适用于一阶常微分方程组的出值问题。 Runge-Kutta方法 对于方程组 四级四阶显示Runge- Kutta公式 若写成分量形式就是 i=1,2,m 为了帮助理解这一公式的计算过程，我们再考虑两个方程的特殊情形： 这时四阶龙格库塔公式具有形式 其中 这是一步法，利用节点 上的值 ， ，由上式顺序计 算 ，然后代入即可求得 节点上的 。 9.9.2 2化高阶方程为一阶方程组化高阶方程为一阶方程组关于高阶微分方程的初值问题，原则上总可以归结为一阶方
47、程组来求解，例如，考察下列m 阶微分方程初始条件为只要引进新的变量即可将m 阶方程化为如下的一阶方程组 初始条件则相应的化为9.3 刚性方程组刚性方程组在求解方程组时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性问题。若线性系统 ，其中，中A的特征值 满足条件0(j=1,N), 且 则称 系统为刚性方程，称s为刚性比 7 7 习题与总结习题与总结1 1数值例题数值例题2 2数值练习数值练习 3 3总结总结1、数值例题数值例题 我们已经学习了很多数值算法，他们的效果到底如何呢？我们已经学习了很多数值算法，他们的效果到底如何呢？下面我们来分析一道例题，看看那些方
48、法，就这个问题，下面我们来分析一道例题，看看那些方法，就这个问题，最能接近真实值最能接近真实值 求初值问题求初值问题的数值解，取的数值解，取h=0.1h=0.1并同精确解并同精确解 比较 （1）用欧拉法来计算这个初值问题（2）用各阶的RungeKutta方法来计算这个初值问题（3）用四阶的Adams定步长，Adams定步长预测校正方 法来 计算这个初值问题。然后将数值结果列成表格：自变量 解析解 欧拉法 后退欧拉法 梯形法 改进的欧拉法 0.11.0048 1.0000 1.0091 1.0048 1.0050 0.21.0187 1.0100 1.0264 1.0186 1.0190 0.3
49、1.0408 1.0290 1.0513 1.0406 1.0412 0.41.0703 1.0561 1.0830 1.0701 1.0708 0.51.1065 1.0905 1.1209 1.1063 1.1071 0.61.1488 1.1314 1.1645 1.1485 1.1494 0.71.1966 1.1783 1.2132 1.1963 1.1972 0.81.2493 1.2305 1.2665 1.2490 1.2500 0.91.3066 1.2874 1.3241 1.3063 1.3072 1.01.3679 1.3488 1.3855 1.3676 1.3685
50、 自变量 解析解 4阶显RK 4阶Adams显法 预测校正法 0.11.0048 1.0048 1.0048 1.0048 0.21.0187 1.0187 1.0187 1.0187 0.31.0408 1.0408 1.0408 1.0408 0.41.0703 1.0703 1.0703 1.0703 0.51.1065 1.1065 1.1065 1.1065 0.61.1488 1.1488 1.1488 1.1488 0.71.1966 1.1966 1.1966 1.1966 0.81.2493 1.2493 1.2493 1.2493 0.91.3066 1.3066 1.30
51、66 1.3066 1.01.3679 1.3679 1.3679 1.3679 2、数值练习、数值练习计算实习题：1.使用改进的Euler算法求初值问题 的数值解，取步长h=0.1,并同精确解 比较 2. 利用4阶Runge-Kutta算法求初值问题 的数值解，取步长：(1)h=0.1;(2)h=0.01.求数值解及精确解 3.使用4阶定步长Adams预测校正算法求初值问题 的数值解，取h=0.1,并同精确解 比较 4.使用4阶定步长Adams预测校正算法求初值问题的数值解，使用h=0.5 5综合题综合题请用已学的数值方法来计算下面初值问题，看看那种精度最好，请用已学的数值方法来计算下面初值
52、问题，看看那种精度最好， 的数值解。的数值解。 总总 结结 本章研究的是常微分方程初值问题的数值解法，我们从单步法、多步法两大类来介绍经典算法，并且对相容性、收敛性、稳定性这些理论进行系统的、全面介绍，为同学们学习偏微分方程数值解法作了较好的铺垫。梯形法、R-K方法、变步长法以及Adams方法，都是一些高精度的方法，但是它们也有缺点，四阶R-K方法必须要求f(x,y)有较好的光滑性，若光滑性差，还不如欧拉法或改进欧拉公式。它的另一个缺点是计算量较大，需要好较多的机器时间（每一步需要计算4次f(x,y)），相比之下汉明方法节省计算量，但它是四步方法不是自开始的需要借助四阶R-K方法提供初始值。
53、对数值方法的分析还涉及到局部截断误差，整体误差，相容性、收敛性、稳定性等概念，他们涉及到步长的选取，如步长不合适，舍入误差恶性增长，结果会完全错误。大家在实际应用时一定要注意这些问题。刚性方程组具有很重要的应用价值，它的理论和解法很多，并且还在进一步发展，我们仅是简单介绍一下，详细介绍请查看参考书目。抽取客观事物的本质特征及其内部联系，用适当的数学语言和工具所建立的一个数学结构什么是数学模型什么是数学模型什么是数学模型什么是数学模型例如著名的Malthus人口模型 又如弱肉强食模型例例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算用欧拉方法，隐式欧拉方法和欧拉中点公式计算的近似解，取步长的近似
54、解，取步长h=0.1，并与精确值比较并与精确值比较解解:欧拉方法的算式为：欧拉方法的算式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉隐式方法在本题可解出方程，不必迭代，公式为：欧拉中点公式是两步法，第一步欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式用欧拉公式，以后用公式 利用 Euler方法求初值问题 解解 此时的 Euler公式为例例1的数值解.此问题的精确解是分别取步长h=0.2 ,0.1 ,0.05,计算结果如下hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.4
55、66320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227例2 .用欧拉公式和梯形公式的预估校正法计算：的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:解: 本例中欧拉公式为：
