《数值分析 第9章 常微分方程数值解法》由会员分享，可在线阅读，更多相关《数值分析 第9章 常微分方程数值解法（71页珍藏版）》请在金锄头文库上搜索。
1、（Numerical Methods for Ordinary Differential Equations ）问题驱动：蝴蝶效应洛伦兹吸引子(Lorenz attractor)是由MIT大学的气象学家Edward Lorenz在1963年给出的，他给出第一个混沌现象蝴蝶效应。 图10.1.1蝴蝶效应示意图洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组 ：该方程组来源于模拟大气对流，该模型除了在天气预报中有显 著的应用之外，还可以用于研究空气污染和全球侯变化。洛伦 兹借助于这个模型，将大气流体运动的强度x与水平和垂直方 向的温度变化y和z联系了起来。参数 称为普兰特数， 是规范 化的瑞
2、利数， 和几何形状相关。洛伦兹方程是非线性方程组， 无法求出解析解，必须使用数值方法求解上述微分方程组。洛 伦兹用数值解绘制结果图10.1.1，并发现了混沌现象。 1 1 引引 言言微分方程数值解一般可分为：常微分方程数值解和偏微分 方程数值解。自然界与工程技术中的许多现象，其数学表达式 可归结为常微分方程（组）的定解问题。一些偏微分方程问题 也可以转化为常微分方程问题来（近似）求解。Newton最早采 用数学方法研究二体问题，其中需要求解的运动方程就是常微 分方程。许多著名的数学家，如 Bernoulli（家族），Euler、 Gauss、Lagrange和Laplace等，都遵循历史传统，
3、研究重要 的力学问题的数学模型，在这些问题中，许多是常微分方程的 求解。作为科学史上的一段佳话，海王星的发现就是通过对常 微分方程的近似计算得到的。本章主要介绍常微分方程数值解 的若干方法。一、初值问题的数值解法1、一阶常微分方程初值问题的一般形式常微分方程的数值解法分为 （1）初值问题的数值解法 （2）边值问题的数值解法(2) 一般构造方法：离散点函数值集合离散点函数值集合 + + 线性组合结构线性组合结构 近似公式近似公式2. 迭代格式的构造(1) 构造思想：将连续的微分方程及初值条件离散为线性方程组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值
4、微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。 (3) 如何保证迭代公式的稳定性与收敛性 ?3. 微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建立求其数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？称 在区域D上对 满足Lipschitz条件是指 :记4、相关定义二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:则上述IVP存在唯一解。只要 在 上连续, 且关于 y 满足 Lipschitz 条件， 即存在与 无关的常数 L 使对任意定义在 上的 都成立，求函数 y(x)
5、在一系列节点a = x0 x1 xn= b 处的近似值 的方法称为微分方程的数值解法。称节点间距 为步长，通常采用等距节点，即取 hi = h (常数)。称为微分方程的数值解。三、初值问题的离散化方法离散化方法的基本特点是依照某一递推公式，值 ,取 。按节点从左至右的顺序依次求出 的近似如果计算 ，只用到前一步的值 ，则称这类方法为单步方法。如果计算 需用到前r步的值 , ,则称这类方法为r步方法。2 2 欧拉方法欧拉方法 /* Eulers Method */* Eulers Method */ 欧拉公式(单步显示公式）：向前差商近似导数记为x0x1亦称为欧拉折线法 /* Eulers po
6、lygonal 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公式于题给初值问题的具体形式为： 其中 。计算结果列于下表： 可用来检验近似解的准确程度。进行计算，数值解已达到了一定的精度。这
7、个初值问题的准确解为 ，从上表最后一列，我们看到取步长 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */向后差商近似导数x0x1)(,()(1101xyxfhyxy+由于未知数 yi+1 同时出现在等式的两边，不能直接 得到，故称为隐式 /* implicit */ 欧拉公式，而前者 称为显式 /* explicit */ 欧拉公式。一般先用显式计算一个初值，再迭代求解。 隐式欧拉法的局部截断误差：即隐式欧拉公式具有 1 阶精度。 梯形公式 /*trapezoid formula */ 显、隐式两种算法的平均注：的确有局部截断误差 ，即梯形公式具有2 阶精度
8、，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。 中点欧拉公式 /* midpoint formula */中心差商近似导数x0x2x1假设 ,则可以导出 即中点公式具有 2 阶精度。方 法 显式欧拉 隐式欧拉 梯形公式中点公式简单精度低 稳定性最好精度低, 计算量大精度提高计算量大精度提高, 显式多一个初值, 可能影响精度Cant you give me a formula with all the advantages yet without any of the disadvantages?Do you think it possibl
9、e?Well, call me greedyOK, lets make it possible. 改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出Step 2: 再将 代入隐式梯形公式的右边作校正，得到1+iy注:此法亦称为预测-校正法 /* predictor-corrector method */ 可以证明该算法具有 2 阶精度，同时可以看到它 是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法 。 改进的欧拉法在实际计算时，可将欧拉法与梯形法则相结合， 计算公式为 应用改进欧拉法,如果序列
10、收敛,它的极限便满足方程改进欧拉法的截断误差因此，改进欧拉法公式具有 2 阶精度例2: 用改进Euler公式求解例1中的初值问题， 取步长 。解：对此初值问题采用改进Euler公式，其具体形式为计算结果列于下表：改进的Euler法Euler法通过计算结果得比较可以看出，改进的Euler方法的计算精度比Euler方法要高。3 3 龙格龙格 - - 库塔法库塔法 /* /* Runge-KuttaRunge-Kutta Method */ Method */建立高精度的单步递推格式 。单步递推法的基本思想是从 ( xi , yi ) 点出发，欧拉法及其各种变形所能达到的最高精度为2阶。以某一斜率沿
11、直线达到 点。 考察改进的欧拉法，可以将其改写为：斜率 一定取K1 K2 的平均值吗？步长一定是一个h 吗 ？将改进欧拉法推广为 ：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+首先希望能确定系数 1、2、p，使得到的算法格 式有2阶精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开Step 2: 将 K2 代入第1式，得到Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较要求 ，则必须有：这里有 个未知 数， 个方程。3 2存在无穷多个解。所有满足上式的格式统称
12、为2阶 龙格 - 库塔格式。注意到， 就是改进的欧拉法。 为获得更高的精度，应该如何进一步推广？).,(),(),(),(.1122112321313312122122111+=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy其中i ( i = 1, , m )，i ( i = 2, , m ) 和 ij ( i = 2, , m; j = 1, , i1 ) 均为待定系数，确定这些 系数的步骤与前面相似。 v考虑一阶常微分方程初值问题将区域a,b进行分划：若则n级显式Runge-Kutta方法二阶Runge-Kutt
13、a方法取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 *
14、/ ：附注：二阶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 )，则称该
