常微分方程数值解法

1、i.常微分方程初值问题数值解法本章讨论常微分方程初值问题数值解法，主要是差分法。解微分方程的所谓差分法的要点如下：首先是区域的离散，即将连续的求解区域离散化成有限个网格点。其次是方程的离散， 例如用差商代替微商，或者对微分方程积分使之变成积分方程，然后数值积分，或者最后得到网格点上的近似解所满足的一个差分方程，解之即得差分解。i.1常微分方程差分法考虑常微分方程初值问题：求函数 u(t)满足竺=f (t,u), 0 : t 汀(i.1a)dtu(0) =u0(i.1b)其中f(t,u)是定义在区域G : 0兰t兰T, u £珀上的函数，u0和T是给定的常数。我们假 设f(t,u)对u

2、满足Lipschitz条件，即存在常数 L使得f (t,uj f(t,u2)兰 L 5 -上,g0,T;比,比亡(亠严)(i.2)这一条件保证了( i.1)的解是适定的，即存在，唯一，而且连续依赖于初值u0。通常情况下，(i.1)的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨 论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler法。为此，首先将求解区域0,T 离散化为若干个离散点：0 =t° : ti ： HI ： t”:：tN =T(i.3)其中tn =hn , h0称为步长。在微积分课程中我们熟知，微商(即导数)是差商的极限。反过来，差商就是微商的近

3、似。在t =t0处，在(i.1a)中用向前差商u(tl)曲0)代替微商du，便得hdtu(tj =u(t°) hf(t°,u(t。)p如果忽略误差项；0，再换个记号，用 口代替u(ti)便得到q -u° =hf(t0，u。)般地，我们有Euler 方法：un4t =Un +hf (tn,Un), n = 0,1,1", N -1 (i.4)从(i.1b)给出的初始值u0出发，由上式可以依次算出tJiltN上的差分解u1| ,uN。F面我们用数值积分法重新导出Euler法以及其它几种方法。为此，在区间tn ,tn 1上积分常微分方程（i.1a），得tn J

4、LU（tn=U（tn）f（t,U（t）dt（ i.5）用各种数值积分公式计算（i.5 ）中的积分，便导致各种不同的差分法。例如，若用左矩形公 式就得到Euler法（i.4）。如果用右矩形公式，便得到下面的：隐式 Euler 方法：Un4i=Un+h f（tn*,Un）, n= 0,1，卅，N1 （ i.6）类似地，如果用梯形公式，就得到改进的Euler方法h2f(tn1，Un1)f 心 1，®,n=0,1川，N -1(i.7)当f（t,U）关于U是非线性函数的时候，不能由（i.6）和（i.7）从Un直接算出Un1，称这一类方法为隐式，通常采用某种迭代法求解。例如，将一般的隐式方法写成

5、Un 1=F（tn,Un,Un J （i.8）则可以利用如下的迭代法由un算出Un 4 :u： = F(tn，Un，Uk + )，k =0,1川0Un+=Un(i.9)关于k的迭代通常只需进行很少几步就可以满足精度要求了。为了避免对隐式方法进行迭代的麻烦，比如说对于改进的Euler方法（i.7），可以采用某种预估法近似算出f （tn 1,un，然后再用（i .7）作校正，这就导致所谓预估校正法。下面给出一个例子:预估:| un = f(tn,un)Un 1 = Unj 2hU nU n 1 _ f(tn-1,Un 1)(i.10)校正:Un1 f 2(Un1 Un)这是一个多步法，即计算节点t

6、n 1上的近似值Un 4时，除了用到前一点的近似值Un之外，还要用到un4，甚至可能用到un/,HI。而用前面的各种 Euler法计算节点tn1上的近似值 比计时，只用到un,因此称之为单步法。F面给出另一个多步法的例子。在区间出心七上积分（i.1a），得U（tn .2） =U（tn）J f （t,U（t）dtLn用Simps on公式（即把被积函数看作二次函数）近似计算积分，便得到Mil ne 法:(i.11)un2 二 un h(fn 4fn 1 fn 2), n =0,1,”|,N -2 3用多步法（i.10 ）或（i.11）计算时，必须先用某种单步法由u0计算出u1，称为造表头。然后再

7、逐次算出U2,U3,|,UN。一般说来，多步法比Euler法等简单的单步法精度要高一些。下面我们讨论一类所谓Ronge-Kutta法。他们是单步法，但是其精度可以与多步法比美。最常用的是下面的标准Ronge-Kutta 法和 Gear 法：Ki = f (tn，Un)hKi、K2 二 f (tn * Nn 2 )标准 Lun ge-Kutta 法:h 2+nc-3K(i.12)9 / 7K4 = f (tnh,UK3)hUn i °(Ki 2& 2K3 K4). 6K1 = f (tn, un )K2 = f (tn 十2，山 +4Gea方法:(i.13)彳 & =

8、f (tn +£,叫 +( ft 岛心 *(1-£)K2)K4 = f (J 卄心 +(-£)心 +(1+£)心) 叫出=山+£(心+2(1岛心+2(1 + £)厶+ ©)从几何上，Ronge-Kutta法可以粗略地解释为：在区间tn,tn1中选取若干个点（可以重复）.1二tn2 -川-k二tn 1，仅仅利用在区间 出乙内可以得到的所有信息，依 次给出函数f （t,u（t）在这些点上尽可能精确的的近似值K1,K2l,Kk，然后把它们组合起来，尽可能精确地近似计算（i.5）中的积分。Ronge-Kutta法一般的构造方式如下。

9、选定常数k ,令K f（tn,un）K2 二 f（tn-hU亦1）K3 = f (tn，二込叫 * -31K1 -32 K2)Kk =f(tn ：kh,Un *Ki 七K2 III ",kKk)Un 1 =6 h( "Ki 上心川K)(i.14)选取这些待定常数:i，'j-m的原则是：将(i.14)在(tn,Un)作Taylor展开，然后按照h的 幂重新整理，使得1213Un 1二比 /2h3h ' ( "5)2!31与微分方程(i.1a)的解u =u(t)在t =tn处的Taylor展式1 2 1 3U(tn 1) =U(tn) 帥可时川(6)有

10、尽可能多的项重合，即要求1 二 fn,2 二仁 3 二 fn ,川这里花二f(tn,Un),仁df (t,u)dt等等。(t,U)虫“叫)(i.18)(i.17)按照(i.14)构造出的都是显式Runge-Kutta方法,在每一个 &的表达式中只出现 Un。如果允许在某一个 Kj的表达式中出现un彳，则可以导出隐式 Runge-Kutta方法。i.2常微分方程组与高阶常微分方程先来考虑下面的常微分方程组初值问题dJ=f1(t,U1,U2l|,Um)dt如=f 2(t,U1,U2,川，Um)dtlllllldJ号=fm(t,U1,U2,川,Um)dtu1(0u；,um(0)=u：利用向量

11、记号，上式可以改写为u (t) =f (t,u).u(0) =Uo上节中各方法都可以直接应用到常微分方程组(i.18)。例如，Euler方法成为UU nhf (tn, Un), n =0,1,川,N -1再来考虑高阶常微分方程dmududmJu(i.19)廿二(t,u,不,川,耐)，t 0 d mJudu.mA (0) = Vm , ： 11 7 (0) = V2 , U (0) = vidtdt这时，可以令/ 1 2m、 /u = (u ,u , ,u ) =(u,dudtdmJu，川討(i.20)f = (f1, f2, II, fm) =(u2,u3, |lLum, f(t,u1,u2J

12、H,um) (i.21)U0 =(V1,V2,|),Vm)(i.22)于是可以把高阶常微分方程(i.19)化成一阶常微分方程组(i.18)。i.3收敛性与稳定性截断误差粗略地说,截断误差可以定义为将微分方程解带入到差分方程后得到的误差，代表了微分方程与差分方程之间的误差。例如，由Taylor展式和微分方程(i.1a)得到2!2! dt其中n是区间tn,tn 1上某个常数。与Euler法ununhf(tn,un)相比较，定义余项df (t,u(t) 为Euler法的截断误差，它关于h是2阶的，记为0(f)。将上节中讨论2! dt t的单步法写成一般形式F(un,un1,h) =0，则可以定义截断

13、误差为F(u(tn),u(tn1),h)。对1于每一个差分法在适当的点(例如 t =tn，tn-1或者(tn tn 1)作Taylor展开，就可以得2到截断误差的阶。对多步法可以类似处理。上节中各方法的截断误差阶分别为：表i.1常微分方程差分法的截断误差阶差分法Euler隐式Euler改进Euler预估校正(i.10)Milne标准 Ronge-KuttaGear截断误差o(h2)o(h2)o(h3)o(h3)0(h5)0(h5)o(h5)相容性 一个差分方法称为相容的，如果其截断误差至少是一阶的。收敛性实用中我们更关心的是近似解的收敛性，即卩h > 0时，是否有un > u(tn

14、)。在适当的条件下，例如步长h足够小，右端函数f和解u足够光滑等等，可以证明以上讨论的各方法都是收敛的，并且有估计式U(tn) -U<Chp(i.23)其中常数C与u和f有关，与h和n无关；阶数p等于截断误差的阶数减去1。稳定性收敛性考虑的是差分方程的精确解与微分方程的精确解之间的误差。但是在计算机上求差分解时，由于计算机字长的限制，不可避免地产生舍入误差。另外，方程中的系数和 初值等等由于测量条件等限制，也会产生数据误差。而稳定性研究的就是舍入误差和数据误差对差分方程计算结果的影响，即计算过程中某一步的 “差之毫厘”会不会导致后面结果的“失之千里”。根据误差来源和误差衡量标准的不同，可

15、以定义形形色色的稳定性。一般说来，隐式方法的稳定性好于同类型的显式方法（例如Euler方法与隐式Euler方法）。另外，为了保证稳定性，常常要求步长h足够小。对于常微分方程数值解来说，最简单常用的是关于初值的稳定性。 以单步法为例。称某差分 法关于初值稳定，如果对同一个微分方程和所有足够小的h，存在常数C，使得从不同初值u0和V。出发的两个差分解 CuJ和：vn /之间的误差满足maxn Un v. EC u。一v°(i.24)其中N二T/h。由于任意一对Vn都可以看作是初值，关于初值的稳定性其实是考察如下问题：假设某一步计算有误差， 而其后的计算不再有误差, 的影响如何。本章所讨论

16、的所有差分方法都是关于初值稳定的。 时，通常指的是关于初值稳定。不难证明，如果一个差分方法是稳定且相容的，则一定收敛，即那么这一步的误差对以后结果不加区别地说差分方法稳定口代）一叫T 0,当nT呛(i.25)du dt可以认为这个模型方程的真解和差分解的性质代表了非线性常微分方程dudtf (t,u)在某因而，我们进一步考察差分方法的所在差分解的实际计算中， 每步都难以避免地舍入误差。谓绝对稳定性，即每一步的舍入误差积累起来，是否会对后面的运算结果产生太坏的影响。考察某一个差分法对非线性常微分方程的绝对稳定性是十分困难的事情，无法得到象对初值稳定性那样的漂亮结论。通常的做法是考虑模型方程(i.

17、26)一个局部的真解和差分解的性质。（例如在（t, u）二（tn,un）处将f （t, u）线性化，将-看作f（tn,un），忽略其他项。）记h二。关于某一个差分法的绝对稳定性的典型结果是：给:u出一个区间（，J，使得对于所有的 h （：），求解（i.25）的这个差分法都是绝对稳定的。这个区间（）称之为这个差分法的绝对稳定域。绝对稳定域越大，h的取值范围就越大，的允许区值就越多，所代表的非线性常微分方程也就越多。表i.2常微分方程差分法的绝对稳定域差分法Euler隐式Euler改进EulerMilne标准R-KGear绝对稳定域（-2，0）(皿,0)(£0)(0,0)(-2.78,0)(-2.78,0)关于绝对稳定性的一个必要条件是：若一个多步法是相容并且稳定的，则它绝对稳定（即 它的绝对稳定域非空）的一个必要条件是：0。这相当于仅当 f : 0时，多步法才可能绝对稳定。而当 f 0时，在差分计算中误差增长很快。那么这时候是否就无法用差分法du近似计算微分方程了呢？