常微分方程初值问题.ppt

第五章 常微分方程初值问题n n引言引言引言引言n n基本概念基本概念基本概念基本概念n nEulerEuler方法及其改进方法及其改进方法及其改进方法及其改进1§1 引言引言一、一、 常微分方程的定解问题与应用常微分方程的定解问题与应用应用：自应用：自然科学领域，如物理；然科学领域，如物理；工程技术问题，如石油勘探工程技术问题，如石油勘探 常微分方程的定解问题主要有常微分方程的定解问题主要有初值问题初值问题和和边值问题边值问题两两大类，我们仅考虑大类，我们仅考虑初值问题初值问题例：例： 马尔萨斯人口模型：假设某特定区域在马尔萨斯人口模型：假设某特定区域在t0时刻的人口时刻的人口p(t0) = p0为已知的，该区域人口的自然增长率为为已知的，该区域人口的自然增长率为α人口的增长与人口的总数成正比，所以口的增长与人口的总数成正比，所以 t 时刻的人口总时刻的人口总数数 p(t)满足如下的微分方程：满足如下的微分方程：一般地，称这样的方程为一般地，称这样的方程为模型方程模型方程2二、二、 常微分方程的解法：常微分方程的解法：①①解析法解析法：给出精确解析解只适合少数简单情况给出精确解析解。

只适合少数简单情况②②近似解法近似解法：给出解的近似表达式如级数法：给出解的近似表达式如级数法,逐步逼近法逐步逼近法③③数值方法数值方法：给出方程在离散点上的近似解它适合计算机：给出方程在离散点上的近似解它适合计算机求解求解,应用广泛应用广泛,具有理论应用价值具有理论应用价值三、常微分方程初值问题的数值方法三、常微分方程初值问题的数值方法单步法单步法Euler方法方法Taylor方法和方法和Runge-Kutta方法方法多步法多步法Adams方法和一般线性多部法方法和一般线性多部法线性多部法的收敛性与稳定性线性多部法的收敛性与稳定性3§2 基本概念基本概念一、常微分方程初值问题的一般提法一、常微分方程初值问题的一般提法问题问题:求函数求函数满足满足其中其中: f (x, y) 为已知函数为已知函数, α是已知值是已知值. (可能是观察值或实验值可能是观察值或实验值)(2)(1)基本条件基本条件:设设①① f (x,y)在在D上连续；上连续；②② f (x,y)在在D上关于变量上关于变量y满足满足Lipschitz连续条件：连续条件：其中其中L为为Lipschitz常数。

常数3)4若若f (x,y)在在D上满足基本条件上满足基本条件,一阶常微分方程初值问一阶常微分方程初值问题题(1),(2)对任意给定的对任意给定的α存在唯一解存在唯一解且在且在[a,b]上连续上连续可微可微.定理定理1关于解关于解y(x)的适定性：的适定性：定义定义1 方程方程(1),(2)的解的解y(x)称为称为适定的适定的,若存在常数若存在常数ε>0，，K > 0,对任意满足条件对任意满足条件常微分方程初值问题：常微分方程初值问题：存在唯一解存在唯一解z(x),且有且有摄动摄动(扰动扰动)误差误差(4)(1),(2) 上各加一上各加一个个摄动摄动(扰动扰动)项项.(1),(2) 上各加一上各加一个个摄动摄动(扰动扰动)项项.5定理定理2①①适定问题的解适定问题的解y(x)连续依赖于连续依赖于(1)式右端的式右端的f(x,y)和初值和初值.或者说解或者说解y(x)关于关于(1)式右端的式右端的f (x, y)和初值稳定和初值稳定. ②②假设假设f (x, y)在在D上满足基本条件上满足基本条件,从而方程从而方程(1),(2)的解的解y(x)存在且适定存在且适定.注注:若若f (x, y)在在D上满足基本条件上满足基本条件,则微分方程则微分方程(1),(2)的的解解y(x)是是适定适定的的.6假设初值问题假设初值问题(1)的解的解y = y(x)唯一存在且足够光滑唯一存在且足够光滑, 对求解区域对求解区域[a,b]做剖分做剖分 1 构造数值解法的基本思想构造数值解法的基本思想在区间在区间[xk, xk+1]上对方程上对方程(1)做积分做积分,则有则有二、二、 初值问题数值解的基本概念初值问题数值解的基本概念 常用等步长常用等步长:, 则有则有(1)，，(2)的准确解记为的准确解记为y(x),称为称为步长步长。

