数值分析第9章常微分方程数值解法.ppt

（Numerical Methods for Ordinary Differential Equations ） 问题驱动：蝴蝶效应 洛伦兹吸引子(Lorenz attractor)是由MIT大学的气象学家E dward Lorenz在1963年给出的，他给出第一个混沌现象蝴 蝶效应。 图10.1.1蝴蝶效应示意图 洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组 ： 该方程组来源于模拟大气对流，该模型除了在天气预报中有显 著的应用之外，还可以用于研究空气污染和全球侯变化。洛伦 兹借助于这个模型，将大气流体运动的强度x与水平和垂直方 向的温度变化y和z联系了起来。参数 称为普兰特数， 是规范 化的瑞利数， 和几何形状相关。洛伦兹方程是非线性方程组， 无法求出解析解，必须使用数值方法求解上述微分方程组。洛 伦兹用数值解绘制结果图10.1.1，并发现了混沌现象。 1 1 引引 言言 微分方程数值解一般可分为：常微分方程数值解和偏微分 方程数值解。自然界与工程技术中的许多现象，其数学表达式 可归结为常微分方程（组）的定解问题。一些偏微分方程问题 也可以转化为常微分方程问题来（近似）求解。Newton最早采 用数学方法研究二体问题，其中需要求解的运动方程就是常微 分方程。许多著名的数学家，如 Bernoulli（家族），Euler、 Gauss、Lagrange和Laplace等，都遵循历史传统，研究重要 的力学问题的数学模型，在这些问题中，许多是常微分方程的 求解。作为科学史上的一段佳话，海王星的发现就是通过对常 微分方程的近似计算得到的。本章主要介绍常微分方程数值解 的若干方法。 一、初值问题的数值解法 1、一阶常微分方程初值问题的一般形式 常微分方程的数值解法分为 （1）初值问题的数值解法 （2）边值问题的数值解法 (2) 一般构造方法： 离散点函数值集合离散点函数值集合 + + 线性组合结构线性组合结构 近似公式近似公式 2. 迭代格式的构造 (1) 构造思想：将连续的微分方程及初值条件离散为线性方程 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。 (3) 如何保证迭代公式的稳定性与收敛性 ? 3. 微分方程的数值解法需要解决的主要问题 (1) 如何将微分方程离散化，并建立求其数值解的迭代公式？ (2) 如何估计迭代公式的局部截断误差与整体误差？ 称 在区域D上对 满足Lipschitz条件是指 : 记 4、相关定义 二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 则上述IVP存在唯一解。 只要 在 上连续, 且关于 y 满足 Lipschitz 条件， 即存在与 无关的常数 L 使 对任意定义在 上的 都成立， 求函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值 的方法称为微分方程的数值解法。 称节点间距 为步长， 通常采用等距节点，即取 hi = h (常数)。 称为微分方程的数值解。 三、初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式， 值 ,取 。 按节点从左至右的顺序依次求出 的近似 如果计算 ，只用到前一步的值 ，则称这 类方法为单步方法。 如果计算 需用到前r步的值 , ,则称这类方法为r步方法。 2 2 欧拉方法欧拉方法 /* Eulers Method */* Eulers Method */ 欧拉公式(单步显示公式）： 向前差商近似导数 记为 x0x1 亦称为欧拉折线法 /* Eulers polygonal arc method*/ 在假设 yi = y(xi)，即第 i 步计算是精确的前提 下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断 误差 /* local truncation error */。 定义 若某算法的局部截断误差为O(hp+1)，则称该 算法有p 阶精度。 定义 欧拉法的局部截断误差： Ri 的主项 /* leading term */ 欧拉法具有 1 阶精度 。 例1: 用欧拉公式求解初值问题 取步长 。 解: 应用Euler公式于题给初值问题的具体形式为： 其中 。 计算结果列于下表： 可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。 这个初值问题的准确解为 ， 从上表最后一列，我们看到取步长 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 x0x1 )(,()( 1101 xyxfhyxy + 由于未知数 yi+1 同时出现在等式的两边，不能直接 得到，故称为隐式 /* implicit */ 欧拉公式，而前者 称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。 隐式欧拉法的局部截断误差： 即隐式欧拉公式具有 1 阶精度。 梯形公式 /*trapezoid formula */ 显、隐式两种算法的平均 注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。 但注意到该公式是隐式公式，计算时不得不用到 迭代法，其迭代收敛性与欧拉公式相似。 中点欧拉公式 /* midpoint formula */ 中心差商近似导数 x0x2x1 假设 ,则可以导出 即中点公式具有 2 阶精度。 方 法 显式欧拉 隐式欧拉 梯形公式 中点公式 简单精度低 稳定性最好精度低, 计算量大 精度提高计算量大 精度提高, 显式 多一个初值, 可能影响精度 Cant you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy OK, lets make it possible. 改进欧拉法 /* modified Eulers method */ Step 1: 先用显式欧拉公式作预测，算出 Step 2: 再将 代入隐式梯形公式的右边作校正，得到1+i y 注:此法亦称为预测-校正法 /* predictor-corrector method */ 可以证明该算法具有 2 阶精度，同时可以看到它 是个单步递推格式，比隐式公式的迭代求解过程 简单。后面将看到，它的稳定性高于显式欧拉法 。 改进的欧拉法 在实际计算时，可将欧拉法与梯形法则相结合， 计算公式为 应用改进欧拉法,如果序列 收敛, 它的极限便满足方程 改进欧拉法的截断误差 因此，改进欧拉法公式具有 2 阶精度 例2: 用改进Euler公式求解例1中的初值问题， 取步长 。 解：对此初值问题采用改进Euler公式， 其具体形式为 计算结果列于下表： 改进的Euler法Euler法 通过计算结果得比较可以看出，改进的Euler方法 的计算精度比Euler方法要高。 3 3 龙格龙格 - - 库塔法库塔法 /* /* Runge-KuttaRunge-Kutta Method */ Method */ 建立高精度的单步递推格式 。 单步递推法的基本思想是从 ( xi , yi ) 点出发， 欧拉法及其各种变形所能达到的最高精度为2阶。 以某一斜率沿直线达到 点。 考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？ 步长一定是一个h 吗 ？ 将改进欧拉法推广为 ： ),( ),( 12 1 22111 phKyphxfK yxfK KKhyy ii ii ii += = += + 首先希望能确定系数 1、2、p，使得到的算法格 式有2阶精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 Step 2: 将 K2 代入第1式，得到 Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ，则必须有： 这里有 个未知 数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。 注意到， 就是改进的欧拉法。 为获得更高的精度，应该如何进一步推广？ ).,( ),( ),( ),( . 112211 23213133 12122 1 22111 + += += += = += mm mmmmim ii ii ii mmii hKhKhKyhxfK hKhKyhxfK hKyhxfK yxfK KKKhyy 其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定系数，确定这些 系数的步骤与前面相似。 v考虑一阶常微分方程初值问题 将区域a,b进行分划： 若 则 n级显式Runge-Kutta方法 二阶Runge-Kutta方法 取n=2 记 由此得 另一方面 为使局部截断误差为 ,应取 改进的Euler方法 v取 中点方法 取 二阶Heun方法 取 二级Runge-Kutta方法不超过二阶 v记 则 因此局部截断误差只能达到 三级Runge-Kutta方法 取n=3 记 v又由于 因此要使局部截断误差为 ，必须 Kutta方法 取 三阶Heun方法 v取 三级Runge-Kutta方法不超过三阶 v完全类似于二级Runge-Kutta方法的分析 只能达到 三级Runge-Kutta方法的局部截断误差 将 和 都展开到 项,易证 四级R-K方法 取n=4 局部截断误 差为O(h5) 四阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： 附注： 二阶Runge-Kutta方法的局部截断误差 只能达到 五阶Runge-Kutta方法的局部截断误差 只能达到 四阶Runge-Kutta方法的局部截断误差 只能达到 三阶Runge-Kutta方法的局部截断误差 只能达到 注： 龙格-库塔法的主要运算在于计算 的值，即计 算 的值。Butcher 于1965年给出了计算量与可达 到的最高精度阶数的关系： 753 可达到的最高精度 642 每步须算Ki 的个数 由于龙格-库塔法的导出基于泰勒展开，故精 太好的解，最好采用低阶算法而将步长h 取小。 度主要受解函数的光滑性影响。对于光滑性不 4 4 单步方法的收敛性与稳定性单步方法的收敛性与稳定性 /* /* ConvergencyConvergency and Stability and Stability */ */ 收敛性 /* Convergency */ 定义 若某算法对于任意固定的 x = xi = x0 + i h，当 h0 ( 同时 i ) 时有 yi y( xi )，则称该算法是 收敛的。 例：就初值问题 考察欧拉显式格式的 收敛性。 解：该问题的精确解为 欧拉公式为 对任意固定的 x = xi = i h ，有 稳定性 /* Stability */ 例：考察初值问题 在区间0, 0.5上 的解.分别用欧拉显、隐式格式和改进的欧拉格式 计算数值解。 0.0 0.1 0.2 0.3 0.4 0.5 精确解改进欧拉法 欧拉隐式欧拉显式 节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000101 3.2000101 1.0000 2.5000101 6.2500102 1.5625102 3.9063103 9.7656104 1.0000 2.5000 6.2500 1.5626101 3.9063101 9.7656101 1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107 What is wrong ?! An Engineer complains: “Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!“ 定义 若某算法在计算过程中任一步产生的误差在 以后的计算中都逐步衰减，则称该算法是绝对稳定 的 /*absolutely stable */。 一般分析时为简单起见，只考虑试验方程 /* test equation */ 常数，可以 是复数 我们称算法A 比算法B 稳定，就是指 A 的绝对稳定 区域比 B 的大。 当步长取为 h 时，将某算法应用于上式，并