最新9常微分方程数值解59382

常微分方程数值解法科学技术中的许多问题，在数学中往往归结为微分方程的求解问题。

例如天文学中研究星体运动，空间技术中研究物体飞行等，都需要求解常微分方程初值问题。

最简单的常微分方程初值问题是座=f (x, y), a £ x £ b,dx ' " (*)y(x o) =y°.定理如果f (x, y)在区域D ={( x,y) a《x〈b, y<°°}内连续，且关于y满足Lipschitz条件，那么初值问题 (*)的解存在且惟一。

因为数值解是求微分方程解y(x)的近似值，所以总假定微分方程的解存在惟一，即初值问题(*)中的f(x, y)满足定理的条件。

除特殊情形外，初值问题(*) 一般求不出解析解，即使有的能求出解析解，其函数表示式也比较复杂，计算量比较大，况且实际问题往往只要求在某一时刻解的函数值，因此，有必要研究初值问题(*)的数值解法。

所谓初值问题(*)的数值解法，就是寻求解y(x)在区间［a,b］上的一系列点x〔：：：x^：： x3 ：：：..• ::：x n:：:■■-上的近似值y i,y2,…，y n,….记h j = x —x口(i =1,2, |||)表示相邻两个节点的间距，称为步长。

求微分方程数值解的主要问题：(1)如何将微分方程也=f (x, y)离散化，并建立求其数值解的递推公式；dx(2)递推公式的局部截断误差、数值数y n与精确解y(x n)的误差估计；(3)递推公式的稳定性与收敛性.、Euler方法Euler法是最早的解决一阶常微分方程初值问题的一种数值方法，虽然它不够精确，很少被采用，但是它在某种程度上反映了数值方法的基本思想和特征。

考虑初值问题% f(x, y), (1.1)dxyMi。

，(1.2) 为了求得它在等距离散点x1 <x2<x n上的数值解，首先将(1.1)离散化。

设h =X i —x — (i =1,2,|||)，将式(1.1)离散化的办法有 Taylor 展开法、数值微分法及数值积分法。

i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。

我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。

通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。

本章讨论常微分方程最常用的近似数值解法--差分方法。

先来讨论最简单的Euler 法。

为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。

在微积分课程中我们熟知，微商（即导数）是差商的极限。

反过来，差商就是微商的近似。

在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

下面我们用数值积分法重新导出 Euler 法以及其它几种方法。

为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。

数值分析--第9章常微分⽅程数值解数值分析--第9章常微分⽅程数值解第九章常微分⽅程数值解法许多实际问题的数学模型是微分⽅程或微分⽅程的定解问题。

如物体运动、电路振荡、化学反映及⽣物群体的变化等。

常微分⽅程可分为线性、⾮线性、⾼阶⽅程与⽅程组等类；线性⽅程包含于⾮线性类中，⾼阶⽅程可化为⼀阶⽅程组。

若⽅程组中的所有未知量视作⼀个向量，则⽅程组可写成向量形式的单个⽅程。

因此研究⼀阶微分⽅程的初值问题=≤≤=0)(),(y a y b x a y x f dx dy ， (9-1)的数值解法具有典型性。

常微分⽅程的解能⽤初等函数、特殊函数或它们的级数与积分表达的很少。

⽤解析⽅法只能求出线性常系数等特殊类型的⽅程的解。

对⾮线性⽅程来说，解析⽅法⼀般是⽆能为⼒的，即使某些解具有解析表达式，这个表达式也可能⾮常复杂⽽不便计算。

因此研究微分⽅程的数值解法是⾮常必要的。

只有保证问题(9-1)的解存在唯⼀的前提下，研究其数值解法或者说寻求其数值解才有意义。

由常微分⽅程的理论知，如果(9-1)中的),(y x f 满⾜条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续；（2）),(y x f 在D 上关于y 满⾜Lipschitz 条件，即存在常数L ，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟⼀的连续解)(x y y =。

在下⾯的讨论中，我们总假定⽅程满⾜以上两个条件。

所谓数值解法，就是求问题(9-1)的解)(x y y =在若⼲点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的⽅法。

),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由n x 到1+n x 的步长。

今后如⽆特别说明，我们总假定步长为常量。