的近似解记为的近似解记为7因此，建立节点处近似值因此，建立节点处近似值yn满足的差分公式满足的差分公式称为称为Euler公式公式. 称为称为梯形公式梯形公式. 若对若对*式右边的积分应用梯形求积公式，则可导出差分公式式右边的积分应用梯形求积公式，则可导出差分公式对对右边的积分应用左矩形公式，则有右边的积分应用左矩形公式，则有8称为称为Euler中点公式中点公式或称或称双步双步Euler公式公式. 若在区间若在区间[xk-1, xk+1]上对方程上对方程(1)做积分做积分,则有则有对对右边的积分应用中矩形求积公式，则得差分公式右边的积分应用中矩形求积公式，则得差分公式因此，求初值问题数值解的基本方法是因此，求初值问题数值解的基本方法是步进法步进法.即逐个节点即逐个节点计算计算, 由由 yi ( i k ) 计算计算yk+1步进法步进法单步法单步法:多步法多步法: 由由yk-j ( j=0,1,…,l-1 )计算计算yk+1共用到共用到l个个值值.即即yk+1，，yk-1,…,yk-l+1称为称为l步法步法仅由仅由yk计算计算yk+19单步法与多步法的区别单步法与多步法的区别:(1)计算方面计算方面: l步方法只用于步方法只用于的计算的计算,即即y0，，y1,…,yl-1的计算要用其它方法。

的计算要用其它方法2)理论分析理论分析: 单步法比单步法比 l > 1的多步法容易分析的多步法容易分析(稳定性稳定性).(3)选步长方面选步长方面:单步法容易改变步长单步法容易改变步长.(4)精度精度:多步法精度高一些多步法精度高一些.单步法与多步法均有单步法与多步法均有显式显式和和隐式方法隐式方法之分之分.在在Euler公式和公式和Euler中点公式中中点公式中,需计算需计算yk+1的已显式表示的已显式表示,称为称为显式公式显式公式,而而梯形公式中梯形公式中,需要计算的需要计算的yk+1隐含在等式两侧隐含在等式两侧,称其为称其为隐式公式隐式公式.计算公式依次可抽象写成计算公式依次可抽象写成:显式单步法显式单步法:隐式单步法隐式单步法:(5)(6)10该式右端项含有该式右端项含有因此若求因此若求需要解方程需要解方程注注:显式多步法显式多步法:隐式多步法隐式多步法:线性多步法线性多步法:注注: (9)关于关于yk-i, fk-i都是线性的都是线性的其中其中是独立于是独立于k和和f 的常数9)是显式是显式(右端不含有右端不含有yk+1)；；则则(9)是隐式的是隐式的7)(8)(9)三、微分方程初值问题的数值解法讨论的问题三、微分方程初值问题的数值解法讨论的问题:1.方法构造方法构造2.误差分析误差分析3.稳定性稳定性11§3 Euler方法方法考虑考虑一、显式一、显式Euler方法方法(折线法折线法)1.显式欧拉公式显式欧拉公式其中其中fk = f (xk, xk+1).(1)(2)(10)设节点为设节点为:a = x0

式14(3) 左矩数值积分法左矩数值积分法将将(1)式两端从式两端从xk到到xk+1积分，得积分，得数值积分采用左矩形公式数值积分采用左矩形公式 , 即即由初始条件亦得由初始条件亦得(10)式式.o类似的类似的152. 几何意义几何意义o 方程方程(1),(2)的解曲线的解曲线 y(x) 过点过点 p0(x0, y0), 具有斜率具有斜率 f0从从 p0出发以出发以 f0 为斜率做直线段，交为斜率做直线段，交 x = x1于点于点 p1(x1, y1), 此此时时 y1 = y0+h0f0, 过过p1(x1,y1)的解曲线如图示，具有斜率的解曲线如图示，具有斜率f1再从再从p1出发，以出发，以f1为斜率做直线段，交为斜率做直线段，交x = x2于点于点p2(x2, y2), 此时此时y2 = y1+h1f1, 过过 p2(x2, y2) 的解曲线如图的解曲线如图 以此类推，最终到达终点以此类推，最终到达终点 pn(xn, yn)，这样得到了一条折线，这样得到了一条折线16 它在点它在点pk的右侧具有斜率的右侧具有斜率fk, 与与(1)过过pk的解曲线相切。

