常微分方程的数值解法

′ = x2 + y2 y
这个一阶微分方程就不能用初等函数及其积 分来表达它的解。 分来表达它的解。
再如， 再如，方程
y ′ = y y (0 ) = 1
虽然有表可查, 的解 y = e x ,虽然有表可查,但对于表 的值, 上没有给出 e x 的值,仍需插值方法来 计算
从实际问题当中归纳出来的微分方程， 从实际问题当中归纳出来的微分方程，通常主要依 靠数值解法来解决。 靠数值解法来解决。本章主要讨论一阶常微分方程 初值问题
x i = x 0 + ih , i = 1, 2, , n
数值解法需要把连续性的问题加以离散化，从而求 数值解法需要把连续性的问题加以离散化， 出离散节点的数值解。 出离散节点的数值解。
对常微分方程数值解法的基本出发点就是离散 化。其数值解法有两个基本特点，它们都采用“步 其数值解法有两个基本特点，它们都采用“ 进式” 进式”，即求解过程顺着节点排列的次序一步一步 地向前推进，描述这类算法， 地向前推进，描述这类算法，要求给出用已知信息
第七章 常微分方程的数值解法
7.1 引言 包含自变量、 包含自变量、未知函数及未知函数的导数或微 分的方程称为微分方程。在微分方程中, 自变量的 分的方程称为微分方程。在微分方程中, 个数只有一个, 称为常微分方程。 个数只有一个, 称为常微分方程。自变量的个数为 两个或两个以上的微分方程叫偏微分方程。 两个或两个以上的微分方程叫偏微分方程。微分方 程中出现的未知函数最高阶导数的阶数称为微分方 程的阶数。如果未知函数y 程的阶数。如果未知函数y及其各阶导数
y i , y i 1 , y i 2 , , y 0 计算 y i + 1
的递推公式。 的递推公式。建立
这类递推公式的基本方法是在这些节点上用数值积 分、数值微分、泰勒展开等离散化方法，对初值问 数值微分、泰勒展开等离散化方法， 题
y ′ = f (x, y) y(x0 ) = y0
进行不同的离散化处理。 中的导数 y ′ 进行不同的离散化处理
本 章 重 点 常微分方程解法回顾 初 值 问 题 讨论一阶常微分方程 y ' = f ( x, y ) 两边积分 通解 y ( x0 ) = y0 x dy 1 2 1 2 = ydy = xdx y + c =- x + c y 2 2 dx y x=0 = 1 y (0) = 1 y2 = x2 + c ∴ y = ± 1 x2
dV
1 = π D 2 dx 4
o
其中x表示高度,直径D是高度x的函数,记为D(x), 其中x表示高度,直径D是高度x的函数,记为D(x), 因此得到如下微分方程初值问题
第七章 常微分方程的数值解法
1 dV = π D 4 dx V (0 ) = 0
2
( x )
( 7.0 )
只要求解上述方程,就可求出体积V与高度x 只要求解上述方程,就可求出体积V与高度x之间的 函数关系,从而可标出容器壁上容积的刻度,但问题 函数关系,从而可标出容器壁上容积的刻度, (7.0)中的函数D(x)无解析表达式 中的函数D(x)无解析表达式, 是(7.0)中的函数D(x)无解析表达式,我们根本无法 求出其解析解. 求出其解析解. 怎样求解不能用解析方法求解的微分方程初制 值问题,就是本章要讨论的重点. 值问题,就是本章要讨论的重点.
在微分方程中, 自变量的个数只有一个, 在微分方程中 自变量的个数只有一个 称为 常微分方程. 常微分方程 自变量的个数为两个或两个以上的微分方程 叫偏微分方程。 叫偏微分方程。
例如
′ y=4x y
2 2 2
一阶线性 二阶非线性
′ y′′-(y) +12xy=0
dy dy + t + y = 0 一阶非线性 dt dt 2 2 2 T T T + 2 + 2 = 0 二阶偏微分方程 2 x y z
H D
0 0.04
0.2 0.11
0.4 0.26
0.6 0.56
0.8 1.04 x
1.0 1.17 1
根据上表的数据,可以拟合出倒 根据上表的数据, 葫芦形状容器的图, 葫芦形状容器的图,建立如图所 示的坐标轴后, 示的坐标轴后,问题即为如何根 据任意高度x标出容器体积V 据任意高度x标出容器体积V的 刻度,由微元思想分析可知 刻度,
V D = π 2
2
H
其中直径D为常数.由于体积V 其中直径D为常数.由于体积V与相对于容器底部的 任意高度H的函数关系明确, 任意高度H的函数关系明确,因此在容器上可以方便 地标出容器刻度,而对于几何形状不是规则的容器, 地标出容器刻度,而对于几何形状不是规则的容器, 比如倒葫芦形状容器壁上如何标出刻度呢? 比如倒葫芦形状容器壁上如何标出刻度呢?下表是 经过测量得到部分容器高度与直径的关系. 经过测量得到部分容器高度与直径的关系.
x2
∫
x
0
e dt , 为计算具体函数值，还需要数值积
t2
分方法。如果需要许多点处的y值，则计算工作量较大。 可见，用求解析的方法来找微分方程的数值解有时是不 适宜的，必须研究微分方程的数值解法来求出其数值解
7.2 数值方法的基本思想 对常微分方程初值问题( 式的数值解法， 对常微分方程初值问题(7.1)式的数值解法，就 式的数值解法 是要算出精确解y(x)在区间[ ] 是要算出精确解y(x)在区间[a,b]上的一系列离散节 y(x)在区间 点 a = x0 < x1 < < x n 1 < x n = b 处的函数值 y ( x 0 ), y ( x1 ), , y ( x n ) 的近似值 y 0 , y 1 , , y n 。相邻两个节点的间距 h = xi + xi 称为步长，步 称为步长， +1 长可以相等，也可以不等。本章总是假定 为定数 为定数， 长可以相等，也可以不等。本章总是假定h为定数， 称为定步长， 称为定步长，这时节点可表示为 定步长
xi p , yi p ( p = 1,2,, k)
步的值， 多步法； ，即前面k步的值，此类方法称为多步法；其代表 即前面 步的值 此类方法称为多步法 是亚当斯法。 是亚当斯法。
7.3 欧拉（Euler）法 欧拉（ ） 7.3.1 Euler公式 公式 欧拉（ 欧拉（Euler）方法是解初值问题的最 ） 简单的数值方法。 简单的数值方法。初值问题
对于初值问题
y ′ = f (x, y) y(x0 ) = y0
的数值解法， 的数值解法，首先要解决的问题就是如何对微分方 程进行离散化，建立求数值解的递推公式。 程进行离散化，建立求数值解的递推公式。递推公 式通常有两类，一类是计算 时只用到x 式通常有两类，一类是计算yi+1时只用到 i+1, xi 和yi， 即前一步的值，因此有了初值以后就可以逐步往下 即前一步的值， 计算，此类方法称为单步法；其代表是龙格—库塔 计算，此类方法称为单步法；其代表是龙格 库塔 龙格 以外， 法。另一类是计算yi+1时，除用到xi+1,xi和yi以外， 另一类是计算y 除用到x 还要用到
1
dy 例 4： = x2 + y2 dx 这是微积分的发明者之一Leibniz在1686年曾经让当时 数学界人士求解的一阶微分方程式，吸引了许多数学家 的注意，大约经过150年的探索，到1838年。刘维尔 （Liouville)在理论上证明了这个微分方程不能用初等 积分法求解，得借助于数值方法 y ′=1-2xy 例5： y(0)=0 的解为 y = e
f ( x, y1 ) f ( x, y 2 ) ≤ L y1 y 2
都成立,则方程( 对R内任意两个 y1 , y 2 都成立,则方程 7.1 )的解 内任意两个 的解 y = y (x ) 在[a, b]上存在且唯一。 ]上存在且唯一。
常微分方程表示方法
2y 例: 方程 xy ′-2y=4x y ′= +4 x 2y f(x,y)= 令： + 4 且给出初值 y(1)=-3 x 就 得到一阶常微分方 程的初值问题 ： 2y dy = f ( x, y) = +4 dx x y (1) = 3
y ′ , y ′′ , , y
(n)
都是一次的,则称它是线性的,否则称为非线性的。 都是一次的,则称它是线性的,否则称为非线性的。
在高等数学中，对于常微分方程的求解，给出 在高等数学中，对于常微分方程的求解， 了一些典型方程求解析解的基本方法， 了一些典型方程求解析解的基本方法，如可分离变 量法、常系数齐次线性方程的解法、 量法、常系数齐次线性方程的解法、常系数非齐次 线性方程的解法等。 线性方程的解法等。但能求解的常微分方程仍然是 有限的，大多数的常微分方程是不可能给出解析解。 有限的，大多数的常微分方程是不可能给出解析解。 譬如
Euler法的求解过程是:从初 法的求解过程是: 法的求解过程是 即点(x ))出发 出发, 始点P 即点 始点 0(即点(x0,y0))出发, 作积分曲线y=y(x)在 作积分曲线y=y(x)在P0点上 y=y(x) 切线 P0 P1 (其斜率为 y′(x0 ) = f (x0 , y0 ) ),与x=x1直线 ),与
可得c=1特解 当x=0时,y=1,可得 时 可得 特解
百度文库
dy = y 2 cos x dy 1 2 = cos xdx = s in x + c dx y y y (0) = 1 1 ∴ y = 可得c=-1特解 当x=0时,y=1,可得 时 可得 特解 1 s in x
dv 例：设降落伞从跳伞塔下落后，所受空气阻力与速度 设降落伞从跳伞塔下落后， = dt 分离变量后得 mg kv m 成正比，并设降落伞离开跳伞塔时（ ） 成正比，并设降落伞离开跳伞塔时（t=0）速度为 dv dt = ∫ 两端积分 求降落伞下落速度与时间的函数关系。 零，求降落伞下落速度与时间的函数关系。 ∫ mg kv m 设降落伞下落速度为v(t),降落伞在空中下落时， 降落伞在空中下落时， 解：设降落伞下落速度为 降落伞在空中下落时 t 1 ln( mg kv ) = + C1 考虑到mg-kv>0 考虑到 同时受到重力P与阻力 的作用，重力大小为mg， 与阻力R的作用 同时受到重力 与阻力 的作用，重力大小为 ，方 k m k t kC 1 向与v一致 阻力大小为kv（ 为比例系数），方向与 一致； 为比例系数）， 向与 一致；阻力大小为 （k为比例系数），方向与 mg kv = e m 即 v相反，从而降落伞所受外力为 相反， 相反 k t mg e kC F=mg-kvm 或 C = v= + Ce k k F=ma 根据牛顿第二定律 (a 为加速度），0得函数v(t)应满足的方程为 mg 为加速度），得函数 应满足的方程为 v t = = 0 代如入上式得 将初始条件 ），得函数 C= dv k m = mg k kv t mg dt v= (1 e m ) 于是所求特解为 = 按题意， 按题意，初始条件为 v t k 0 = 0
y ′ = f (x, y) y(x0 ) = y0
( 7.1 )
在区间a b上的数值解法 在区间a ≤ x ≤ b上的数值解法。 可以证明,如果函数在带形区域 R=a≤x≤b, 可以证明, -∞＜y＜∞}内连续，且关于y满足李普希兹 内连续，且关于y (Lipschitz)条件，即存在常数L(它与x,y无关) (Lipschitz)条件，即存在常数L(它与x,y无关)使 条件 L(它与x,y无关
y ′ = f (x, y) y(x0 ) = y0
的解y=y(x)代表通过点 的解y=y(x)代表通过点 ( x0 , y 0 ) 的一条称之为 y=y(x) 微分方程的积分曲线。 微分方程的积分曲线。积分曲线上每一点 ( x, y ) 的切线的斜率 y′(x) 等于函数 f ( x, y ) 在 这点的值。 这点的值。
Tel： ： 86613747 E-mail： ： lss@zjtcm.net 授课: 授课 68 学分： 学分：4
第七章 常微分方程的数值解法
问题提出 倒葫芦形状容器壁上的刻度问题. 倒葫芦形状容器壁上的刻度问题.对于如图所示 圆柱形状容器壁上的容积刻度, 圆柱形状容器壁上的容积刻度,可以利用圆柱体体积 公式