数值微分法

yn+1 = yn + hf (xn， n )， = 0 1 y n ， ， yn+1 = yn + hf (xn+1， n+1)， = 0 1 y n ， 。 ，
从
x0 处的初值 y0 开始，按（6.1.2）可逐步计算以后各点上的值。称 开始， ）可逐步计算以后各点上的值。 显式Euler。由于（6.1.3）式的右端隐含有待求函数值 yn +1 ， （6.1.2）式为显式 ）式为显式 。由于（ ）
表 6-2
xn
0.1 1.0959 1.0954
0.2 1.1841 1.1832
0.3 1.2662 1.2649
0.4 1.3434 1.3416
0.5 1.4164 1.4142
0.6 1.4860 1.4832
0.7 1.5525 1.5492
0.8 1.6153 1.6165
yn
y(xn)
6.1.2 局部误差和方法的阶
初值问题（ 初值问题（6.1.1）的单步法可以写成如下统一形式 ）
yn+1 = yn + h(xn，xn+1，yn，yn+1，h)，
（6.1.6）
有关。 其中 与 f 有关。若 中不含 y n + 1 则方法是显式的，否则是隐式的，所 ，则方法是显式的，否则是隐式的， 以一般显式单步法表示为 （6.1.7） ） yn+1 = yn + h xn，yn，h 。
n +1 n +1 n +1
由函数y 由函数y(xn+1)在xn处的Taylor展开式： 处的Taylor展开式 展开式：
( xn +1 xn )2 y( xn +1 ) = y( xn ) + ( xn +1 xn ) y′( xn ) + y′(ξ ) n +1 < ξ < xn ) (x 2
y(xn + h) y(xn ) ≈ f (xn， (xn ))， y h y(xn + h) y(xn ) ≈ f (xn+1， (xn+1))。 y h
令 的近似值，将上面两个近似写成等式， yn 为 y ( xn ) 的近似值，将上面两个近似写成等式，整理后得 （6.1.2） ） （6.1.3） ）
（6.1.5） ）
也可以表示为下列平均化的形式
yp = yn + hf ( xn， n ) y ， 1 ( yp + yq ) 。 2 yq = yn + hf (xn+1 yp ) ， ， yn+1 =
例6.2 取h=0.1，用改进的 ，用改进的Euler方法解 方法解
2x ， y y(0) =1 。 y′ = y
不能逐步显式计算， 式为隐式 隐式Euler公式或后退 公式或 公式。 不能逐步显式计算，称（6.1.3 ）式为隐式 公式 后退Euler公式。如果 公式 梯形公式。 将（6.1.2）和（6.1.3）两式作算术平均，就得梯形公式。 ） ）两式作算术平均，就得梯形公式
yn+1 = yn +
h [ f (xn，yn ) + f (xn+1，yn+1)]，n = o， 。 （6.1.4） 1 ， ） 2
yn+1 = 0.1xn + 0.9yn + 0.1 。
同理，用隐式 同理，用隐式Euler方法有 方法有
1 yn+1 = (0.1xn+1 + yn + 0.1)。 1.1
用梯形公式有
yn+1 =
1 (0.1xn + 0.95yn + 0.105)。 1.05
三种方法及准确解 别是 1.4 × 10
yn+1 = yn + hf ( xn， n )， y yn+1 = yn + h [ f (xn，yn ) + f (xn+1，yn+1)]，n = 0， 。 1 ， 2
称该公式为改进的 公式。 称该公式为改进的Euler公式。它显然等价于显式公式为 改进的 公式
yn+1 = yn +
h [ f (xn , yn ) + f (xn+1, yn + hf (xn , yn ))] ， 2
主要问题
如何将微分方程离散化， 如何将微分方程离散化，并建立求其 数值解的递推公式； 数值解的递推公式； 递推公式的局部截断误差， 递推公式的局部截断误差，数值解与 精确解的误差估计； 精确解的误差估计； 递推公式的稳定性分析。 递推公式的稳定性分析。
6.1 Euler 方法
6.1.1 Euler 方法及其有关的方法
d2y S qx y( x ) + ( x l )， = 2 dx EI ( x ) 2 EI ( x ) y（0） y ( L) = 0。 =
针对实际问题建立的数学模型， 针对实际问题建立的数学模型，要找出模型解的解析表达式往往 是困难的，甚至是不可能的。因此， 是困难的，甚至是不可能的。因此，需要研究和掌握微分方程的数值 解法，即计算解域内离散点上的近似值的方法。 解法，即计算解域内离散点上的近似值的方法。本章讨论常微分方程 数值解的基本方法和理论。 数值解的基本方法和理论。
以上两个例子是常微分方程初值问题， 以上两个例子是常微分方程初值问题，下面是一个两点边值问 题的例子。 题的例子。 设一跟长为L的矩形截面的梁 两端固定。 是弹性模量 的矩形截面的梁， 是弹性模量， 是端 设一跟长为 的矩形截面的梁，两端固定。E是弹性模量，S是端 点作用力， （ ）是惯性矩， 是均匀荷载强度 梁的桡度y（ ） 是均匀荷载强度， 点作用力，I（x）是惯性矩，q是均匀荷载强度，梁的桡度 （x）满 足如下方程
梯形公式也是隐式公式。 梯形公式也是隐式公式。以上公式都是由 yn 去计算
yn +1 ，故称它们为单步法。 故称它们为单步法。
方法、 例6.1 取h=0.1，用Euler方法、隐式 ， 方法 隐式Euler方法和梯形方法解 方法和梯形方法解
y ′ = x y + 1， y ( 0 ) = 1。
解 本题有 并代入h=0.1得 得 并代入 方法， 方法 ） f (x，y) = x y +1 y0 = 1 如果用 ， 。如果用Euler方法，由（6.1.2）
微分方程的数值解
设方程问题的解y(x)的存在区间是 的存在区间是[a,b]，令a= x0< 设方程问题的解 的存在区间是 ， 其中h 如是等距节点h=(b-a)/n , h x1<…<xn =b,其中 k=xk+1-xk , 如是等距节点 其中 称为步长。 称为步长。 y(x)的解析表达式不容易得到或根本无法得到，我 的解析表达式不容易得到或根本无法得到， 的解析表达式不容易得到或根本无法得到 们用数值方法求得y(x)在每个节点 k上y(xk)的近似值， 在每个节点x 的近似值， 们用数值方法求得 在每个节点 的近似值 用yk表示，即 yk≈y(xk)，这样 0 , y1 ,...,yn称为微分方程的 表示， ，这样y 数值解。 数值解。
），改进的 解 按（6.1.5），改进的 ），改进的Euler方法解 方法解
yn+1 = yn + h( yn
2xn )， yn
h 2x 2x yn+1 = yn + ( yn n ) + ( yn+1 n+1 )，n = 0， 。 1 ， 2 yn yn+1
得计算结果如表6-2。 由 y0 = 1，h = 0.1 得计算结果如表 。该初值问题的准确解为 y ( x ) = 1+ 2 x 。
dN = αN ( t ) β N 2 ( t ) dt N ( t 0 ) = N 0。
是电容器上的带电量， 为电容 为电容， 为电阻 为电阻， 为电源的 设Q是电容器上的带电量，C为电容，R为电阻，E为电源的 是电容器上的带电量 电动势， 电动势，描述该电容器充电过程的数学模型是
dQ Q (t ) = E ， RC dt Q ( t 0 ) = Q 0。
准确解 1 1.004837 1.018731 1.040818 1.070320 1.106531
对于隐式公式，通常采用估计 校正技术 即先用显式公式计算， 校正技术， 对于隐式公式，通常采用估计-校正技术，即先用显式公式计算，得到 预估值，然后以预估值作为隐式公式的迭代初值， 预估值，然后以预估值作为隐式公式的迭代初值，用隐式公式迭代一次得 到校正值，称为预估 校正技术。例如，用显式Euler公式作预估，用梯形公 预估-校正技术 公式作预估， 到校正值，称为预估 校正技术。例如，用显式 公式作预估 式作校正， 式作校正，即
表6-1
xn
0 0.1 0.2 0.3 0.4 0.5
Euler方法 方法 1 1.000000 1.010000 1.029000 1.056100 1.090490
隐式Euler方法 方法 隐式 1 1.009091 1.026446 1.051315 1.083013 1.120921
梯形法 1 1.004762 1.018549 1.040633 1.070096 1.106278
2
y( x) = x + e x
的数值结果如表6-1所示。 的数Βιβλιοθήκη Baidu结果如表 所示。从表中看 所示 分 。
方法和隐式Euler方法的误差 y ( x n ) y n 到，在 xn = 0.5 处，Euler方法和隐式 方法和隐式 方法的误差
2 和 1.6 × 10 ，而梯形方法的误差却是
2.5 ×10 4
第六章
常微分方程数值解法
引言
科学技术与工程问题常常需要建立微分方程形式的数学模型， 科学技术与工程问题常常需要建立微分方程形式的数学模型， 下面是这类问题的例子。 下面是这类问题的例子。
α 设N（t）为某物种的数量， 为该物种的的出生率与死亡率之 （ ）为某物种的数量， β 为生物的食物供给及它们所占空间的限制， 差， 为生物的食物供给及它们所占空间的限制，描述该物种增长 率的数学模型是
定义6.1. 2 误差为
如果给果给定的算法的截断
Tn +1 = O(h p +1 )
则称该算法具有p阶精度。如果 Tn +1 = g ( xn ,y ( xn ))h p +1 + O(h p + 2 ) 则非零项g ( xn , y ( xn ))h p +1称为为局部截断误差主。
欧拉法的局部截断误差： 欧拉法的局部截断误差：T = y ( x ) y
考虑一阶常微分方程初值的问题： 考虑一阶常微分方程初值的问题： y' = f ( x, y) y( x0 ) = y0 设f（x,y）是连续函数，这样初值问题的解是存在唯一的，而且连续依 （ ）是连续函数，这样初值问题的解是存在唯一的， 赖于初始条件。 赖于初始条件。 为了求得离散点上的函数值，将微分方程的连续问题（ 为了求得离散点上的函数值，将微分方程的连续问题（6.1.1）进行 ） 离散化。一般是引入点列{ 离散化。一般是引入点列 xn }，这里 x n = x n 1 + hn , n = 1, 2 ,...。称 hn ， 1， 为步长，经常考虑定长的情形， 为步长，经常考虑定长的情形，即 hn = h ， x n = x 0 + nh ， n = 0， 。 为初始问题（ 处的值， 记 y ( x n ) 为初始问题（6.1.1）的问题准确解 y ( x ) 在 x n 处的值， ） 用均差近似代替（ 用均差近似代替（6.1.1）的导数得 ）
在例6.1中 由于 （ ， ） 是线性的， 在例 中，由于f（x，y）对y是线性的，所以对隐式公式也可以方便地计 是线性的 但是， （ ， ） 的非线性函数时， 算 y n + 1 。但是，当f（x，y）是y的非线性函数时，如 y ′ = 5 x + 3 y ，其隐式 的非线性函数时 Euler公式为yn+1 = yn + h(5 xn +1 + 3 yn +1 )。显然，它就不是很方便用隐式 公式为 显然，它就不是很方便用隐式Euler方法 方法 解出进行递推计算，此时，可用预测-校正的方法计算 解出进行递推计算，此时，可用预测 校正的方法计算 y n + 1 。
(
)
定义 6.1.1
从初值y ( x 0 ) = y 0出发，由单步法显式或隐式逐步计算，
得x n + 1的值y n + 1 , 则e n + 1 = y ( x n + 1 ) y n + 1 称为在点x n + 1上的整体截断误差。如 果第n步在点x n的值计算没有误差，即y n = y ( x n ), 由单步法计算出 y n + 1 , 则Tn + 1 = y ( x n + 1 ) y n + 1 , 称为点x n + 1上的局部截断误差。