因此，的解曲线相切因此，Euler法的几何意义就是用折线法的几何意义就是用折线 作为解曲线的近作为解曲线的近似折线折线称为称为欧拉折线欧拉折线，所以欧拉方法又称为，所以欧拉方法又称为折线法折线法特点：特点：简单，精度低简单，精度低.17二、隐式二、隐式Euler方法和梯形方法方法和梯形方法(Euler方法的改进方法的改进)忽略忽略项项,用用分别近似分别近似可得可得隐式欧拉方法隐式欧拉方法:1. 隐式隐式 Euler方法方法或或(13)将将 y (xk)在在x = xk+1点进行点进行Taylor展开展开18(2)隐式隐式 Euler方法方法(13)是关于是关于yk+1的方程，若求的方程，若求yk需需要解方程要解方程(13)说明说明: (1)隐式隐式 Euler方法也可用向后差商方法也可用向后差商,即用即用近似微分近似微分或者用右矩数值求积公式来建立或者用右矩数值求积公式来建立.192.梯形方法梯形方法取平均取平均,得得忽略忽略项项,用用得梯形方法得梯形方法:与与说明说明: 梯形方法也是隐式方法梯形方法也是隐式方法,要求要求,需解方程需解方程(14).由由分别近似分别近似(14)也可由梯形求也可由梯形求积公式导出积公式导出203.Euler方法、隐式方法、隐式Euler方法、梯形方法等单步法的对应关系方法、梯形方法等单步法的对应关系显式单步法显式单步法显式显式 Euler方法方法隐式单步法隐式单步法隐式隐式Euler方法方法梯形方法梯形方法(隐式隐式)21在实际计算中希望有较好的在实际计算中希望有较好的 , 用较少的迭代步用较少的迭代步(次数次数), 取取得得称称 为为 的的m次迭代次迭代最常用的方法之一是先用显式最常用的方法之一是先用显式Euler方法所得的方法所得的 为为 量较大量较大 ,往往取往往取 作为作为 来用来用,更精确的更精确的 (足够精度的足够精度的 ). 而在实际计算中而在实际计算中 的计算的计算改进改进。

初始值初始值), 再由梯形方法改进一次，即是再由梯形方法改进一次，即是预估－校正预估－校正Euler方法方法(或称为或称为改进改进Euler方法方法)22三、预估－校正三、预估－校正Euler方法方法(显式Euler公式)(梯形公式) 或写成 (15)(16)此算法是此算法是单步递推格式单步递推格式，比隐式公式的迭代求解过，比隐式公式的迭代求解过程简单稳定性高于显式欧拉法稳定性高于显式欧拉法注：注：23分别用显式分别用显式Euler方法，梯形方法和预估－校正方法，梯形方法和预估－校正Euler方方法解初值问题法解初值问题解：解：取取 h =0.1，，Euler方法为方法为:例例解析解：解析解：y = e-x + x24梯形方法为：梯形方法为：预估－校正预估－校正Euler方法：方法：25 EulerEuler方法方法方法方法 梯形方法梯形方法梯形方法梯形方法预估－校正方法预估－校正方法预估－校正方法预估－校正方法 真值真值真值真值0.00.01.0000001.0000000.00.01.0000001.0000000.00.01.0000001.0000000.00.01 10.10.11.0000001.0000004.8×104.8×10－－－－3 31.0047621.0047627.5×107.5×10－－－－5 51.0050001.0050001.6×101.6×10－－－－4 41.0048371.0048370.20.21.0100001.0100008.7×108.7×10－－－－3 31.0185941.0185941.4×101.4×10－－－－4 41.0190251.0190252.9×102.9×10－－－－4 41.0187311.0187310.30.31.0290001.0290001.2×101.2×10－－－－2 21.0406331.0406331.9×101.9×10－－－－4 41.0412181.0412184.0×104.0×10－－－－4 41.0408181.0408180.40.41.0561001.0561001.4×101.4×10－－－－2 21.0700961.0700962.2×102.2×10－－－－4 41.0708001.0708004.8×104.8×10－－－－4 41.0703201.0703200.50.51.0904901.0904901.6×101.6×10－－－－2 21.1062781.1062782.5×102.5×10－－－－4 41.1070761.1070765.5×105.5×10－－－－4 41.1065311.1065310.60.61.1314411.1314411.7×101.7×10－－－－2 21.1485371.1485372.7×102.7×10－－－－4 41.1494041.1494045.9×105.9×10－－－－4 41.1488121.1488120.70.71.1782971.1782971.8×101.8×10－－－－2 21.1962951.1962952.9×102.9×10－－－－4 41.1972101.1972106.2×106.2×10－－－－4 41.1965851.1965850.80.81.2304671.2304671.9×101.9×10－－－－2 21.2490191.2490193.0×103.0×10－－－－4 41.2499751.2499756.5×106.5×10－－－－4 41.2493291.2493290.90.91.2874201.2874201.9×101.9×10－－－－2 21.3062641.3062643.1×103.1×10－－－－4 41.3072281.3072286.6×106.6×10－－－－4 41.3065701.3065701.01.01.3486781.3486781.9×101.9×10－－－－2 21.3675731.3675733.1×103.1×10－－－－4 41.3685141.3685146.6×106.6×10－－－－4 41.3678791.367879数值例子表明：梯形方法和预估－校正数值例子表明：梯形方法和预估－校正Euler方法比显式方法比显式Euler方法有更好的精度。

