1、第九章 常微分方程初值问题的数值解法,科学研究和工程技术中的问题往往归结为求某个常微分方程的定解问题.常微分方程的理论指出，很多方程的定解虽然存在，但可能十分复杂难于计算，也可能不能用简单的初等函数表示,因此常求其能满足精度要求的近似解.,1 引 言,常微分方程的数值解法常用来求近似解，由于它提供的算法能通过计算机便捷地实现，因此近年来得到迅速的发展和广泛的应用。常微分方程数值解法的特点是：对求解区间进行剖分，然后把微分方程离散成在节点上的近似公式或近似方程，最后结合定解条件求出近似解.因此数值解法得到的近似解是一个离散的函数表.,1,设方程问题的解 y(x) 的存在， 令a= x0 x1 x
2、n 其中 hk =xk+1 - xk , 如是等距节点 xk=a+kh , h 称为步长。 y(x) 的解析表达式不容易得到或根本无法得到，我们用数值方法求得 y(x)在每个节点 xk 上 y(xk) 的近似值，用 yk 表示，即yk y(xk)，这样 y0 , y1 ,., yn 称为微分方程的数值解。,如方法中要计算yn+1 时只用到前一点的值 yn 则称为单步法. 又如需要前k点值 yn , yn-1, yn-k+1 来计算 yn+1 称其为k步法. 其次我们还得研究算法的局部误差和阶, 数值解与精确解间 的误差及收敛性,还有数值算法的稳定性.,微分方程的数值解：,2,基本知识复习 一阶
3、常微方程的初值问题：,其中f (x,y)是已知函数，(2)是定解条件，常微分方程的理论指出：当f (x,y)定义在区域G=(axb,y)，若存在正的常数L，使： f (x,y1) f (x,y2) Ly1- y2|,对任意的 x a, b 和任意的y1, y2 成立，则称 f (x , y)对 y 满足利普希茨(Lipschitz)条件, L称利普希茨常数. 若f ( x, y) 在区域G连续，关于y 满足利普希茨条件，则上述一阶常微分方程的初值问题的解存在且唯一. 我们以下的讨论，都在满足上述条件下进行.,3,一阶常微分方程组和高阶常微分方程,一阶常微分方程组常表述为：,写成向量形式为,4,
4、还有高阶常微分方程定解 问题，如二阶定解问题：,可化作一阶常微方程组的 定解问题：,本章着重讨论一阶常微方程初值问题的各种数值解法. 这些解法都可以写成向量形式而用于一阶常微分方程组的初值问题.这样，也就解决了高阶方程的 定解问题.,5,返回,2 简单的数值方法与基本概念,对问题,最简单而直观的方法是欧拉方法.欧拉方法在精度要求不高时，仍不失为一实用方法.同时通过欧拉方法的讨论，容易弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.,欧拉方法 (Euler) 的导出,取步长为 h，节点 xn=x0+ nh, n = 0,1,2,，其中 x0=a 又设 y (x)为上述问题的解.,对上
5、式两端积分，积分区间取xn, xn+1,6,记g (x)=f (x, y (x)对上述 (*) 积分 用左矩形公式,得到：,舍去 Rn，并令 yn= y(xn),得到：,这就是欧拉公式.,7,欧拉公式还可由泰勒展开得到. 假定y (x)二阶连续可导，把 y(xn+1)在 xn 点展开：,舍去h2项，并令 yn= y( xn) ，注意到y(xn) = f (xn, yn), 也得到,欧拉公式有明显的几何意义，如图9-1(P337)，过点(x0, y0)的曲线是解 y(x).欧拉方法是在(x0, y0)作 y(x) 的切线，它与直线x=x1 交于 (x1, y1)，过(x1, y1)作过此点的积分
6、曲线的切线，又与x=x2交于(x2, y2)，如此下去，得到一条折线，欧拉方法就是用这条折线近似地代替曲线y(x)，故欧拉方法有时也称欧拉折线法.,8,推出：,不难看出， (2.5) 式中公式两端都含有 yn+1, 一般情况不能由 yn 的值计算 yn+1, 而需要解方程， 故称为 欧拉隐式公式,隐式方法一般用迭代方式求解: 设由欧拉方法得:,代如(2,5)式得:,如此有:,对(*)式的右端用右矩形公式代替积分，可得到,欧拉隐式公式和欧拉中点公式,(后退的欧拉公式),9,推出：,称为欧拉中点公式(二步法).容易看出,中点公式计算 yn+1时,不仅需要 yn 的值，还需要 yn-1的值.,又如把
7、积分区间改为xn-1, xn+1，并对右端用中矩形公式代替积分，可得到：,利用利普希茨条件: 由 (2.5) 减 (2.6) 式得:,由此只要 hL1 上述迭代就可收敛到 y n+1,10,凡是从已知(或已算出)的 y0,y1,yn 能直接从公式算出 yn+1的公式称显式，否则称隐式公式. 在计算 yn+1 时，只需要 yn 的值，则称公式为单步法；若除 yn 之外还需要以前的 yn-1 等多个值，则称多步法公式.,局部截断误差和方法的阶,初值问题的单步法可用一般形式表示成:,其中 称为增量函数, 与 f (x, y) 有关, 当含有yn+1 时, 方法是隐式的. 故显式为,11,对以上三个公
8、式. 用公式计算yn+1时产生的局部截断误差:,对欧拉公式：,即: 当 y(xn) = yn 时,对欧拉隐式公式：,为显式单步法(2.10)的局部截断误差,定义: 设 y(x) 是初值问题的准确解, 假定 yn+1以前的 yn 值都 没有误差，称:,12,对欧拉中点公式的局部截断误差：,定义: 凡是局部截断误差为O(hp+1)的方法称为具有 p阶精度的.这样欧拉方法和欧拉隐式方法称为一阶方法，而欧拉中点方法称为二阶方法. 一般讲高阶方法比低阶方法要精确.,梯形公式及其预估校正法,对积分式(*)的右端用梯形积分公式近似代替，得,13,把 y(xn+1)记为 yn+1，并略去余项 -h3y()/1
9、2得到梯形公式,它是单步隐式方法，也是一个二阶方法. 又是隐式方法，一般求yn+1需要解方程，常用迭代法解此方程，初值y(0)n+1由显式的欧拉公式给出,称作改进的欧拉法, 完整的公式写为:,14,它称为梯形公式的预估校正法，用低精度的显式欧拉公式给出下一步 yn+1的初值 y(0)n+1，再用较高精度的隐式方法梯形公式进行迭代,一般迭代一两次就能达到精度要求.迭代的收敛与否与步长 h 有关系，这里迭代函数是:,所以 h 要适当地小,才能得到收敛序列:,15,或将(2.7)(2.8)相减可得:,于是有,所以当h选取充分小,且使,则当,说明迭代法(2.8)是收敛的.,如迭代不收敛，则应减小步长h
10、 重新计算,16,例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算,的近似解，取步长h=0.1，并比较精确解 y=e x ,解:欧拉方法的算式为：,欧拉隐式方法在本题可解出方程，不必迭代，公式为：,欧拉中点公式是两步法，第一步yn用欧拉公式，以后用公式,17,例2 用欧拉公式和梯形公式的预估校正法计算：,的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为:,解: 本例中欧拉公式为：,计算结果见表9-1(P338),18,梯形公式只校正一次的格式(改进的欧拉法)为:,结果列入表9.2(P344),19,返回,20,3 龙格-库塔(Runge-Kutta)法,欧拉方法是显式
11、的一步法, 使用方便, 但精度较低.我们将构造出高精度的显式一步法：龙格-库塔法, 简称R-K法.,前面所得的改进的欧拉法可表示为:,显式龙格-库塔(Runge-Kutta)法的一般形式,增量函数,由欧拉法的公式： yn+1= yn+h f (xn ,yn) n=0,1,2,为此对一般公式:,选择函数 (xn ,yn,h), 一种方法是用若干个点的函数值的线性组合代替 (xn ,yn, h), 如：,决定其精度的是函数 f (xn , yn). 改进成上述增量函数, 就可 提高公式的精度.,21,其中ci , i , ij 是待定参数, aj和bjl满足,以上方法称为 r 级R-K法, 选择c
12、i , i , ij, 可能使以上方法为r 阶方法.,显然欧拉法就是一阶R-K法.,二级R-K法的形式是:,22,当令y(xn) = yn 时, 此时,由二元函数的泰勒展开：,其中所有的偏导数都是它们在点(xn,y(xn)的值, 下同,又由于：,23,Taylor展开式,(3.6)与上式相减,得局部截断误差,24,只要 c2 , 2 , 21, 满足以上方程, 就得到一个二阶的R-K法. 这是一个不定方程, 有无穷多解. 比如：,(1)取 c1 = c2 = 1/2, 2 = 21 =1 得,这实际上是(2.8)公式,即梯形公式的预估-校正公式只迭代一次的形式,通常称为改进的欧拉法.,25,(
13、2) 取c1=0, c2=1, 2 = 21 = 1/2 得:,这公式又称中点公式. 我们还可以构造其他的二阶R-K法.,三阶与四阶 显式R-K 法,用类似的方法可以确定三级和四级R-K法的参数, 构造出 三阶和四阶的R-K法. 同样方法的形式也是不唯一的.,26,常用的三阶R-K公式,n=1,2,同样经过复杂的数学演算,可得各种四阶R-K公式, 其局部误差为 Tn+1 = O(h5),常用的有:,27,当r=3时 得 三阶显式 R-K方法,局部误差,将 K1 K2 按二元泰勒展开, 使 Tn+1 = O(h4)可得未知参数所 满足的方程, 但其解也不是唯一的(P348),28,例3 用经典的
14、四阶R-K法计算 ,取步长为0.2, 且与准确值比较.,计算结果列入表9-3：可见 即使用h=0.2计算,也比一阶 和二阶方法精度好得多(P349) ,29,30,变步长的R-K方法,步长的选择决定了方法的精度,和计算中舍入误差积累的可能性,选择步长需考虑二个问题:,怎样衡量和检验计算结果的精度. 如何依据所获得的精度处理步长.,考虑四阶R-K公式,从 xn 以 h 为步长求出一个近似值 yn+1(h),由局部误差为 Tn+1 = O(h5) 故有:,再以 h/2 为步长计算二步到 xn+1 节点出的近似值 yn+1(h/2),则有:,31,由此可得:,得事后估计式,这样我们就可以通过检查步长
15、折半前后的计算结果偏差,与预先给定的精度来判断选取合适的步长, 这样的处理步长的方法称为变步长方法.,32,返回,4 单步法的 收敛性与稳定性,收敛性与相容性,数值解法的基本思想就是通过某种离散手段将初值问题 转化为差分方程,则,称为整体局部误差,定义:单步法当任意固定 x = xn ,取 h = (x-x0) / n 时,有,则称该方法是收敛的,33,定理1 假设单步法(4.1)具有p阶精度, 且增量(x,y,h)关于y 满足利普希茨条件 (x,y1,h) (x,y2,h) L y1- y2| 又设初值y0是准确的, 即y0=y(x0), 则其整体误差为,证明: 设,故:,34,(4.2)
16、.,反复递推可得:,注意到,得估计式,故当y0=y(x0), 则其整体误差为,依据定理1, 判断单步法的收敛性, 可归结为验证增量函数是否满足利普希茨条件.,35,对于欧拉方法,因为其增量函数就是 就是 f(x,y) ,所以当 f (x, y)关于 y 满足利普希茨条件,它即为收敛的.,对于改进的欧拉方法,因为其增量函数,所以,由 f(x, y)关于y满足利普希茨条件, 得:,设限定,则关于y的利普希茨常数为,36,类似不难验证其他R-K方法的收敛性. 如四阶R-K法的增量函数关于y的利普希茨常数为:,定理1表明当(4.1)至少为一阶时单步法收敛,且当y(x)是初值 问题(1.1)(1.2)的
17、解,(4.1)具有 p 阶精度时,则有展开式:,37,所以当p1的充要条件是,又,由此得如下定义,定义: 若单步法(4.1)的增量函数满足,则称单步法(4.1)与初值问题(1.1)(1.2)相容,易得: P阶精度方法(4.1)当p1时,与(1.1)(1.2)相容 反之相容方法至少是1阶的,方法(4.1)收敛的充分必要条件是此方法是相容的,结合定理1可知:,38,前面所讨论的收敛性都是在数值方法本身的计算是准确的假设条件下下进行的,由于计算方法都是递推的.初值y0就可能含有误差,在计算y1,y2,yn时又会发生舍入误差得y*n,称差值 n= y*n yn 第n步数值解的扰动.,绝对稳定与绝对稳定
18、域,稳定性与具体的方程有关,为比较各方法的稳定性,常选用 试验方程讨论稳定性问题,试验方程常取为,定义: 如一种数值方法在节点值yn上有扰动n,与以后各节点值ym (mn)上产生的偏差均不超过n, 即 m n 则称 该方法是稳定的. 否则称此方法是不稳定的.,39,定义: 单步法(4.1)用于试验方程(4.8),如果得到的解 yn+1=E(h)yn ,满足E(h)1, 则称方法(4.1)是绝对稳定的. 在=h的复平面上, 使E(h)1的变量围成的区域, 称为 绝对稳定区域,它与实轴的交称为绝对稳定区间.,单步法(4.1)用于试验方程(4.8)得到的解形式yn+1=E(h)yn 故得扰动方程 n
19、+1=E(h)n,Euler法的绝对稳定区域,从而,绝对稳定区域越大,h可选大些,方法的适应性就越强.,后退的Euler 方法稳定性,得,只要 R(h) 0, 则有 E(h) 1,故后退的Euler 方法是无条件稳定性的,但收敛性要求 0 h 1,h仍受限制.,41,梯形公式的稳定性,所以梯形公式是无条件绝对收敛的.,42,四阶R-K方法的绝对稳定区域,43,所以,所以绝对稳定区域为:,44,返回,45,5 线性多步法,单步法只利用前一步的结果,只要给出初值,就能开始计算但也因为它只利用前一步的值,为了提高精度就要计算一些非结点处的函数值,增加了计算量. R-K法就是通过这一途径提高精度的.下
20、面介绍的线性多步法,在求yn+1时,不仅用到yn的值, 还用到前若干步的yn-1,yn-k的值,这些值都是已知的,因此可在计算量增加不多的情况下提高精度.形成线性多步法.,用待定系数法构造线性多步法,线性k步法的一般形式是：,46,其中i , i 都是常数, 0+00,y n+i为y(xn+i) 近似值, fn+j= f (xn+j, yn+j), xn+i = x0+ih,由(5.1)可看出要计算yn+k,要利用它前面的k个值yn,yn+1,yn+k-1, 又因为(5.1)关于yn+j和 f n+j都是线性组合,所以这一类方法都称为线性k步法. 如果k=0, 则称显式k步法; 如果k0, 则
21、称隐式k步法,定义: 设y(x)是初值问题的准确解,线性多步法在 xn+k 上的局部误差为,若,称(5.1)是p阶方法.p1, 则称(5.1)方法 与初值问题(1.1)是相容的.,把 y(xn+ih) 和 y(xn+ih)作 Taylor展开：,代入(5.1)得,其中,满足,47,q=2,3,若选择j,j,使C0=C1=Cp=0, Cp+10, 则就是一种p阶的线性k步方法, Tn+k= Cp+1hp+1y(p+1)(xn) +O(hp+2 ) 称第一项为局部截断误差的主项, Cp+1称为误差常数.,由相容性定义, p1, 即有: C0=C1=0 得:,所以(5.1)与微分 (1.1)相容的充
22、要 条件是(5.6)成立.,48,当 k= 1 时, 若1 =0 则由(5.6)得: 0=0 = 1, 此时即为 欧拉法,由(5.4)可得 C2=1/2 0 所以方法为1阶精度. 局部误差为:,对 k= 1 时, 若1 0 此时方法为隐式公式. 为了确定 0 ,0, 1 可由C0=C1=C2 =0 解得0=1, 0 = 1= 1/2, 得梯形法,由(5.4)得 C3=-1/12 0 故p=2, 为二阶方法,对 k2 ,的多步法可由(5.4) 确定 j,j ,且由(5.5)可得 局部误差.,49,阿当姆斯显式与隐式公式,形如:,的k步法,称为阿当姆斯(Adams)方法. k =0为显式方法, k
23、 0为隐式方法, 称为阿当姆斯显式与隐式公式. 也称Adams- Bashforth公式和Adams-Monlton公式. 这可由方程(1.1)两端积分 求得.,对(1.1)两端在xn+k-1, xn+k上积分,50,用节点xn+k, xn+k-1, xn, 或 xn+k-1, xn 的F(x)的k次或k-1插值多项式 (x)代替F(x)求积分, 即得线性多步公式.,选择节点xn+3, xn+2, xn+1, xn 作插值节点作函数 f (x, y (x)的 三次插值多项式：,可得,就是四步4阶阿达姆斯(显式)公式,其局部截断误差为,51,也可由(5.4)中 C0=C1=Cp=0 推出. 对比
24、(5.1)与(5.7)式,得:此时0= 1= k-2= 0,k-1=1. 由(5.4)显有C0=0. 令: C1=Ck+1=0 则可得 0 ,1, k (若 k =0,则令C0= C1= =Ck=0来求得0 ,1, k-1 ),例如k=3时, 由C1= C2=C3= C4= 0,根据(5.4)得:,52,若3 =0(显式), 则由前三个方程解得: 0 =5/12 1 =-16/12 2 =23/12 可得 k=3 时的阿当姆斯显式公式,再由(5.4)可得: C4=3/8, 所以(5.8)是三阶方法. 且,若3 0 (隐式), 则由四个方程解得: 0 =1/24 1 =-5/24 2 =19/2
25、4 3 =3/8 可得 k=3 时的阿当姆斯隐式公式,53,(5.9)是四阶方法. 局部误差是,类似可以求得常用的阿当姆斯显式或隐式公式,k p 阿当姆斯显式公式 (P364 表9-5) Cp+1,54,k p 阿当姆斯隐式公式 (P364 表9-6) Cp+1,55,米尔尼(Milne)方法与辛普森(Simpson)方法,考虑另一个k=4的显式公式,其中i 是待定常数,可根据使公式的阶尽可能的高这一条件 来确定其数值.,根据(5.4)得, C0=0 令C1= C2=C3= C4= 0, 得:,可得:0 =0 1 = 8/3 2 =-4/3 3 =8/3,56,57,故所求的公式为,称为米尔尼
26、(Milne)公式, 为4阶方法, 它的局部截断误差为：,米尔尼(Milne)方法也可用积分得:,类似还可以构造二步四阶的辛普森(Simpson)公式：,局部截断误差为,右端积分用辛普森求积公式可得辛普森方法:,汉明方法 辛普森公式是二步法中阶数最高的,但它的稳定性较差. 可考虑另一类三步法公式,58,0,1,2,1 ,2,3为常数,为使公式为四阶的, 若取1,=1,则 可的辛普森公式, 若取1,=0, 利用泰勒展开由(5.4), 令 C0=C1= C2=C3= C4= 0 得,解之得： 0= - 1/8 2,=9/8 1 =-3/8 2 =6/8 3 =3/8,得称为汉明(Hamming)公
27、式,59,60,局部截断误差为：,它是一个四阶隐式三步法.,且可知 C5 =-1/40,线性多步法要用到二个以上的初始值才能开始计算,不能“自启动”.一般应选择同阶或高阶的单步法算出前几个值,再使用线性多步法.由于线性多步法每一步实际上只计算一次函数值,又能获得较高的精度,所以应用比较广泛.,61,例 隐式三步法,解: k=3,0=1=0,2=-1,3=1,解之,0=1/24,1=-5/24,2=19/24,3=9/24.,公式为：,62,63,局部截断误差为,它是一个线性三步四阶隐式公式,称为四阶阿达姆斯内插公式,应用十分普遍.,此外常用的还有线性四步四阶的显式阿达姆斯公式,常称为阿达姆斯外
28、推公式：,64,预估校正法,对于隐式线性多步法, 计算时要进行迭代, 计算量大, 通常 采用显式方法先给出yn+k 的初始近似, 记为y p n+k 称为P(Predictor) 预估, 再计算 fn+k 的函数值E(Evaluation) , 继而以隐式公式计 算yn+k 称为C(Corrector)校正及M(Modifier)改进,阿达姆斯四阶预测-校正公式的PECE模式,65,阿达姆斯公式的PMECME模式,由四阶阿达姆斯公式的截断误差,对于PECE有,二式相减,66,由此得事后误差估计：,于是加上两个改进步(M),最后再求一次函数值(E)供下一次预测用,就形成阿达姆斯法的PMECME模
29、式,67,修正米尔尼-汉明预测-校正法(PMECME),利用米尔尼公式(5.11)进行预测,再用哈明公式(5.15)进行校正,中间利用两公式局部截断误差的主项系数的比例进行改进两次,就得到PMECME模式,此模式称修正米尔尼-汉明预测-校正法,有时侯并不一定要按系数公式(5.4)确定多步法(5.1)的 系数 i , i ,可直接用泰勒展开进行构造.,例 解初值问题的显式二步法,试确定系数i , i使方法阶数尽量的高,并求局部截断误差,解:,68,69,为使方法阶数尽量高,令,解得:,此时为三阶公式, 且所求局部误差为:,例 证明存在的一个值,使线性多步法是四阶的,解: 只需证明局部误差 T n
30、+1= O(h5) 由泰勒展开可得:,70,71,当 = 9 时, T n+1= O(h5). 故方法是四阶的.,72,返回,6 一阶常微分方程组和高阶方程(P373),一阶常微分方程组可写成向量形式,然后所有适用于一阶常微分方程的方法都可以移植到一阶常微分方程组.,考察一阶方程组,若采用记号,则上述问题可表示为,73,1.用欧拉公式为：,例如以含两个方程的常微分方程组为例：,74,75,2.用经典的四阶R-K公式为：,76,3.用阿达姆斯PEC模式计算为：,例如m阶微分方程,初始条件,引进新变量,则(6.4)可为如下一阶方程组,高阶常微分方程,原则上可以化为一阶方程组来求解,初始条件为,77,78,以两阶常微方程为例：,则可令z=y,化为一阶方程组求解：,完全套用前述的方程组各解法即可解决.,返回,
