对两边积分, 得到原方程的解为 : 其中 c 为任意常数. dy dy = y dt y = dt, y(t) = c e t, 从例 8.1 看出, 仅根据常微分方程一般无法得到唯一的解. 要确定唯一解, 还需在一些自变量点上给出未知函数的值, 称为边界条件. 一种边界条件设置方法是给出 t =

Size: px

Start display at page:

Download "对两边积分, 得到原方程的解为 : 其中 c 为任意常数. dy dy = y dt y = dt, y(t) = c e t, 从例 8.1 看出, 仅根据常微分方程一般无法得到唯一的解. 要确定唯一解, 还需在一些自变量点上给出未知函数的值, 称为边界条件. 一种边界条件设置方法是给出 t ="

浚弃玩鄂
7 years ago
Views:

1 第八章常微分方程初值问题的解法在科学与工程问题中, 常微分方程描述物理量的变化规律, 应用非常广泛. 本章介绍最基本的常微分方程初值问题的解法, 主要针对单个常微分方程, 也讨论常微分方程组的有关技术. 8.1 引言行分析. 本节介绍常微分方程以及初值问题的基本概念, 并对常微分方程初值问题的敏感性进问题分类与可解性很多科学与工程问题在数学上都用微分方程来描述, 比如, 天体运动的轨迹机器人控制化学反应过程的描述和控制以及电路瞬态过程分析, 等等. 这些问题中要求解随时间变化的物理量, 即未知函数 y(t),t 表示时间, 而微分方程描述了未知函数与它的一阶或高阶导数之间的关系. 由于未知函数是单变量函数, 这种微分方程被称为常微分方程 (ordinary differential equation, ODE), 它具有如下的一般形式 1 : g(t, y, y,, y (k) ) = 0, (8.1) 其中函数 g: R k+2 R. 类似地, 如果待求的物理量为多元函数, 则由它及其偏导函数构成的微分方程称为偏微分方程 (partial differential equation, PDE). 偏微分方程的数值解法超出了本书的范围, 但其基础是常微分方程的解法. 在实际问题中, 往往有多个物理量相互关联, 它们构成的一组常微分方程决定了整个系统的变化规律. 我们先针对单个常微分方程的问题介绍一些基本概念和求解方法, 然后在第 8.5 节讨论常微分方程组的有关问题. 如公式 (8.1), 若常微分方程包含未知函数的最高阶导数为 y (k), 则称之为 k 阶常微分方程. 大多数情况下, 可将常微分方程 (8.1) 写成如下的等价形式 : y (k) = f(t, y, y,, y (k 1) ), (8.2) 其中函数 f: R k+1 R. 这种等号左边为未知函数的最高阶导数 y (k) 的方程称为显式常微分方程, 对应的形如 (8.1) 式的方程称为隐式常微分方程. 通过简单的变量代换可将一般的 k 阶常微分方程转化为一阶常微分方程组. 例如对于方程 (8.2), 设 u 1 (t) = y(t), u 2 (t) = y (t),, u k (t) = y (k 1), 则得到等价的一阶显式常微分方程组为 : u 1 = u 2 u { 2 = u 3. (8.3) u k = f(t, u 1, u 2,, u k ) 本书仅讨论显式常微分方程, 并且不失一般性, 只需考虑一阶常微分方程或方程组. 例 8.1 ( 一阶显式常微分方程 ): 试用微积分知识求解如下一阶常微分方程 : [ 解 ] 采用分离变量法进行推导 : y = y. 1 为了表达式简洁, 在常微分方程中一般省略函数的自变量, 即将 y(t) 简记为 y,y (t) 简记为 y, 等等.

2 对两边积分, 得到原方程的解为 : 其中 c 为任意常数. dy dy = y dt y = dt, y(t) = c e t, 从例 8.1 看出, 仅根据常微分方程一般无法得到唯一的解. 要确定唯一解, 还需在一些自变量点上给出未知函数的值, 称为边界条件. 一种边界条件设置方法是给出 t = t 0 时未知函数的值 : y(t 0 ) = y 0. 在合理的假定下, 从 t 0 时刻对应的初始状态 y 0 开始, 常微分方程决定了未知函数在 t > t 0 时的变化情况, 也就是说这个边界条件可以确定常微分方程的唯一解 ( 见定理 8.1). 相应地, 称 y(t 0 ) = y 0 为初始条件, 而带初始条件的常微分方程问题 : 为初值问题 (initial value problem, IVP). { y = f(t, y), t t 0 y(t 0 ) = y 0. (8.4) 定理 8.1: 若函数 f(t, y) 关于 y 满足李普希兹 (Lipschitz) 条件, 即存在常数 L > 0, 使得对任意 t t 0, 任意的 y 与 y, 有 : f(t, y) f(t, y ) L y y, (8.5) 则常微分方程初值问题 (8.4) 存在唯一的解. 一般情况下, 定理 8.1 的条件总是满足的, 因此常微分方程初值问题的解总是唯一存在的. 为了更清楚地理解这一点, 考虑 f(t, y) 的偏导数希兹条件 (8.5), 因为 f y f(t, y) f(t, y ) = f (t, ξ) (y y ), y 存在, 则它在求解区域内可推出李普其中 ξ 为介于 y 和 y 之间的某个值. 设 L 为 f y (t, ξ) 的上界,(8.5) 式即得以满足. 对公式 (8.4) 中的一阶常微分方程还可进一步分类. 若 f(t, y) 是关于 y 的线性函数, f(t, y) = a(t)y + b(t), (8.6) 其中 a(t), b(t) 表示自变量为 t 的两个一元函数, 则对应的常微分方程为线性常微分方程, 若 b(t) 0, 则为线性齐次常微分方程. 例 8.1 中的方程属于线性齐次常系数微分方程, 这里的常系数是强调 a(t) 为常数函数问题的敏感性对常微分方程初值问题, 可分析它的敏感性, 即考虑初值发生扰动对结果的影响. 注意这里的结果 ( 解 ) 是一个函数, 而不是一个或多个值. 由于实际应用的需要, 分析常微分方程初值问题的敏感性时主要关心 t 时 y(t) 受影响的情况, 并给出有关的定义. 此外, 考虑到常微分方程的求解总与数值算法交织在一起以及历史的原因, 一般用稳定不稳定等词汇说明问题的敏感性. 定义 8.1: 对于常微分方程初值问题 (8.4), 考虑初值 y 0 的扰动使问题的解 y(t) 发生偏差的情形. 若 t 时 y(t) 的偏差被控制在有界范围内, 则称该初值问题是稳定的 (stable), 否则该初值问题是不稳定的 (unstable). 特别地, 若 t 时 y(t) 的偏差收敛到零, 则称该初值问题是渐进稳定的 (asymptotically stable). 关于定义 8.1, 说明两点 :

3 渐进稳定是比稳定更强的结论, 若一个问题是渐进稳定的, 它必然是稳定的. 对于不稳定的常微分方程初值问题, 初始数据的扰动将使 t 时的结果误差无穷大. 因此为了保证数值求解的有效性, 常微分方程初值问题具有稳定性是非常重要的. 例 8.2 ( 初值问题的稳定性 ): 考察如下模型问题的稳定性 : { y = λy, t t 0 y(t 0 ) = y 0. (8.7) [ 解 ] 易知此常微分方程的准确解为 :y(t) = y 0 e λ(t t 0). 假设初值经过扰动后变为 y 0 + Δy 0, 对应的扰动后解为所以扰动带来的误差为 y (t) = (y 0 + Δy 0 )e λ(t t 0), Δy(t) = Δy 0 e λ(t t 0). 根据定义 8.1, 需考虑 t 时 Δy(t) 的值, 它取决于 λ. 易知, 若 λ 0, 则原问题是稳定的, 若 λ > 0, 原问题不稳定. 而且当 λ < 0 时, 原问题渐进稳定. 图 8-1 分三种情况显示了初值扰动对问题 (8.7) 的解的影响, 从中可以看出不稳定稳定渐进稳定的不同含义. 对例 8.2 中的模型问题, 若考虑参数 λ 为一般的复数, 则问题的稳定性取决于 λ 的实部, 若 Re(λ) 0, 则问题是稳定的, 否则不稳定. 例 8.2 的结论还可推广到线性常系数常微分方程, 即根据 f(t, y) 中 y 的系数可确定初值问题的稳定性. 对于一般的线性常微分方程 (8.6), 由于方程中 y 的系数为关于 t 的函数, 仅能分析 t 取某个值时的局部稳定性. 例 8.3 ( 局部稳定性 ): 考察如下常微分方程初值问题的稳定性 : { y = 10ty, t 0 y(0) = 1. (8.8) [ 解 ] 此常微分方程为线性常微分方程, 其中 y 的系数为 a(t) = 10t. 当 t 0 时,a(t) 0, 在定义域内每个时间点上该问题都是局部稳定的. 事实上, 方程 (8.8) 的解析为 y(t) = e 5t2, 初值扰动 Δy 0 造成的结果误差为 Δy(t) = Δy 0 e 5t2. 这说明初值问题 (8.8) 是稳定的. 对于更一般的一阶常微分方程 (8.4), 由于其中 f(t, y) 可能是非线性函数, 分析它的稳定性非常复杂. 一种方法是通过泰勒展开用一个线性常微分方程来近似它, 再利用线性常微分方程稳定性分析的结论了解它的局部稳定性. 具体的说, 在某个解函数 y (t) 附近用一阶泰勒展开近似 f(t, y), (a) (b) (c) 图 8-1 (a) λ > 0 对应的不稳定问题, (b) λ = 0 对应的稳定问题, (c) λ < 0 对应的渐进稳定问题. f(t, y) f(t, y ) + f y (t, y ) (y y ) 则原微分方程被局部近似为 ( 用符号 z 代替 y):

4 z = f y (t, y ) (z y ) + f(t, y ) 这是关于未知函数 z(t) 的一阶线性常微分方程, 可分析 t 取某个值时的局部稳定性. 因此, 对于具体的 y f (t) 和 t 的取值, 常微分方程初值问题 (8.4) 的局部稳定性取决于 (t, y y ) 的实部的正负号. 应注意的是, 这样得到的关于稳定性的结论只是局部有效的. 实际遇到的大多数常微分方程初值问题都是稳定的, 因此在后面讨论数值解法时这常常是默认的条件. 8.2 简单的数值解法与有关概念大多数常微分方程都无法解析求解 ( 尤其是常微分方程组 ), 只能得到解的数值近似. 数值解与解析解有很大差别, 它是解函数在离散点集上近似值的列表, 因此求解常微分方程的数值方法也叫离散变量法. 本节先介绍最简单的常微分方程初值问题解法欧拉法 (Euler method), 然后给出数值解法的稳定性和准确度的概念, 最后介绍两种隐格式解法欧拉法数值求解常微分方程初值问题, 一般都是步进式的计算过程, 即从 t 0 开始依次算出离散自变量点上的函数近似值. 这些离散自变量点和对应的函数近似值记为 : t 0 < t 1 < < t n < t n+1 < y 0, y 1, < y n, y n+1, 其中 y 0 是根据初值条件已知的. 相邻自变量点的间距为 n = t n+1 t n, 称为步长. 数值解法通常使用形如 y n+1 = G(y n+1, y n, y n 1,, y n k ) (8.9) 的计算公式, 其中 G 表示某个多元函数. 公式 (8.9) 是若干个相邻时间点上函数近似值满足的关系式, 利用它以及较早时间点上函数近似值可算出 y n+1. 若公式 (8.9) 中 k = 0, 则对应的解法称为单步法 (single-step method), 其计算公式为 : y n+1 = G(y n+1, y n ). (8.10) 否则, 称为多步法 (multiple-step method). 另一方面, 若函数 G 与 y n+1 无关, 即 : y n+1 = G(y n, y n 1,, y n k ), 则称为显格式方法 (explicit method), 否则称为隐格式方法 (implicit method). 显然, 显格式方法的计算较简单, 只需将已得到的函数近似值代入等号右边, 则可算出 y n+1. 欧拉法是一种显格式单步法, 对初值问题 (8.4) 其计算公式为 : y n+1 = y n + n f(t n, y n ), n = 0, 1, 2,. (8.11) 它可根据数值微分的向前差分公式 ( 第 7.7 节 ) 导出. 由于 y = f(t, y), 则 y (t n ) = f(t n, y(t n )) y(t n+1) y(t n ) n, 得到近似公式 y(t n+1 ) y(t n ) + n f(t n, y(t n )), 将其中的函数值换为数值近似值, 则得到欧拉法的递推计算公式 (8.11). 还可以从数值积分的角度进行推导, 由于 t n+1 y(t n+1 ) = y(t n ) + y (s)ds = y(t n ) + f(s, y(s))ds, t n t n 用左矩形公式近似计算其中的积分 ( 矩形的高为 s = t n 时被积函数值 ), 则有 t n+1 y(t n+1 ) y(t n ) + n f(t n, y(t n )),

5 将其中的函数值换为数值近似值, 便得到欧拉法的计算公式. 例 8.4 ( 欧拉法 ): 用欧拉法求解初值问题 { y = t y + 1 y(0) = 1 求 t = 0.5 时 y(t) 的值, 计算中将步长分别固定为 0.1 和 [ 解 ] 在本题中,f(t, y) = t y + 1, t 0 = 0, y 0 = 1, 则欧拉法计算公式为 : y n+1 = y n + (t n y n + 1), n = 0, 1, 2, 当步长 h=0.1 时, 计算公式为 y n+1 = 0.9y n + 0.1t n + 0.1; 当步长 h=0.05 时, 计算公式为 y n+1 = 0.95y n t n 两种情况的计算结果列于表 8-1 中, 同时也给出了准确解 y(t) = t + e t 的结果. h=0.1 表 8-1 欧拉法计算例 8.4 的结果. h=0.05 t n y n y(t n ) t n y n t n y n 从计算结果可以看出, 步长取 0.05 时, 计算的误差较小. 在常微分方程初值问题的数值求解过程中, 步长 n, (n = 0, 1, 2, ) 的设置对计算的准确性和计算量都有影响. 一般地, 步长越小计算结果越准确, 但计算步数也越多 ( 对于固定的计算区间右端点 ), 因此总计算量就越大. 在实际的数值求解过程中, 如何设置合适的步长达到准确度与效率的最佳平衡是很重要的一个问题数值解法的稳定性与准确度在使用数值方法求解初值问题时, 还应考虑数值方法的稳定性. 实际的计算过程中都存在误差, 若某一步的解函数近似值 y n 存在误差, 在后续递推计算过程中, 它会如何传播呢? 会不会恶性增长, 以至于淹没准确解? 通过数值方法的稳定性分析可以回答这些问题. 首先给出稳定性的定义. 定义 8.2: 采用某个数值方法求解常微分方程初值问题 (8.4), 若在节点 t n 上的函数近似值存在扰动 δ n, 由它引起的后续各节点上的误差 δ m (m>n) 均不超过 δ n, 即 δ m δ n, (m > n), 则称该方法是稳定的. 在大多数实际问题中, 截断误差是常微分方程数值求解中的主要计算误差, 因此我们忽略舍入误差. 此外, 仅考虑稳定的常微分方程初值问题. 考虑单步法的稳定性, 需要分析扰动 δ n 对 y n+1 的影响, 推导 δ n+1 与 δ n 的关系式. 以欧拉法为例, 先考虑模型问题 (8.7), 并且设 Re(λ) 0. 此时欧拉法的计算公式为 2 : y n+1 = y n + λy n = (1 + λ)y n, 由 y n 上的扰动 δ n 引起 y n+1 的误差为 : δ n+1 = (1 + λ)δ n, 要使 δ n+1 的大小不超过 δ n, 则要求 1 + λ 1. (8.12) 2 对于稳定性分析以及后面的一些场合, 由于只考虑一步的计算, 将步长 n 记为.

常微分方程的数值解法 - Numerical solution of ordinary differential equation

常微分方程的数值解法 - Numerical solution of ordinary differential equation 常微分方程的数值解法 Numerical solution of ordinary differential equation 张晓平 2018 年 12 月 17 日武汉大学数学与统计学院 Table of contents 1. 一般概念 2. 欧拉方法 3. 龙格 - 库塔方法 (Runge-Kutta method) 1 一般概念一般概念 1. 常微分方程的求解问题在实践中经常遇到, 但我们只知道一些特殊类型的常微分方程的解析解