微分方程初值问题的数值解法精品文稿.ppt

《微分方程初值问题的数值解法精品文稿.ppt》由会员分享，可在线阅读，更多相关《微分方程初值问题的数值解法精品文稿.ppt（41页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、微分方程初值问题的数值解法第1页，本讲稿共41页 引言引言 初值问题的数值解法初值问题的数值解法:求初值问题的解在一系列节点的值求初值问题的解在一系列节点的值 y(xn)的近似值的近似值 yn 的方法的方法.本章数值解法的特点本章数值解法的特点:都是采用都是采用“步步进式进式”,即求解过程顺着节点排列的次序一步步向前推进即求解过程顺着节点排列的次序一步步向前推进.基本知识基本知识:(1(1)定理定理1:1:如果函数如果函数 f(x,y)在区域在区域 上连续上连续,且关于且关于 y 满足满足Lipschitz条件条件常微分方程初值问题:求未知函数求未知函数 y=y(x).第2页，本讲稿共41页此

2、时此时Lipschitz条件显然成立条件显然成立.故常用故常用 在在D上连续有界来代上连续有界来代替替 f(x,y)关于关于 y 满足满足Lipschitz条件条件.注注:如无特别说明如无特别说明,总假设总假设(1)(1)的解存在唯一且足够光滑的解存在唯一且足够光滑.在在f(x,y)对变量对变量 y 可微的情形下可微的情形下,若偏导数若偏导数 连续有界连续有界,则可则可取取L为为除了要保证除了要保证(1)(1)有唯一解外有唯一解外,还需保证微分方程本身是稳定的还需保证微分方程本身是稳定的,即即(1)(1)的解连续依赖于初始值和函数的解连续依赖于初始值和函数 f(x,y).也就是说也就是说,当初

3、当初始值始值 y0 及函数及函数 f(x,y)有微小变化时有微小变化时,只能引起解的微小变化只能引起解的微小变化.(其中其中L 称为称为LipschitzLipschitz常数常数),),则对任何则对任何 ,初值问题初值问题(1)(1)在在 a,b 上存在唯一连续可微解上存在唯一连续可微解 y=y(x).定理定理2:2:如果函数如果函数 f(x,y)在区域在区域 上关于上关于 y 满足满足Lipschitz条件条件,则则(1)是稳定的是稳定的.第3页，本讲稿共41页单步迭代单步迭代:计算计算 yn+1时仅用时仅用 yn;初值问题(1)与下列积分方程的解等价:初值问题的数值解就是求一系列节点上函

4、数 y=y(x)的近似值 .称为步长.一般取等步长 h.多步迭代多步迭代:计算计算 yn+1时除用时除用 yn 外外,还要用到还要用到 yn-1,yn-2,;k 步迭步迭代要用到代要用到 yn-1,yn-2,yn-k+1.显式单步迭代:隐式单步迭代:(2(2)第4页，本讲稿共41页一、Euler方法及其改进 将a,b n 等分,记 微分法:积分法:积分项利用矩形公式计算 1.显式Euler方法()第5页，本讲稿共41页Taylor公式推导:Euler公式几何意义:P1P2Pk也称折线法 P0 xy第6页，本讲稿共41页2.梯形法 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法:取

5、 先用Euler法求出初值 ,即 ,将其代入梯形公式的右端,使之转化为显式公式,即 注:当 f(x,y)关于y满足Lipschitz条件且步长h 满足 直至满足:若采用梯形公式计算若采用梯形公式计算()()中的积分项中的积分项,则有则有类似地,可得()第7页，本讲稿共41页时时,迭代格式迭代格式()收敛收敛.3.改进的Euler方法 把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式),则得改进的Euler方法:或也称为预估-校正法.有时为了方便,预估-校正格式也写成下面形式:第8页，本讲稿共41页二、单步法的局部截断误差及精度 Def 1:先假设 ,再估计误差这种

