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

常微分方程的数值解法分为 （1）初值问题的数值解法 （2）边值问题的数值解法 一、初值问题的数值解法
1、 1、一阶常微分方程初值问题的一般形式
y′ = f ( x, y ), y (a) = y0
a≤ x≤b
(1)
2. 迭代格式的构造
(1) 构造思想 ： 将连续的微分方程及初值条件离散为线性方程 构造思想： 组加以求解。由于离散化的出发点不同， 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分） 值方法。基本方法有：有限差分法（数值微分）、有限体积法 数值积分） 有限元法（函数插值）等等。 （数值积分）、有限元法（函数插值）等等。 (2) 一般构造方法： 一般构造方法： 离散点函数值集合 + 线性组合结构 → 近似公式
二、初值问题解的存在唯一性
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
dy = f ( x, y) dx y(a ) = y0 x ∈ [a , b ]
f ( x, y ) 在[ a, b ] × R1上连续 且关于 y 满足 上连续, 只要
y ( xi ) − yi
0.009901 0.018462 0.024153 0.026320 0.025250 0.021852 0.017213 0.012262 0.007626 0.003642 0.000422 0.002053
这个初值问题的准确解为 y ( x ) = 1 (1 + x 2 ) ， 可用来检验近似解的准确程度。 可用来检验近似解的准确程度。 从上表最后一列， 从上表最后一列，我们看到取步长 h = 0.1 进行计算，数值解已达到了一定的精度。 进行计算，数值解已达到了一定的精度。
三、初值问题的离散化方法
离散化方法的基本特点是依照某一递推公式， 离散化方法的基本特点是依照某一递推公式， 按节点从左至右的顺序依次求出 y(xi ) 的近似 值 yi (i = 1, ... , n) ,取 y0 = η 。 取
如果计算 yi +1 ，只用到前一步的值 yi ，则称这 类方法为单步方法。 类方法为单步方法。 单步方法 需用到前r步的值 如果计算yi +1 需用到前 步的值yi , yi −1 ,L, yi −r +1 ,则称这类方法为 步方法。 则称这类方法为r步方法 则称这类方法为 步方法。
4、相关定义 、 记 D = {( x, y ) a ≤ x ≤ b, y ( x ) − δ ≤ y ≤ y ( x) + δ } 在区域D上对 满足Lipschitz条件是指 条件是指 称 f ( x, y ) 在区域 上对 y 满足 条件是指:
∃L > 0
s .t .
f ( x , y1 ) − f ( x , y 2 ) ≤ L y1 − y 2 , ∀ x ∈ [ a , b ], y1 , y 2 ∈ [ y ( x ) − δ , y ( x ) + δ ]
yi
1.000000 0.980000 0.941584 0.888389 0.825250 0.757147 0.688354 0.622018 0.560113 0.503642 0.452911 0.407783
y ( xi )
0.990099 0.961538 0.917431 0.863069 0.800000 0.735294 0.671141 0.609756 0.552486 0.500000 0.452489 0.409836
Ri = y ( xi +1 ) − yi +1
h2 = − y′′( xi ) + O(h3 ) 2
即隐式欧拉公式具有 1 阶精度。 阶精度。
梯形公式 /*trapezoid formula */
— 显、隐式两种算法的平均 隐式两种算法的平均
h yi +1 = yi + [ f ( xi , yi ) + f ( xi +1 , yi +1 )] (i = 0, ... , n − 1) 2
同时出现在等式的两边， 由于未知数 yi+1 同时出现在等式的两边，不能直接 得到，故称为隐式 欧拉公式， 得到，故称为隐式 /* implicit */ 欧拉公式，而前者 称为显式 欧拉公式。 称为显式 /* explicit */ 欧拉公式。
一般先用显式计算一个初值， 迭代求解。 一般先用显式计算一个初值，再迭代求解。 求解 隐式欧拉法的局部截断误差： 隐式欧拉法的局部截断误差： 欧拉法的局部截断误差
§1 引 言
微分方程数值解一般可分为： 微分方程数值解一般可分为：常微分方程数值解和偏微分 方程数值解。自然界与工程技术中的许多现象， 方程数值解。自然界与工程技术中的许多现象，其数学表达式 可归结为常微分方程（组）的定解问题。一些偏微分方程问题 可归结为常微分方程（ 的定解问题。 也可以转化为常微分方程问题来（近似）求解。 也可以转化为常微分方程问题来（近似）求解。Newton最早采 最早采 用数学方法研究二体问题， 用数学方法研究二体问题，其中需要求解的运动方程就是常微 分方程。许多著名的数学家， ），Euler、 分方程。许多著名的数学家，如 Bernoulli（家族）， （家族）， 、 Gauss、Lagrange和Laplace等，都遵循历史传统，研究重要 、 和 等 都遵循历史传统， 的力学问题的数学模型，在这些问题中， 的力学问题的数学模型，在这些问题中，许多是常微分方程的 求解。作为科学史上的一段佳话， 求解。作为科学史上的一段佳话，海王星的发现就是通过对常 微分方程的近似计算得到的。 微分方程的近似计算得到的。本章主要介绍常微分方程数值解 的若干方法。 的若干方法。
记为
y1
亦称为欧拉折线法 亦称为欧拉折线法
/* Euler’s polygonal arc method*/
yi+1 = yi +h f (xi , yi ) (i = 0, ...,n−1)
在假设 yi = y(xi)，即第 i 步计算是精确的前提 ， 称为局部截断 下，考虑的截断误差 Ri = y(xi+1) − yi+1 称为局部截断 误差 /* local truncation error */。 。
中点欧拉公式 /* midpoint formula */
求函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值
y i ≈ y ( xi ) (i = 1, ... , n )
的方法称为微分方程的数值解法。 的方法称为微分方程的数值解法。
y1 ,L , yn 称为微分方程的数值解。 称为微分方程的数值解 数值解。
为步长， 称节点间距 hi = xi +1 − xi (i = 0, ... , n −1) 为步长， 通常采用等距节点， 常数)。 通常采用等距节点，即取 hi = h (常数 。 等距节点 常数
图10.1.1蝴蝶效应示意图 蝴蝶效应示意图
洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组： 洛伦兹方程是大气流体动力学模型的一个简化的常微分方程组：
dx dt = −σ x + σ y dy = rx − y − xz dt dz dt = −bz + xy
2
/* leading term */
=
h2 2
′′( x i ) + O ( h 3 ) y
欧拉法具有 1 阶精 度。
例1: 用欧拉公式求解初值问题
y′ = −2 xy 2 y ( 0) = 1
（0 ≤ x ≤ 1.2）
取步长 h = 0.1。 应用Euler公式于题给初值问题的具体形式为： Euler公式于题给初值问题的具体形式为 解: 应用Euler公式于题给初值问题的具体形式为：
定义
若某算法的局部截断误差为O(hp+1)，则称该 若某算法的局部截断误差为 ， 定义 算法有p 阶精度。 算法有 阶精度。
欧拉法的局部截断误差： 欧拉法的局部截断误差： Ri 的主项
Ri = y( xi +1 ) − yi +1 = [ y( xi ) + hy′( xi ) + h2 y′′( xi ) + O(h3 )]− [ yi + hf ( xi , yi )]
注：的确有局部截断误差 Ri = y ( xi +1 ) − yi +1 = O ( h3 ) ， 即梯形公式具有 阶精度，比欧拉方法有了进步。 即梯形公式具有2 阶精度，比欧拉方法有了进步。 具有 但注意到该公式是隐式公式， 但注意到该公式是隐式公式，计算时不得不用到 隐式公式 迭代法，其迭代收敛性与欧拉公式相似。 迭代法，其迭代收敛性与欧拉公式相似。
该方程组来源于模拟大气对流， 该方程组来源于模拟大气对流，该模型除了在天气预报中有显 著的应用之外，还可以用于研究空气污染和全球侯变化。 著的应用之外，还可以用于研究空气污染和全球侯变化。洛伦 兹借助于这个模型，将大气流体运动的强度x与水平和垂直方 兹借助于这个模型，将大气流体运动的强度 与水平和垂直方 r 称为普兰特数， 向的温度变化y和 联系了起来 联系了起来。 向的温度变化 和z联系了起来。参数σ 称为普兰特数， 是规范 b 化的瑞利数， 和几何形状相关。洛伦兹方程是非线性方程组， 化的瑞利数， 和几何形状相关。洛伦兹方程是非线性方程组， 无法求出解析解，必须使用数值方法求解上述微分方程组。 无法求出解析解，必须使用数值方法求解上述微分方程组。洛 伦兹用数值解绘制结果图10.1.1，并发现了混沌现象。 伦兹用数值解绘制结果图 ，并发现了混沌现象。
（Numerical Methods for Ordinary Differential Equations ）
问题驱动： 问题驱动：蝴蝶效应 洛伦兹吸引子(Lorenz attractor)是由 是由MIT大学的气象学家 大学的气象学家E 洛伦兹吸引子 是由 大学的气象学家 dward Lorenz在1963年给出的，他给出第一个混沌现象 年给出的， 在 年给出的 他给出第一个混沌现象——蝴 蝴 蝶效应。 蝶效应。
3. 微分方程的数值解法需要解决的主要问题
(1) 如何将微分方程离散化，并建立求其数值解的迭代公式？ 如何将微分方程离散化，并建立求其数值解的迭代公式 数值解的迭代公式？ (2) 如何估计迭代公式的局部截断误差与整体误差？ 如何估计迭代公式的局部截断误差与整体误差？ (3) 如何保证迭代公式的稳定性与收敛性 如何保证迭代公式的稳定性与收敛性?
Lipschitz 条件， 条件， 即存在与 x, y 无关的常数 L 使
| f ( x, y1 ) − f ( x, y2 ) | ≤ L| y1 − y2Baidu Nhomakorabea|
都成立， 对任意定义在 [ a, b ] 上的 y1 ( x ) , y2 ( x ) 都成立， 则上述IVP存在唯一解。 存在唯一解。 则上述 存在唯一解
§2
欧拉方法 /* Euler’s Method */
y( x1 ) − y( x0 ) y′( x0 ) ≈ h
欧拉公式(单步显示公式）： 欧拉公式 单步显示公式）： 单步显示公式
向前差商近似导数
x0
x1
y ( x1 ) ≈ y ( x 0 ) + hy ′( x 0 ) = y0 + h f ( x 0 , y0 )
欧拉公式的改进： 欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */ 向后差商近似导数
y′( x1 ) ≈ y( x1 ) − y( x0 ) h
y ( x1 ) ≈ y 0 + h f ( x 1 , y ( x1 ))
x0
x1
yi +1 = yi + hf ( xi +1 , yi +1 )( i = 0, …, n − 1)
yi +1 = yi − 2hxi yi2 ( i = 0,1,...,11) y (0) = 1 计算结果列于下表： 其中 xi = 0.1i。 计算结果列于下表：
i
1 2 3 4 5 6 7 8 9 10 11 12
xi
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2