方法有更好的精度2627四、单步法的局部截断误差、整体截断误差四、单步法的局部截断误差、整体截断误差问题：问题：满足基本条件满足基本条件:①① f (x,y)在在D上连续；上连续； ②② f(x,y)在在D上关于变量上关于变量y 满足满足Lipschitz条件条件.281.局部截断误差的定义局部截断误差的定义设单步法为（数值方法公式）：设单步法为（数值方法公式）：定义定义2 设设y(x)是方程是方程(1),(2)的准确解，称的准确解，称为单步法为单步法(17)在在xk+1点的点的局部截断误差局部截断误差(方法误差方法误差)17)(18)29定义定义3若对充分小的若对充分小的h > 0 , 常数常数c 独立于独立于h(与与h无关无关),称称(17)是是 p 阶方法，即阶方法，即O(hp+1)判断某种方法的阶数往往通过局部截断误差的阶数来确判断某种方法的阶数往往通过局部截断误差的阶数来确定，而局部截断误差的阶容易由公式来确定定，而局部截断误差的阶容易由公式来确定说明：说明：(19)设设y(x)是方程是方程(1)(2)的准确解，的准确解，yk, k=0,1,…,n是单步是单步法法(17)的数值解，称的数值解，称 ek = y(xk) - yk,是单步法是单步法(18)在在xk点的点的整体截断误差整体截断误差。

成立成立定义定义4302.局部截断误差与整体截断误差之间的关系局部截断误差与整体截断误差之间的关系若单步法若单步法(17)的局部截断误差是的局部截断误差是p+1阶的，即阶的，即c1独立于独立于h(不依赖于不依赖于h或与或与h无关无关)，，且函数且函数Φ(x, u, v, h)在区域在区域定理定理4上关于上关于u,v 满足条件满足条件则单步法则单步法(17)是是 p 阶阶方法20)31推论推论: 当当f 在在D上满足基本条件时上满足基本条件时,单步法单步法(17)的阶由的阶由局部截断局部截断误差的阶误差的阶来确定 结论：结论：一般情况下，若局部截断误差是一般情况下，若局部截断误差是 p+1阶的，则单步法是阶的，则单步法是p 阶方法说明：说明：可用可用Taylor展开法估计展开法估计的阶，即的阶，即 的主项的主项高阶项高阶项则对应的单步法是则对应的单步法是p 阶方法由于由于(21)323. 定理的应用（讨论各种方法的阶数）定理的应用（讨论各种方法的阶数） 显式显式Euler方法是一阶方法方法是一阶方法事实上事实上,当当(1),(2)的解的解y(x)二阶连续可导时二阶连续可导时, f (x, y)关于关于y满足满足Lipschitz条件，条件，Lipschitz常数为常数为L，其局部截断误差，其局部截断误差为为:将将y(xk+1)在在xk点展开可得点展开可得局部截断误差有界，存在局部截断误差有界，存在M，使得，使得记记于是于是33 反复递推得：反复递推得：反复递推得：反复递推得：34显式显式显式显式EulerEuler方法的整体截断误差与方法的整体截断误差与方法的整体截断误差与方法的整体截断误差与h h同阶。

同阶35n 隐式隐式Euler方法是一阶方法方法是一阶方法事实上，局部截断误差为事实上，局部截断误差为:将将 f (x, y)在在y用微分中值定理，有用微分中值定理，有则则将方程的解作将方程的解作Taylor展开展开36因此因此37n 梯形方法是二阶方法梯形方法是二阶方法将将 f (xk+1, yk+1)在在xkTaylor展开，有展开，有则则将方程的解作将方程的解作Taylor展开展开将将 f (x, y)在在y用微分中值定理，有用微分中值定理，有38因此因此39n 预估－校正预估－校正Euler方法是二阶方法方法是二阶方法事实上，局部截断误差为：事实上，局部截断误差为：40二阶方法二阶方法即得整体截断误差即得整体截断误差:结论：结论：一般情况下，若局部截断误差是一般情况下，若局部截断误差是 p+1阶的，则单步法阶的，则单步法是是p 阶方法方方 法法  显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大41五、五、 单步法解的唯一性、收敛性和稳定性单步法解的唯一性、收敛性和稳定性u 唯一性唯一性若若f (x, y)在在D上满足基本条件，上满足基本条件，f (x, y)关于关于y的的Lipschitz常数常数为为L时，只要时，只要hkL<1，则由定理，则由定理1可知可知(13)确定了唯一的确定了唯一的yk+1。