6、误差称为单步迭代法在 xk+1处的局部截断误差.Def 2:若某种数值方法的局部截断误差为 ,则称该数值方法的精度为P 阶的.注:通常情况下,P 越大,h 越小,则截断误差越小,数值方法越精确.设 10.Euler方法是一阶方法.第9页，本讲稿共41页所以Euler方法为一阶方法.而 20.梯形法是二阶方法.Taylor展开 第10页，本讲稿共41页将 代入上式,得 而代入上式得:当h充分小时,若 ,则可选取 h,使得第11页，本讲稿共41页故梯形法的精度为2.同样可以证明改进的Euler法也是二阶方法.梯形法的局部截断误差为:从而第12页，本讲稿共41页例1:取步长取步长 h=2/10,2/

7、20,2/30,2/40,分别用欧拉法、改进的欧拉法分别用欧拉法、改进的欧拉法和梯形法求解和梯形法求解.解解:记记 f(x,y)=y x y2,xk=k h (k=0,1,2,n)(1).Euler法法:yk+1=yk+h(yk xk yk2)(k=0,1,n)y0=1当当 h=2/10时时,n=10.由由Euler公式可得公式可得:k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.021130.891690.783788第13页，本讲稿共41页(2).改进的Euler法:k01234yk+11.19121.3438

8、41.423481.419051.3473k56789yk+11.237261.114240.994151 0.884751 0.788666(3).梯形法(计算过程略)第14页，本讲稿共41页 n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差 0.1059 0.0521 0.0342 0.0256Euler法误差:改进的Euler法误差:n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差 0.0123 0.0026 0.0011 5.9612e-004第15页，本讲稿共41页预预-校方法校方法,h=0.2时时误差最大值误差最大值:0.0123

9、欧拉方法欧拉方法,h=0.2时时误差最大值误差最大值:0.1059解析解:第16页，本讲稿共41页三、Runge-Kutta 方法1、Taylor 级数法 设初值问题 有解 y(x),由Tayler公式得:令当当 时时,有有 .此时此时为为 p 阶阶Taylor方法方法.p=1时即为时即为Euler公式公式.称之为Taylor级数法.其中例2:取步长 h=0.1,用一阶、二阶和四阶Taylor方法求解下列初值问题第17页，本讲稿共41页解:(1)一阶Taylor法k01234yk+11.11.2211.370081.557791.80046(2)二阶Taylor法k01234yk+11.111

10、.246891.421751.652631.97088第18页，本讲稿共41页(3)四阶Taylor法k01234yk+11.11111.249961.428481.666441.99942第19页，本讲稿共41页记由得称为xk,xk+1上的平均斜率.故2、Runge-Kutta方法只要对K*提供不同的算法,就会得出不同的计算公式.如取则得改进的Euler公式,它是利用xk,xk+1两点的斜率值K1,K2 的算术平均值作为K*,精度比Euler法高.则得Euler公式;取第20页，本讲稿共41页Runge-Kutta法的基本思想:设法在xk,xk+1内多预报几个点的斜率,再将它们的加权平均值作

11、为平均斜率K*一般显式Runge-Kutta公式为:其中 为待定参数,且 .称为r 级Runge-Kutta方法计算公式.注注:式中待定参数的确定式中待定参数的确定:先将先将式右端在式右端在(xk,yk)处展成处展成h的的幂级数幂级数(即将即将 yk+1 展成展成 h 的幂级数的幂级数);再将再将 y(xk+1)作作Taylor 级数展级数展开开;最后比较两式中最后比较两式中hk(k=0,1,2,)的系数的系数,以确定出所有待定参以确定出所有待定参数数.第21页，本讲稿共41页即可得即可得 p 个方程个方程,从而确定出待定参数从而确定出待定参数.代入表达式即可得到计代入表达式即可得到计算公式算

12、公式.如果要求两个表达式的前如果要求两个表达式的前p+1项完全重合项完全重合,即局部截断即局部截断误差达到误差达到 ,则称则称式为式为 p 阶阶 r 级级的的Runge-Kutta方法方法.常常用的是用的是 r=2,3,4 级的级的R-K方法方法,且适当选取参数使得且适当选取参数使得 p=r.如要求:Runge-Kutta方法的推导(以r=2为例):当当r=2 时时第22页，本讲稿共41页则记又第23页，本讲稿共41页(1)常用的二阶Runge-Kutta方法:预估-校正算法(2)这是一个四个参数三个方程的非线性方程组.它有一个自由度.称满足上述方程组的一族公式为二级二阶Runge-Kutta

