计算方法第三章

梯形方法关于yn+1是隐式的，而Euler方法是显 式的。一般情形下不容易从上式解出yn+1，因而可将 上式与Euler公式联合使用，即
( yn0)1 yn hf ( xn , yn ) ( k 1) h ( ) yn 1 yn [ f ( xn , y n ) f ( xn 1 , ynk 1 )]( k 0,1,2,; n 0,1,2,) 2
§3.2 泰勒展开法与龙格-库 塔（Runge–Kutta）方法
•问题：利用泰勒展开法推导高阶单步的 求解常微分方程初值问题的数值方法。 •从提高截断误差阶的阶数入手。
一、Taylor 方法
假定初值问题的解y（x）及函数f（x，y） 是充分光滑的，则:
h2 h p ( p) y（xn1） y（xn） hy ' ( xn ) y ' ' ( xn ) y ( xn ) O(h p1 ) 2! p!
的数值解法。 常微分方程初值问题的数值解是求上述初 值问题的解y(x)在区间[a，b]中的点列 xi xi 1 hi i 1， n） （ 2， 上的近似值 yi .以下设 hi 不变，记为h-步长 。
假设解y(x)在区间[a,b]上是存在而且唯一的， 并且具有充分的光滑度，因此，要求f(x,y)也充 分光滑。初值问题的解析解(理论解）用 y( xn ) 表示, 数值解法的精确解用 yn 表示。 常微分方程数值解法一般分为： （1）一步法：在计算yn 1 时，只用到xn 1 ,xn和 y， n 即前一步的值。 （2）多步法：计算 yn 1 时，除用到 xn 1,xn 和 y n 以外，还要用 xn p 和 yn p ( p 1,2 k ; k 0) ,即前 k步的值。 （3）显式格式与隐式格式。
将x=xn+1处该直线上的函数值做为y(xn+1)的近似 值，则有Euler公式。这实质上是在每个小区间上利 用折线来代替曲线的结果，故Euler法又称Euler折线 法。
二、梯形法
在式 y( xn1） y( xn ) x f（t，y（t））dt 中,将积分用 梯形公式来代替，则有
n
xn1
稳定性 Remark：由于稳定性问题比较复杂，通常的 做法是将满足李普希兹条件的微分方程模型 化。设f y==常数，此时微分方程为线性 方程 y= y。为保证微分方程的稳定性， 假定<0。讨论某方法的稳定性，就是讨论 该方法对模型方程的稳定性。
稳定性结论
0 Euler法的稳定性条件是： h
在 yn y ( xn ) 的假定下，
k
1
hf ( xn , yn ) hf ( xn , y ( xn )) hy ' ( xn )
而
k2 hf ( xn h, yn k1 )
hf ( xn h, y ( xn ) k1 ) 2 h f ( xn , y ( xn )) h f ( xn , y ( xn )) k1 f ( xn , y ( xn )) O(h ) x y
而
h2 y ( xn 1） y（xn ) hy （xn） y ' ' ( x n ) O( h 3 ) 2!
因此有
Rn y（x n 1 ） y n 1 O(h 3 )
故Euler预估－校正方法为的局部截断误差阶为O(h3)。
截断误差
定义3：若一个方法的局部截断误差为 O(h ） ， 则称该方法为p阶方法,或称该方法具有p阶精度。 Remark：Euler方法是一阶方法，梯形法和 Euler预估－校正法是二阶方法。
( ( ) 有： ynk 1) yn 1 1 h f ( xn 1 , ynk 1 ) f ( xn 1 , yn 1 ) 1
2 y 1 ( ) hL ynk 1 yn 1 2
2 1 f ( ) h ( xn 1 , )( ynk 1 yn 1 )
hf ( xn , y ( xn )) h f ( xn , y ( xn )) f ( xn , y ( xn )) f ( xn , y ( xn )) O(h 3 ) y x
2
hy ' ( x n ) h 2 y ' ' ( xn ) O(h 3 )
h h3 y ( xn1 ) y ( xn ) [ f ( xn , y ( xn )) f ( xn1 , y ( xn1 )) f ' ' ( , y ( )) ( xn , xn 1 ) 2 12
从而得到梯形公式：
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 ) 2 n 0,1,2,
第三章
常微分方程初值问题的数值解法
§3.0 概述
§3.1 欧拉法与梯形法
§3.2 泰勒展开法与龙格-库塔 （Runge–Kutta）方法
§3.3 线性多步法
§3.4 数值算例
§3.0 概述
本章着重讨论一阶常微分方程初值问题
dy f（x，y） dx y（x0） y0
a xb
§3.1 欧拉法与梯形法
一、欧拉(Euler)法
1， 设节点为 xn x0 nh（n 0，2,3） ，得欧拉方法计算公 式为： y y hf（x ，y ）（n 0，2，） 1，3， n 1 n n n
下面通过几种常用的方法来推导该公式。
1、泰勒展开法
假设在 xn附近把y（x）做Taylor展开，有：
h3 Rn y ( xn 1 ) yn 1 y( n ) 12
n [ xn , xn 1 ]
用泰勒展开法推导Euler预估－校正 格式的局部截断误差
改写Euler预估－校正公式为：
1 1 yn 1 yn k1 k2 2 2 k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2
定义1：n y（xn） yn 为 xn某一数值方法在xn处 的整体截断误差(不考虑舍入误差的影响)。 定义2：对单步法，在 yn y（xn） 的假设下， Rn y（xn1） yn1 称为在 xn处的局部截断误差。
Remark1: Euler法的局部截断误差为:
h2 Rn y ( xn 1 ) yn 1 y( n ) n [ xn , xn 1 ] 2 Remark2: 梯形方法的局部截断误差为:
( yn0)1 yn hf ( xn , yn ) h ( yn 1 yn [ f ( xn , y n ) f ( xn 1 , yn0)1 )]( n 0,1,2,) 2
第一式称为预估公式，第二式称为校正公式。
四、方法的误差估计、收敛性和稳定性
h2 h p ( p1) y( xn ) hf ( xn , y( xn )) f ' ( xn , y( xn )) f ( xn , y( xn )) O(h p1 ) 2! p!
O (h p 1 ) 当n 充分小时，略去余项 ，则有p阶计算公式
h2 h p ( p 1) ( x n , y( xn )) yn1 yn hf ( xn , yn ) f ' ( xn , yn ) f 2 p! y 0
梯形法是绝对稳定的。 Euler预估－校正格式的稳定性条件是：
1 1 h ( h ) 2 1 2
2
对非线性方程，应视 f ，此时将是变化的。 y f 如果步长h固定， y 的变化将引起h的变化， f h h 此时，若 y 属于绝对稳定区域，则认为对 此方程而言，方法是稳定的。
p 1
整体截断误差与局部截断误差的关系
定理：如果f(x，y)满足李普希兹（Lipschitz）条件
f ( x, y1 ) f ( x, y2 ) L y1 y2
且局部截断误差有界： 1 2 | Rn | h M 2 (n 1,2, ) 2 则Euler法的整体截断误差n满足估计式： hM 2 ( ba ) L ( ba ) L n e 0 (e 1) 2L 其中L为李普希兹常数，b-a为求解区间长度，
其中，
y'n f ( xn , yn )
y' 'n f x ( xn , yn ) f ( xn , yn ) f y ( xn , y n )
从而得到Euler公式。
3、数值微分法
y( xn1 ) y ( xn ) xn 1 xn y( xn ) f ( xn , y ( xn ))
4、几何方法
过点(xn,yn)作以f(xn ,yn)为斜率的直线方程：
y yn f ( xn , yn )( x xn )
由
取h的线性部分,并用 y n 表示 y ( xn ) 的近似值,得
2、数值积分法 从xn 到 xn +h对等式 y´(t )=f(t,y(t)) 进行积分得到
y ( xn h) y ( xn )
xn1 xn
f (t , y (t )) dt
再利用左矩形公式，得
y( xn h) y( xn ) hf ( xn , y( xn ))
f h L, 且 L 1 y 2
可以证明,当f(x,y)满足Lipschitz条件,即： (L为Lipschitz常数)时,上述数列收敛。
证明：由 和
h [ ( ) ynk1] yn [ f ( xn , y n ) f ( xn1 , ynk 1 )] 1 2
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 ) ] 2