第6章 常微分方程数值解法6.1 6.1 引引 言言6.2 6.2 欧拉方法欧拉方法6.3 6.3 龙格龙格库塔方法库塔方法1详细课资6.1 6.1 引引 言言 微分方程数值解一般可分为：微分方程数值解一般可分为：常微分方程数值解常微分方程数值解和和偏微分偏微分方程数值解方程数值解。自然界与工程技术中的许多现象，其数学表达式。自然界与工程技术中的许多现象，其数学表达式可归结为常微分方程（组）的定解问题。一些偏微分方程问题可归结为常微分方程（组）的定解问题。一些偏微分方程问题也可以转化为常微分方程问题来（近似）求解。也可以转化为常微分方程问题来（近似）求解。Newton最早采最早采用数学方法研究二体问题，其中需要求解的运动方程就是常微用数学方法研究二体问题，其中需要求解的运动方程就是常微分方程。许多著名的数学家，如分方程。许多著名的数学家，如 Bernoulli（家族），（家族），Euler、Gauss、Lagrange和和Laplace等，都遵循历史传统，研究重要等，都遵循历史传统，研究重要的力学问题的数学模型，在这些问题中，许多是常微分方程的的力学问题的数学模型，在这些问题中，许多是常微分方程的求解。作为科学史上的一段佳话，海王星的发现就是通过对常求解。作为科学史上的一段佳话，海王星的发现就是通过对常微分方程的近似计算得到的。本章主要介绍常微分方程数值解微分方程的近似计算得到的。本章主要介绍常微分方程数值解的若干方法。的若干方法。2详细课资1 1、常微分方程与解、常微分方程与解为为n n阶常微分方程阶常微分方程。如果函数如果函数 在区间在区间a,ba,b内内n n阶可导，称方程阶可导，称方程满足方程的函数满足方程的函数称为微分方程的称为微分方程的解解。则则如如为任意常数）为任意常数）一般称为方程的一般称为方程的通解通解。为方程的解。为方程的解。如果如果则有则有为方程满足定解条件的解。为方程满足定解条件的解。一一、初值问题的数值解法初值问题的数值解法3详细课资方程的通解方程的通解满足定解条件的解满足定解条件的解微分关系（方程）微分关系（方程）解的图示解的图示4详细课资本教材重点讨论定解问题本教材重点讨论定解问题( (初值问题）初值问题）定解条件（初始条件）定解条件（初始条件）是否能够找到定解问题的解取决于是否能够找到定解问题的解取决于仅有极少数的方程可以通过仅有极少数的方程可以通过“常数变易法常数变易法”、“可分可分离变量法离变量法”等特殊方法求得初等函数形式的解，绝大等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解。部分方程至今无法理论求解。如如等等等等5详细课资2 2、数值解的思想、数值解的思想（1 1）将连续变量）将连续变量 离散为离散为（2 2）用代数的方法求出解函数）用代数的方法求出解函数 在在 点的近似值点的近似值* *数学界关注数学界关注工程师关注工程师关注如果找不到解函数如果找不到解函数数学界还关注：数学界还关注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振动性解的振动性解的周期性解的周期性解的稳定性解的稳定性解的混沌性解的混沌性6详细课资 求函数求函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值 的方法称为微分方程的数值解法。的方法称为微分方程的数值解法。称节点间距称节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。称为微分方程的数值解。称为微分方程的数值解。所谓数值解法：所谓数值解法：7详细课资称称 在区域在区域D上对上对 满足满足Lipschitz条件条件是指是指:记记3 3、相关定义、相关定义8详细课资(2) 一般构造方法：一般构造方法： 离散点函数值集合离散点函数值集合离散点函数值集合离散点函数值集合 + + 线性组合结构线性组合结构线性组合结构线性组合结构 近似公式近似公式近似公式近似公式4、 迭代格式的构造迭代格式的构造(1) 构构造造思思想想：将将连连续续的的微微分分方方程程及及初初值值条条件件离离散散为为线线性性方方程程组组加加以以求求解解。由由于于离离散散化化的的出出发发点点不不同同，产产生生出出各各种种不不同同的的数数值值方方法法。基基本本方方法法有有：有有限限差差分分法法（数数值值微微分分）、有有限限体体积积法法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。 9详细课资(3) 如何保证迭代公式的稳定性与收敛性如何保证迭代公式的稳定性与收敛性?5、微分方程的数值解法需要解决的主要问题、微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建立求其如何将微分方程离散化，并建立求其数值解的迭代公式？数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？如何估计迭代公式的局部截断误差与整体误差？10详细课资二、初值问题解的存在唯一性二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:则上述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上连续, 且关于且关于 y 满足满足 Lipschitz 条件，条件， 即存在与即存在与 无关的常数无关的常数 L 使使对任意定义在对任意定义在 上的上的 都成立，都成立，11详细课资三三、初值问题的离散化方法初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式，离散化方法的基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右的顺序依次求出按节点从左至右的顺序依次求出 的近似的近似 如果计算如果计算 ，只用到前一步的值，只用到前一步的值 ，则称这则称这类方法为类方法为单步方法单步方法。如果计算如果计算 需用到前需用到前r步的值步的值 , ,则称这类方法为则称这类方法为r步方法步方法。12详细课资6.2 Euler6.2 Euler方法方法第一步：连续变量离散化第一步：连续变量离散化第二步：用直线步进第二步：用直线步进EulerEuler格式格式1 1、EulerEuler格式格式13详细课资l 18 18世纪最杰出的数学家之一，世纪最杰出的数学家之一，1313岁岁时入读巴塞尔大学，时入读巴塞尔大学，1515岁大学毕业，岁大学毕业，1616岁获得硕士学位。岁获得硕士学位。l 17271727年年-1741-1741年（年（2020岁岁-34-34岁）在彼岁）在彼得堡科学院从事研究工作，在分析学、得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、俄国政府要求，解决了不少地图学、造船业等实际问题。造船业等实际问题。l 2424岁晋升物理学教授。岁晋升物理学教授。l 17351735年（年（2828岁）右眼失明。岁）右眼失明。14详细课资l 1741 1741年年 - 1766- 1766（3434岁岁-59-59岁）任德国科学院物理数学所所长，岁）任德国科学院物理数学所所长，任职任职2525年。在行星运动、刚体运动、热力学、弹道学、人口学、年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。微分方程、曲面微分几何等研究领域均有开创性的工作。l 17661766年应沙皇礼聘重回彼得堡，在年应沙皇礼聘重回彼得堡，在17711771年（年（6464岁）左眼失明。岁）左眼失明。l EulerEuler是数学史上最多产的数学家，平均以每年是数学史上最多产的数学家，平均以每年800800页的速度页的速度写出创造性论文。他去世后，人们用写出创造性论文。他去世后，人们用3535年整理出他的研究成果年整理出他的研究成果7474卷。卷。 15详细课资在在假假设设 yi = y(xi)，即即第第 i 步步计计算算是是精精确确的的前前提提下下，考考虑虑的的截截断断误误差差 Ri = y(xi+1) yi+1 称称为为局局部部截截断断误差误差 /* local truncation error */。定义定义2.2 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.12、欧拉法的局部截断误差、欧拉法的局部截断误差16详细课资 欧拉法的局部截断误差：欧拉法的局部截断误差：Ri 的的主项主项/* leading term */欧拉法具有欧拉法具有 1 阶精阶精度度。17详细课资例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题取步长取步长 。 解解: : 应用应用EulerEuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 其中其中 。计算结果列于下表：计算结果列于下表： 18详细课资 19详细课资可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ，从上表最后一列，我们看到取步长从上表最后一列，我们看到取步长20详细课资3、 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数x0x1)(,()(1101xyxfhyxy+ + 21详细课资由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接同时出现在等式的两边，不能直接得到，故称为得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者欧拉公式，而前者称为称为显式显式 /* explicit */ 欧拉公式。欧拉公式。22详细课资一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。23详细课资 梯形公式梯形公式 /*trapezoid formula */*trapezoid formula */ 显、隐式两种算法的显、隐式两种算法的平均平均注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。24详细课资中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数x0x2x1假设假设 , 则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。25详细课资方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度26详细课资 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+ +ny27详细课资28详细课资注注: :此法亦称为此法亦称为预测预测- -校正法校正法 /* predictor-corrector method */* predictor-corrector method */可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显式欧拉法。29详细课资在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 应用改进欧拉法应用改进欧拉法, ,如果序列如果序列 收敛收敛, ,它的极限便满足方程它的极限便满足方程30详细课资因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度31详细课资例例2:2: 用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。解：解：对此初值问题采用改进对此初值问题采用改进EulerEuler公式，公式， 其具体形式为其具体形式为计算结果列于下表：计算结果列于下表：例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题32详细课资改进的改进的Euler法法Euler法法33详细课资通过计算结果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。34详细课资欧拉法误差概述35详细课资6.3 龙格库塔方法 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方析这两个公式的特点，从而探索一条构造高精度方法的途径法的途径. 36详细课资改进欧拉法改进欧拉法37详细课资38详细课资39详细课资 40详细课资41详细课资 42详细课资三阶龙格-库塔方法三阶龙格-库塔方法是用三个值 k1, k2, k3 的线性组合要使三阶龙格-库塔方法具有三阶精度，必须使其局部截断误差为 O(h4)将 k1, k2, k3 代入 yn+1 的表达式中，在 (xn, yn) 处用二元泰勒公式展开，与 y(xn+1) 在 xn 处的泰勒展开式比较43详细课资类似二阶龙格-库塔方法的推导过程，8 个待定系数 c1, c2, c3, a2, a3, b21, b31, b32 应满足： 8 个个未未知知参参数数，6 个个方方程程，有有无无穷穷多多组组解解三阶龙格库塔公式三阶龙格库塔公式44详细课资四阶四阶Runge-Kutta方法方法45详细课资附注：附注：二阶二阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 五阶五阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 四阶四阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 三阶三阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 46详细课资附注：附注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 的值，即计的值，即计算算 的值。的值。Butcher 于于1965年给出了计算量与可年给出了计算量与可达到的最高精度阶数的关系：达到的最高精度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精库塔法的导出基于泰勒展开，故精太好的解，最好采用太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数的光滑性影响。对于光滑性不度主要受解函数的光滑性影响。对于光滑性不47详细课资感谢您的听讲！感谢您的听讲！48详细课资