13、方法.为使局部截断误差为为使局部截断误差为 ,比较上述两式右端同次幂系数比较上述两式右端同次幂系数,应应取取第24页，本讲稿共41页注:二级Runge-Kutta方法的精度最高是二阶的,不可能达到三阶.要提高计算方法的阶,就必须增加预报点.常用的三阶Runge-Kutta方法(r=3):(1)Heun(休恩)方法 中间点方法 (3)三阶Kutta方法 第25页，本讲稿共41页(1)三阶Heun方法 标准(经典)四阶Runge-Kutta方法 (2)常用的四阶Runge-Kutta方法(r=4):第26页，本讲稿共41页(2)称为Gill(吉尔)方法 注:从理论上讲,可以构造任意高阶的计算方法.

14、但事实上,精度的阶数与预报点的个数之间并非等量关系.预报点的个数 r123456789r 10精度的阶数123445667 r-2一般情况下,四阶Runge-Kutta方法已可满足精度要求.第27页，本讲稿共41页例3:用经典Runge-Kutta方法求解下列初值问题(取 h=0.1)解:标准Runge-Kutta公式为:计算结果见下表.为比较在相同计算量条件下近似解的精度,表中列出了Euler法(h=0.025)和改进的Euler法(h=0.05)在相应节点上的计算结果.第28页，本讲稿共41页xiEuler法h=0.025改进Euler法h=0.05经典R-K法h=0.1准确解0.11.1

15、114391.1153801.1155121.1155130.21.2552091.2639141.2642081.2642080.31.4346671.4490891.4495761.4495760.41.6535171.6747561.6754731.6754740.51.9158491.9451711.9461621.9461640.62.2261782.2650402.2663542.2663560.72.5894852.6395612.6412552.6412580.83.0112713.0744793.0766193.0766230.93.4976063.5761443.57880

16、43.5788091.04.0551924.1515734.1548394.154845注:用表中每种方法计算 yi 都需要计算四次 f 的值,即它们的计算量基本相等.第29页，本讲稿共41页四、四、单步单步法的进一步讨论法的进一步讨论收敛性、相容性与稳定性收敛性、相容性与稳定性注:由定义可知,数值方法的收敛性并不涉及计算过程的舍入误差,只与方法的截断误差有关.若格式收敛,则整体截断误差必趋于零.Def:(整体截断误差)称 为某一数值方法在点为某一数值方法在点 xk 处的整体截断误差处的整体截断误差.它不仅与它不仅与 xk 有关有关,也也与与xk-1,xk-2,x1,x0 有关有关.则称该单步

17、法收敛.Def:对满足解存在唯一性条件的初值问题(1),如果一个显式单步法(3)产生的近似解对于任一固定的 ,均有 1.收敛性第30页，本讲稿共41页由于 ,且 关于 y 满足Lipschitz条件,得 则存在常数 c 0 使得 且单步法中函数 关于 y 满足Lipschitz条件,则 定理1:若初值问题的一个单步法的局部截断误差为 记 证:由局部截断误差的定义知第31页，本讲稿共41页故 从而有 故 若 y(x0)=y0,则e0=0,由不等式 得 第32页，本讲稿共41页设单步法为 注注:定理表明定理表明,数值方法的整体截断误差比局部截断误差低一阶数值方法的整体截断误差比局部截断误差低一阶.

18、收敛的方法至少是一阶方法收敛的方法至少是一阶方法.在该定义条件下在该定义条件下,Euler方法是一方法是一阶的阶的,预估预估-校正方法是二阶校正方法是二阶.当当f(x,y)关于关于 y 也满足也满足Lipschitz条件条件,r 级级Runge-Kutta方法中的方法中的 关于关于 y 也满足也满足Lipschitz条件条件,故定理中的条件得到满足故定理中的条件得到满足,解的收敛性得到保证解的收敛性得到保证.由于由于R n,h0(h0),且且 xn为任意点为任意点,故该式相当于用近似方程故该式相当于用近似方程 当x=xn+1固定时,所以有 2.相容性第33页，本讲稿共41页通过在通过在 x=x