第8章常微分方程的数值解法8.4单步法的收敛性与稳定性8.4.1相容性与收敛性上面所介绍的方法都是用离散化的方法，将微分方程初值问题化为差分方程初值问题求解的.这些转化是否合理？即当h →∞时，差分方程是否能无限逼近微分方程，差分方程的解n y 是否能无限逼近微分方程初值问题的准确解()n y x ，这就是相容性与收敛性问题.用单步法(8.3.14)求解初值问题(8.1.1)，即用差分方程初值问题100(,,)()n n n n y y h x y h y x y ϕ+=+⎧⎨=⎩(8.4.1)的解作为问题(8.1.1)的近似解，如果近似是合理的，则应有()()(,(),)0 (0)y x h y x x y x h h hϕ+--→→(8.4.2)其中()y x 为问题(8.1.1)的精确解.因为0()()lim ()(,)h y x h y x y x f x y h→+-'==故由(8.4.2)得lim (,,)(,)h x y h f x y ϕ→=如果增量函数(,(),)x y x h ϕ关于h 连续，则有(,,0)(,)x y f x y ϕ=(8.4.3)定义8.3如果单步法的增量函数(,,)x y h ϕ满足条件(8.4.3)，则称单步法(8.3.14)与初值问题(8.1.1)相容.通常称(8.4.3)为单步法的相容条件.满足相容条件(8.4.3)是可以用单步法求解初值问题(8.1.1)的必要条件.容易验证欧拉法和改进欧拉法均满足相容性条件.一般地，如果单步法有p 阶精度（1p ≥），则其局部截断误差为[]1()()(,(),)()p y x h y x h x y x h O h ϕ++-+=上式两端同除以h ，得()()(,,)()p y x h y x x y h O h hϕ+--=令0h →，如果(,(),)x y x h ϕ连续，则有()(,,0)0y x x y ϕ'-=所以1p ≥的单步法均与问题(8.1.1)相容.由此即得各阶龙格－库塔法与初值问题(8.1.1)相容.定义8.4一种数值方法称为是收敛的，如果对于任意初值0y 及任意固定的(,]x a b ∈，都有lim () ()n h y y x x a nh →==+其中()y x 为初值问题(8.1.1)的精确解.如果我们取消局部化假定，使用某单步法公式，从0x 出发，一步一步地推算到1n x +处的近似值1n y +.若不计各步的舍入误差，而每一步都有局部截断误差，这些局部截断误差的积累就是整体截断误差.定义8.5称111()n n n e y x y +++=-为某数值方法的整体截断误差.其中()y x 为初值问题(8.1.1)的精确解，1n y +为不计舍入误差时用某数值方法从0x 开始，逐步得到的在1n x +处的近似值（不考虑舍入误差的情况下，局部截断误差的积累）.定理8.1设单步法(8.3.14)具有p 阶精度，其增量函数(,,)x y h ϕ关于y 满足利普希茨条件，问题(8.1.1)的初值是精确的，即00()y x y =，则单步法的整体截断误差为111()()p n n n e y x y O h +++=-=证明由已知，(,,)x y h ϕ关于y 满足利普希茨条件，故存在0L >，使得对任意的12,y y 及[,]x a b ∈，00h h <≤，都有1212(,,)(,,)x y h x y h L y y ϕϕ-≤-记1()(,(),)n n n n y y x h x y x h ϕ+=+，因为单步法具有p 阶精度，故存在0M >，使得1111()p n n n R y x y Mh ++++=-≤从而有111111111()()()(,(),)(,,)()(,(),)(,,)n n n n n n n p n n n n n n p n n n n n n e y x y y x y y y Mh y x h x y x h y h x y h Mh y x y h x y x h x y h ϕϕϕϕ+++++++++=-≤-+-≤++--≤+-+-1(1)p nMh hL e +≤++反复递推得11111101110(1)(1)1(1)(1)(1)(1)1(1)p p n n n p n n p n e Mh hL Mh hL e hL hL Mh hL e hL Mh hL e hL+++-+++++⎡⎤≤++++⎣⎦⎡⎤≤+++++++⎣⎦+-≤++因为00()y x y =，即00e =，又(1)n h b a +≤-，于是ln(1)1()(1)(1)b a b a hL n L b a h h hL hL e e --++-+≤+=≤所以()11()p L b a p n M e h e O h L -+⎡⎤≤-=⎣⎦推论设单步法具有p (1p ≥)阶精度，增量函数(,,)x y h ϕ在区域G ：, , 0a x b y h h ≤≤-∞<<+∞≤≤上连续，且关于y 满足利普希茨条件，则单步法是收敛的.当(,)f x y 在区域:,D a x b y ≤≤-∞<<+∞上连续，且关于y 满足利普希茨条件时，改进欧拉法，各阶龙格－库塔法的增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，因而它们都是收敛的.关于单步法收敛的一般结果是：定理8.2设增量函数(,,)x y h ϕ在区域G 上连续，且关于y 满足利普希茨条件，则单步法收敛的充分必要条件是相容性条件(8.4.3).8.4.2稳定性稳定性与收敛性是两个不同的概念，收敛性是在假定每一步计算都准确的前提下，讨论当步长0h →时，方法的整体截断误差是否趋于零的问题.而稳定性则是讨论舍入误差的积累能否对计算结果有严重影响的问题.定义8.6若一种数值方法在节点值n y 上有一个大小为δ的扰动，于以后各节点()m y m n >上产生的偏差均不超过δ，则称该方法是稳定的.我们以欧拉法为例进行讨论.假设由于舍入误差，实际得到的不是n y 而是n n n y y δ=+，其中n δ是误差.由此再计算一步，得到1(,)n n n n y y hf x y +=+把它与不考虑舍入误差的欧拉公式相减，并记111n n n y y δ+++=-，就有[]1(,)(,)1(,)n n n n n n y n nh f x y f x y hf x δδηδ+⎡⎤=+-=+⎣⎦其中y f f y∂=∂.如果满足条件1(,)1y n hf x η+≤，(8.4.4)则从n y 到1n y +的计算，误差是不增的，可以认为计算是稳定的.如果条件(8.4.4)不满足，则每步误差将增大.当0y f >时，显然条件(8.4.4)不可能满足，我们认为问题本身具有先天的不稳定性.当0y f <时，为了满足稳定性要求(8.4.4)，有时h 要很小.一般的，稳定性与方法有关，也与步长h 的大小有关，当然也与方程中的(,)f x y 有关.为简单起见，通常只考虑数值方法用于求解模型方程的稳定性，模型方程为y y λ'=(8.4.5)其中λ为复数.一般的方程可以通过局部线性化转化为模型方程，例如在(,)x y 的邻域(,)(,)(,)()(,)()x y y f x y f x y f x y x x f x y y y '==+-+-+略去高阶项，再作变量替换就得到u u λ'=的形式.对于模型方程(8.4.5)，若Re 0λ>，类似以上分析，可以认为方程是不稳定的.所以我们只考虑Re 0λ<的情形，这时不同的数值方法可能是数值稳定的或者是数值不稳定的.当一个单步法用于试验方程y y λ'=，从n y 计算一步得到1()n n y E h y λ+=(8.4.6)其中()E h λ依赖于所选的方法.因为通过点(,)n n x y 试验方程的解曲线（它满足,()n n y y y x y λ'==）为[]exp ()n n y y x x λ=-，而一个p 阶单步法的局部截断误差在()n n y x y =时有1111()()p n n n T y x y O h ++++=-=，所以有1exp()()()p n n y h E h y O h λλ+-=(8.4.7)这样可以看出()E h λ是h e λ的一个近似值.由(8.4.6)可以看到，若n y 计算中有误差ε，则计算1n y +时将产生误差()E h λε，所以有下面定义.定义8.7如果(8.4.6)式中，()1E h λ<，则称单步法(8.3.14)是绝对稳定的.在复平面上复变量h λ满足()1E h λ<的区域，称为方法(8.3.14)的绝对稳定区域，它与实轴的交称为绝对稳定区间.在上述定义中，规定严格不等式成立，是为了和线性多步法的绝对稳定性定义一致.事实上，()1E h λ=时也可以认为误差不增长.(1)欧拉法的稳定性欧拉法用于模型方程(8.4.5)，得1(1)n n y h y λ+=+，所以有()1E h h λλ=+.所以绝对稳定条件是11h λ+<，它的绝对稳定区域是h λ复平面上以(1,0)-为中心的单位圆，见图8.3.而λ为实数时，绝对稳定区间是(2,0)-.Im()h λRe()h λ2-1-O 图8.3欧拉法的绝对稳定区域(2)梯形公式的稳定性对模型方程，梯形公式的具体表达式为11()2n n n n h y y y y λλ++=++，即11212n nh y y h λλ++=-，所以梯形公式的绝对稳定区域为12112h h λλ+<-.化简得Re()0h λ<，因此梯形公式的绝对稳定区域为h λ平面的左半平面，见图8.4.特别地，当λ为负实数时，对任意的0h >，梯形公式都是稳定的.Im()h λRe()h λO 图8.4梯形公式的绝对稳定区域(3)龙格－库塔法的稳定性与前面的讨论相仿，将龙格－库塔法用于模型方程(8.4.5)，可得二、三、四阶龙格－库塔法的绝对稳定区域分别为211()12h h λλ++<23111()()126h h h λλλ+++<2341111()()()12624h h h h λλλλ++++<当λ为实数时，二、三、四阶显式龙格－库塔法的绝对稳定区域分别为20h λ-<<、2.510h λ-<<、 2.780h λ-<<.例8.5设有初值问题21010101(0)0xy y x x y ⎧'=-≤≤⎪+⎨⎪=⎩用四阶经典龙格－库塔公式求解时，从绝对稳定性考虑，对步长h 有何限制？解对于所给的微分方程有2100,(010)1f x x y xλ∂==-<≤≤∂+在区间[0,10]上，有201010max ||max51t x x λ<<==+由于四阶经典龙格－库塔公式的绝对稳定区间为 2.7850h λ-<<，则步长h 应满足00.557h <<.。