常微分方程初值问题

( k 1) yn 1 yn 1 .
Euler预估-校正公式（改进Euler公式）：
( p) yn 1 yn hf ( xn , yn ), (c) h ( p) yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )]. 2
2 xn 2 xn1 0.1 ( p) yn 1 yn yn 1 ( p ) , (n 0,1,...,9). yn 2 yn yn 1
x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Euler法y 1.000000 1.000000 1.191818 1.277438 1.358213 1.435133 1.508966 1.580338 1.649783 1.717779 1.784770
hL L为Lipschitz常数。如果步长 h满足 0 1, 2
(k ) { y 则产生的序列 n 1} (k 0,1, 2...) 收敛。
dy 2、将方程 y ' f ( x, y ) 两端在 dx
[ xn , xn 1 ] 上积分，

xn 1 xn
dy dx dx
yn 1 yn hf ( xn 1 , yn 1 ).
隐式Euler公式 1阶方法
隐格式，需迭代求解：
yn 1( k 1) yn hf ( xn 1 , yn 1( k ) ),
y
(0) n 1
可由显式Euler公式求出.
* 定理：设函数 f(x,y)对变量 y满足 Lipschitz条件，
改进欧拉法的预-校公式：
2 yk 1 yk 2hxk yk h 2 2 yk 1 yk (2 xk yk 2 xk 1 y k 1 ) 2
计算得：y1 = 0.99, y2 = 0.9614, y3 = 0.9173. 精确解： y(0.1)=0.99, y(0.2)=0.9614, y(0.3)=0.9173. 可见改进欧拉公式比欧拉公式精度高.
例6.2：用 Euler公式和改进 Euler公式分别求下列 初值问题的数值解（取步长 h=0.1计算到 y3）：
y ' 2 xy 2 . y (0) 1
解：由欧拉公式 yn+1 = yn + hf(xn,yn) = yn - 2hxnyn2 计算如下： y1 = y0 - 2hx0y02 = 1-2·0·0·12 = 1 y2 = y1 - 2hx1y12 = 1-2·0.1·0.1·12 = 0.98 y3 = y2 - 2hx2y22 = 1-2·0.1·0.2·0.982 = 0.9416.
yn 1 yn h(1k1 2 k2 ) k1 f ( xn , yn ) k f ( x h, y hk ) n n 1 2
适当选取 1 , 2 , , , 使局部截断误差
3
Tn1 y( xn1 ) yn1 O(h ).
(*)
为使解存在唯一，一般要对 f 加限制条件， 一个常用的是，f 对 y 满足Lipschitz条件：
f ( x, y1 ) f ( x, y2 ) L y1 y2 .
(**)
数值解：求(*)式的解 y y(x) 在一系列离散节点
y ( xi ) 上的函数值的近似值 yi .
假设对区间[a,b]做等距分割：
由二元函数Taylor展式，
k2 f ( xn , yn ) hf x( xn , yn ) hk1 f y( xn , yn ) O(h2 )
2 y( xn ) h f ( x , y ) f ( x , y ) f ( x , y ) O ( h ) n n y n n x n n
梯形公式 2阶方法
用迭代法求解：
(0) yn 1 yn hf ( xn , yn ), ( k 1) h (k ) y y [ f ( x , y ) f ( x , y n 1 n n n n 1 n 1 )]. 2
( k 1) (k ) when yn y , 1 n
对右端的定积分可用数值积分方法计算. 若用左矩形公式，有：
y( xk 1 ) y( xk ) ( xk 1 xk ) f ( xk , y( xk )),
yk 1 yk hf ( xk , yk ).
显式单步法：
y n 1 y n h ( xn , y n , h ),
第6章 常微分方程
在工程和科学计算中，所建立的各种常微分 方程的初值或边值问题，除很少几类的特殊 方程能给出解析解，绝大多数方程很难甚至 不可能给出解析解，其主要原因在于积分工 具的局限性。 因此，转向用数值方法求常微分方程的数值解。
考虑一阶常微分方程初值问题：
dy f ( x, y), x [a, b] dx 定解条件 y(a) y0 .
线性化，
y( xk 1 ) y( xk ) hy '( xk ),
yk 1 yk hf ( xk , yk ).
• 数值积分方法 将 y’=f(x,y) 在[xk,xk+1]积分，有
y( xk 1 ) y( xk )
xk 1
xk
xk 1 dy dx f ( x, y( x))dx, xk dx
~
~
Tn1 y(xn1 ) yn1 称为点 xn1上的局部截断误差.
误差的累积和传播
定义：若某算法的局部截断误差为 O(hp+1) ， 则称该算法为 p 阶方法。 定理：假设单步法具有 p阶精度，且 f(x,y)关于 y 满足 Lipschitz条件(**)，则整体截断误差为 以显式Euler法为例：
改进Euler法y 1.000000 1.095909 1.184097 1.266201 1.343360 1.416402 1.485956 1.552514 1.616475 1.678166 1.737867
精确解 1.000000 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051
en 1 y ( xn 1 ) yn 1
y ( xn ) yn h f ( xn , y ( xn )) f ( xn , yn ) Tn 1
(1 hL) y( xn ) yn Tn 1
(1 hL) en T , T max T j
yk 1 yk hf ( xk , yk ), k 0,1,, n 1.
（显式）Euler公式
• 利用Taylor展开 设f(x,y)充分光滑，将y(xk+1)在xk处作Taylor展开，
h2 y ( xk 1 ) y ( xk ) hy '( xk ) y '( ), 2
2x dy y y 例6.1：考虑初值问题 dx y (0) 1
分别用Euler法和改进Euler法求解，并与精确解
y 1 2x 进行比较。
解：取 h=0.1，在 [0,1]上计算， 2 xn ), (n 0,1,,9) Euler法： yn1 yn 0.1( yn yn 改进欧拉法： 2 xn ( p) y n1 yn 0.1( yn ), yn
Tn 1 y ( xn 1 ) y n 1
h2 y ( xn ) hf ( xn , y ( xn )) y ''( n ) yn hf ( xn , yn ) 2 yn y ( xn ) 2
O(h ).
p
h y ''( n ) O h 2 . 2
其中函数 与 f(x,y) 有关，称为增量函数. 定义：从初值 y ( x0 ) y0 出发，逐步计算得 xn1的值 yn1,
en1 y( xn1 ) yn1 称为点 xn1上的整体截断误差.
假设到第n步的计算没有误差，即 yn y( xn )，计算 y n 1 ,
(b a) xi a ih, i 0,1,, n, h . n
6.1 Euler法
•微分算子离散化：利用向前差商公式逼近微商
y( xk 1 ) y( xk ) y '( xk ) f ( xk , y( xk )), h
y( xk 1 ) y( xk ) hf ( xk , y( xk )),
yn1 y( xn ) (1 2 )hy( xn )
3 h [2 f x ( xn , yn ) 2 f ( xn , yn ) f y ( xn , yn )] O (h ) 2
与Taylor展开比较，有：
1 2 1 2 1 /2 1 /2 2
6.3 Runge-Kutta法
假设 f(x,y)充分光滑，由Taylor展开， h2 h3 y ( xn 1 ) y ( xn ) hy '( xn ) y "( xn ) y '''( ) 2! 3! 其中， y '( x) f ( x, y ), y ''( x) f x( x, y ) f y( x, y ) f ( x, y ). 所以，可以构造格式 h2 yn1 yn hf ( xn , yn ) f x( xn , yn ) f y( xn , yn ) f ( xn , yn ) 2! 局部截断误差为 O(h3 ) ，是二阶方法； 如何计算偏导数项？
1、二阶R-K公式 改进Euler公式可写为：
yn 1 yn 1 2 f ( k1 k 2 ) (n 0,1, 2,...). k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2
受改进Euler法的启示，设

xn 1 xn
f ( x , y )dx
右端积分用梯形公式计算，
h h3 y ( xn 1 ) y ( xn ) f ( xn , y ( xn )) f ( xn 1 , y ( xn 1 )) f '''( ) 2 12
h yn1 yn [ f ( xn , yn ) f ( xHale Waihona Puke Baidu1 , yn 1 )] 2
e0 0, (1 x)n enx
(1 hL) n 1 e0 T /(hL)
e( n 1) hLT /(hL)
一阶方法
en 1 O(h).
6.2 隐式Euler法和梯形公式
1、采用向后差商公式近似微商
y ( xn 1 ) y ( xn ) h y ' ( xn 1 ) y ' ' ( n ) h 2 2 h y ( xn1 ) y ( xn ) hf ( xn 1 , y ( xn 1 )) y ' ' ( n ) 2
j
en 1 (1 hL)n 1 e0 (1 hL)n (1 hL) 1 T
n 1 1 (1 hL ) (1 hL) n 1 e0 T 1 (1 hL) n 1 (1 hL ) (1 hL) n 1 e0 T hL