第五章 常微分方程初值问题数值解法

y=y(x) Pn
x0
x1
xi
xi+1
xn
这样,从x0逐个算出 x1 , x2 , xn 对应的数值解
y1 , y2 , yn
从图形上看,就获得了一条近似于曲线y=y(x) 的折线 PP 1 2P 3
P n.
Euler格式算法实现 （1）计算步骤 ba ① 输入 a, b, h, f ( x, y), x0 , y0 ,并计算出 N ; h ② 使用Euler格式进行计算
向前推进，描述这类算法，要求给出用已知信息 差分格式 yn , yn1 , yn2 , 计算 yn1 的递推公式. 建立这类递推公式的基本方法是在这些节点上用数 值积分、数值微分、泰勒展开等离散化方法，对初 值问题 y f ( x, y ) y ( x0 ) y 0 中的导数 y 进行不同的离散化处理.
输出 x1, y1 n
n=N ? n+1 n
y 结束
x1 x0 y1 y0
例5.1 用Euler格式解初值问题
2x y y y y(0) 1
取步长h=0.1.
(0 x 1)
Euler格式的具体形式为
yn1 yn hf ( xn , yn ) 2 xn yn 0.1( yn ) yn 0.2 xn 1.1yn yn
y x y
2
2
这个一阶微分方程就不能用初等函数及其积分来 表达它的解.
再如，方程
y y y (0) 1
的解 y e x ,虽然有表可查,但对于表上没 有给出 e x 的值,仍需插值方法来计算
从实际问题当中归纳出来的微分方程，通常主要依 靠数值解法来解决.本章主要讨论一阶常微分方程初 值问题
5.2.2 后退的Euler格式
对 y '( x) f ( x, y( x)),在点xn+1处为
y '( xn1 ) f ( xn1, y( xn1 ))
y ( xn 1 ) y ( xn ) y ( xn 1 ) y ( xn ) 取 y '( xn 1 ) xn 1 xn h
则有 y( xn1 ) y( xn ) hf ( xn1, y( xn1 )) 设y(xn)的近似值yn已知，用它带入上式右端进行计
算，并取计算结果yn+1作为y(xn+1)的近似值,则可得
yn1 yn hf ( xn1 , yn1 )
( 5.3 )
后退的Euler格式
隐式方程(5.3) yn1 yn hf ( xn1 , yn1 )常用迭代法求解 (0) (0) 设用Euler格式 yn 给出迭代初值 y hf ( x , y ) y 1 n n n n 1
长可以相等，也可以不等.本章总是假定h为定数，
称为定步长，这时节点可表示为
,
xn x0 nh, n 0,1, 2,
数值解法需要把连续性的问题加以离散化，从而求 出离散节点的数值解.
对常微分方程数值解法的基本出发点就是离散 化.其数值解法有个基本特点，它们都采用“步进 式”，即求解过程顺着节点排列的次序一步一步地
对R内任意两个 y1 , y2 都成立,则方程( 5.1 )的解 y y ( x) 在a, b上存在且唯一.
数值方法的基本思想 对常微分方程初值问题(5.1)式的数值解法，就 是要算出精确解y(x)在区间a,b上的一系列离散节
点 x1 x2
xn xn1 处的函数值 y( x1 ), y( x2 ), y( xn ), y( xn1 ), 的近似值 y1, y2 , , yn , yn1, 相邻两个节点的间距 h xi 1 xi 称为步长，步
从而有
h2 h3 y( xn1 ) yn 1 y ''( xn ) y '''( xn ) 2 3 hf y ( xn 1 , )( y( xn 1 ) yn 1 ) 移项整理得 1 h2 h3 y( xn1 ) yn1 ( y ''( xn ) y '''( xn ) 1 hf y ( xn1 , ) 2 3 1 由 1 x 2 x3 可知 1 x 1 2 1 hf y ( xn1 , ) (hf y ( xn 1 , )) 1 hf y ( xn1 , )
以 f ( xn , yn )为斜率作直线 y yn f ( xn , yn )(x xn )
当 x xn1 时,得 yn1 yn f ( xn , yn )(x x1 xn ) 取 y( xn1 ) yn1
P1 P0 P1
Pi+1 Pn Pi Pi+1 Pi
程的阶数.如果未知函数y及其各阶导数
y, y,, y
( n)
都是一次的,则称它是线性的,否则称为非线性的.
在《常微分方程》中，对于常微分方程的求解，
给出了一些典型方程求解析解的基本方法，如可分
离变量法、常系数齐次线性方程的解法、常系数非
齐次线性方程的解法等.但能求解的常微分方程仍 然是有限的，大多数的常微分方程是不可能给出解 析解. 譬如
当 x x2 时,得 y2 y1 f ( x1 , y1 )(x2 x1 ) 由此获得了P2的坐标.
P1 P0 P1
Pi+1 Pn Pi Pi+1 Pi
y=y(x) Pn
x0
x1
xi
xi+1
xn
重复以上过程,就可获得一系列的点:P1,P1,…,Pn.对 已求得点 Pn ( xn , y n )
y f ( x, y ) y ( x0 ) y 0 的数值解法，首先要解决的问题就是如何对微分方 程进行离散化，建立求数值解的递推公式. 递推公式通常有两类，一类是计算yn+1时只用到xn+1, xn和yn，即前一步的值，因此有了初值以后就可以逐 步往下计算，此类方法称为单步法. 另一类是计算yn+1时，除用到xn+1,xn和yn以外，还要 用到 xn1, yn1, xn2 , yn2 , xnk , ynk , 即前面k步的值，此 类方法称为多步法.
即有 y '( x) f ( x, y( x))
在点xn处为 y '( xn ) f ( xn , y( xn ))
y( xn1 ) y( xn ) y( xn1 ) y( xn ) 取 y '( xn ) xn1 xn h
y( xn1 ) y( xn ) hf ( xn , y( xn )) yn1 yn 设y(xn)的近似值yn已知，用它带入上式右端进行计 算，并取计算结果yn+1作为y(xn+1)的近似值,则可得
计算公式的精度 常以Taylor展开为工具来分析计算公式的精度. 为简化分析，假定yn是准确的，即在 yn y( xn ) 的前提下估计误差 y( xn1 ) yn1 Euler格式的局部截断误差 由 从而
局部截断误差
f ( xn , yn ) f ( xn , y( xn )) y '( xn ) y( xn1 ) yn1 y( xn1 ) ( yn hf ( xn , yn )) y( xn1 ) y( xn ) hy '( xn ) h2 y ''( xn ) 2
(k ) 如果迭代过程收敛，则极限值 yn 1 lim yn 1 k
必满足隐式方程(5.3),从而获得后退的Euler格式的解.
后退的Euler格式的局部截断误差(假设 yn y( xn ) )
yn1 yn hf ( xn1 , yn1 ) y( xn ) h f ( xn1 , yn1 )
用它带入(5.3)的右端，使之转化为显式，直接计算得
(1) ( 0) yn y hf ( x , y 1 n n1 n1 )
然后再用 y
(1) n 1 ( 2) n1
带入(5.3)的右端，又有 (1) y yn hf ( xn1, yn 1 )
如此反复进行迭代，得
( k 1) (k ) yn y hf ( x , y 1 n n1 n1 ), k 0,1,
y1 y0 hf ( x0 , y0 ) x1 x0 h
③ 输出 x1 , y1，并使 x1 x0 , y1 y0 转到 ② 直至n > N 结束.
（ ） 格 式 的 流 程 图
开始 输入 a,b,h,x0, y0,f(x,y) N=(b-a)/h;n=1
2 Euler
y0+hf(x0,y0 ) y1 x0+h x1
y f ( x, y ) ( 5.1 ) y ( x ) y 0 0 在区间a≤x≤b上的数值解法. 可以证明,如果函数在带形区域 R={a≤x≤b,-∞＜y＜∞} 内连续，且关于y满足Lipschitz条件，即存在常数 L(它与x,y无关)使 f ( x, y1 ) f ( x, y2 ) L y1 y2
对于初值问题
§5.2 Euler方法 5.2.1 Euler格式
Euler方法是解初值问题的最简单的数值方法.
初值问题
y f ( x, y ) y ( x0 ) y 0 的解y=y(x)称为它的积分曲线.积分曲线上每一点 ( x, y ) 处的切线的斜率 y ( x) 等于函数 f ( x, y ) 在这 点的值.
由f ( xn1, yn1 ) f ( xn1, y( xn1 )) f y ( xn1,)( yn1 y( xn1 ))
f ( xn1, y( xn1 )) y '( xn1 )（在xn点Taylor展开） h2 y '( xn ) hy ''( xn ) y '''( xn ) ... 2 3 h 2 因此yn 1 y ( xn ) hy '( xn ) h y ''( xn ) y '''( xn ) 2 hf y ( xn 1 , )( yn 1 y ( xn 1 )) 2 3 h h y ( xn1 ) y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 3!