19、n 处求解近似方程而获得原方程的近似解处求解近似方程而获得原方程的近似解.因此因此,必必须要求当须要求当h0 时时,近似方程应逼近于原方程近似方程应逼近于原方程.来代替 因此因此,要使要使 h0 时时,近似方程的极限状态为原微分方程近似方程的极限状态为原微分方程,需且需且只需下列极限成立只需下列极限成立:由于 由于假设 是连续函数,故上式可表示为 Def:如果当如果当h0时时,近似方程逼近微分方程近似方程逼近微分方程,则称数值公式则称数值公式与原微分方程相容.相容的充要条件:第34页，本讲稿共41页事实上:Remark:可以证明若单步法的阶大于或等于1,则单步法与微分方程相容;反之,如果单步法

20、与微分方程相容,且 关于 y 满足Lipschitz条件,则单步法至少为一阶方法.(h0)(1)若单步法的阶大于或等于1,由 知 即单步法与微分方程相容.故有 (2)如果单步法与微分方程相容,且 关于y 满足Lipschitz条件,则第35页，本讲稿共41页关于 单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的,但在实际计算时通常会有舍入误差及其积累,数值求解微分方程的过程是一个递推公式,必须考 即与微分方程相容的单步法至少为一阶方法.Remark:在定理条件下,Euler方法、预估-校正方法以及Runge-Kutta方法都与原微分方程相容.中连续,且关于变量y 满

21、足Lipschitz条件,则单步法收敛的充要条件为相容性条件成立.Th1.设增量函数 在区域 3.稳定性第36页，本讲稿共41页 如果数值方法在计算过程中舍入误差的积累越来越大,得不到有效控制,则称其是不稳定的;反之如果计算结果对初始数据的误差及计算过程中的误差不敏感,即舍入误差不增长,则称相应的算法是稳定的.数值方法的稳定性有各种定义,这里仅考虑绝对稳定性概念.虑误差积累能否得到控制.Remark:从上面的定义可知,单步法是绝对稳定的,与模型方程 设某数值方法在节点 xn 处对初值问题的数值解为 yn,实际计算得到的近似解为 ,称为第 n 步数值解得扰动.=-Def:若某种数值方法在计算 y

22、n 时有扰动 但在计算后面的 ym(m n)由 引起的误差 满足 则称该数值方法是绝对稳定的.设设 f(x,y)关于关于 y 满足满足Lipschitz条件条件,下面仅对典型方程下面仅对典型方程(模型方模型方程程)进行讨论进行讨论.(其中其中为复常数为复常数,且且 )第37页，本讲稿共41页(1)Euler显式公式:几个常用公式的稳定性要求误差不增加,即 必须中的复数中的复数以及所用步长以及所用步长 h 有关有关.若对复平面上的某个区域若对复平面上的某个区域G,当当 时单步法绝对收敛时单步法绝对收敛,则称则称G为单步法的绝对收敛域为单步法的绝对收敛域,G与与实轴的交集称为绝对稳定区间实轴的交集

23、称为绝对稳定区间.显然绝对稳定域越大显然绝对稳定域越大,数值方法数值方法的绝对稳定性越好的绝对稳定性越好.将Euler显式公式用于模型方程 ,则有设 yn 有误差 ,参与运算的量为 由此引起的 yn+1 有误差 ,则实际得到近似 yn+1的量为 ,即故有第38页，本讲稿共41页故当当为实数时为实数时,得绝对稳定区间为得绝对稳定区间为(-2,0).即 是保持绝对稳定性对步长 h 所加的限制.因此Euler法的绝对稳定域为(Euler法是条件稳定的):(2)梯形公式:将梯形公式用于模型方程 ,则有由于设第39页，本讲稿共41页当当x=Re(h)0时时,上式右端总小于上式右端总小于1.故梯形法的绝对稳定区域故梯形法的绝对稳定区域为为:Re(h)0,即左半平面即左半平面.(3)四阶经典R-K方法:四阶经典R-K方法用于模型方程 ,则有第40页，本讲稿共41页故四阶经典R-K方法的绝对稳定域为:扰动满足:第41页，本讲稿共41页