同理，只要同理，只要hkL<2，，(14)确定了唯一的确定了唯一的yk+1.u 收敛性收敛性(以以(14)式为例说明式为例说明)当当时时,以以y 为变量的函数为变量的函数:同样，从任意初始值出发同样，从任意初始值出发,迭代迭代(13)、、 (14)，都收敛，都收敛到到 事实上，事实上，42在在 上关于上关于y 满足满足Lipschitz条件条件,且且L－－常数为常数为 ,由由压缩不动点定理压缩不动点定理得方程得方程有唯一不动点有唯一不动点 , 而且从而且从 出发出发,迭代迭代(15)都收敛到都收敛到见非线性方程数见非线性方程数值解的定理值解的定理343u 稳定性稳定性例：例：考察初值问题考察初值问题 在区间在区间[0, 0.5]上的解分别用欧拉显、隐式格式和改进的欧拉格式计算数值解分别用欧拉显、隐式格式和改进的欧拉格式计算数值解0.00.00.10.10.20.20.30.30.40.40.50.5精确解精确解精确解精确解改进欧拉法改进欧拉法改进欧拉法改进欧拉法 欧欧欧欧拉拉拉拉隐式隐式隐式隐式欧拉欧拉欧拉欧拉显式显式显式显式 节点节点节点节点 x xi i 1.00001.0000   2.00002.0000 4.0000 4.0000   8.00008.0000 1.6000 1.6000   10101 1    3.20003.2000   10101 1 1.00001.00002.50002.5000   1010   1 1 6.25006.2500   1010   2 21.56251.5625   1010   2 23.90633.9063   1010   3 39.76569.7656   1010   4 41.00001.00002.50002.50006.25006.25001.56261.5626   10101 13.90633.9063   10101 19.76569.7656   10101 11.00001.00004.97874.9787   1010   2 22.47882.4788   1010   3 31.23411.2341   1010   4 46.14426.1442   1010   6 63.05903.0590   1010   7 744定义定义若某算法在计算过程中任一步产生的误差在以后若某算法在计算过程中任一步产生的误差在以后的计算中都的计算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的。

讨论数值方法的稳定性讨论数值方法的稳定性, 通常仅限于典型的试验方程通常仅限于典型的试验方程 其中其中 是复数在复平面上是复数在复平面上,当方法稳定时要求变量当方法稳定时要求变量 h的取的取值范围称为方法的值范围称为方法的绝对稳定域绝对稳定域, 它与实轴的交集称为它与实轴的交集称为绝对稳绝对稳定区间定区间. 45将将Euler方法应用于方程方法应用于方程y  =  y, 得到得到 设在设在计算计算yn时产生误差时产生误差 n,计算值计算值 yn = yn +  n, 则则 n将对以后将对以后各节点值计算产生影响各节点值计算产生影响. 记记 ym = ym +  m, m   n, 由上式可知由上式可知误差误差 m满足方程满足方程 可见可见,若要若要| m|<| n|,必须且只须必须且只须|1+ h|<1 ,因此因此Euler法的绝对稳定域为法的绝对稳定域为|1+ h|<1, 绝对绝对稳定区间是稳定区间是-2

方法计算出现严重错误46类似前面分析类似前面分析,可知绝对稳定区域为可知绝对稳定区域为由于由于Re( )<0,所以此不等式对任意步长所以此不等式对任意步长h恒成立恒成立,这是隐这是隐式公式的优点式公式的优点.解出解出yn+1得得 隐式单步方法可类似讨论隐式单步方法可类似讨论.如对梯形公式使用试验方程如对梯形公式使用试验方程,有有 47结论1.单步显式方法的稳定性与步长密切相关, 在一种步长下是稳定的差分公式,取大一点步长就可能是不稳定的.2.收敛性是反映差分公式本身的截断误差对数值解的影响;稳定性是反映计算过程中舍入误差对数值解的影响.只有即收敛又稳定的差分公式才有实用价值.一些常用方法的绝对稳定区间一些常用方法的绝对稳定区间方方 法法方法的阶数方法的阶数稳稳 定定 区区 间间EulerEuler方法方法梯形方法梯形方法改进改进EulerEuler方法方法1 12 22 2(-2 , 0)(-2 , 0)(- (-  , 0) , 0)(-2 , 0)(-2 , 0)48作业作业1．思考题．思考题1中的（中的（a））—（（f））2．习题．习题349。