微分方程初值问题的数值解法课件.ppt

包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数,8.1 引 言,8.1.1 微分方程知识回顾,都是一次的,则称它是线性的,否则称为非线性的。,在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法：一阶：可分离变量法、齐次方程、一阶线性方程以及伯努利方程 高阶：可降阶方程、常系数齐次线性方程的解法、常系数非齐次线性方程的解法 但能求解的常微分方程仍然是有限的，大多数的常微分方程是求不出解析解的。,这个一阶微分方程就不能用初等函数及其积分来表达它的解。,8.1.1 微分方程知识回顾,例如,从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题,(8.1),在区间,可以证明,如果函数在带形区域 R=axb,-y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使,对R内任意两个 都成立,则方程(8.1)的解 在a,b上存在且唯一。,上的数值解法。,常微分方程初值问题(8.1)式的数值解法，首先要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值:相邻两个节点的间距 称为步长，本章总是假定h为定数，称为定步长，这时节点可表示为,8.1.2 数值方法的基本思想,1、数值解法需要把连续性的问题加以离散化，从而求出离散节点处的数值解。,描述这类算法，要求给出用已知信息 计算 的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数 进行不同的离散化处理。,2、数值解法的基本特点是采用“步进式”：即求解过程按照递推公式顺着节点排列的次序一步一步地向前推进。,8.1.2 数值方法的基本思想,递推公式通常有两类，一类是计算yi+1时只用到xi+1,xi 和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法，其代表是龙格库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法，其代表是亚当斯法。,8.1.2 数值方法的基本思想,一、Euler方法及其改进,将a,b n 等分,记,微分法:,积分法:,积分项利用矩形公式计算,1.显式Euler方法,(),Taylor公式推导:,p1,p2,pi,pi+1,pn-1,pn,切线p0p1的斜率为f(x0,y0),欧拉公式的几何意义：,Euler法的求解过程是:从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上切线(其斜率为),与直线x=x1相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值）,这样就获得了P1点的坐标。,当 时,得,重复以上过程,就可获得一系列的点:p1,p2,pn，相应的可求出y1,y2,yn,取,从图形上看,就获得了一条近似于曲线y=y(x)的折线,通常取(常数),则Euler法的计算格式,i=0,1,n(7.2),2.梯形法,称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法:,取,先用Euler法求出初值,即,将其代入梯形公式的右端,使之转化为显式公式,即,注:当 f(x,y)关于y满足Lipschitz条件且步长h 满足,直至满足:,若采用梯形公式计算()中的积分项,则有,类似地,可得,(),时,迭代格式()收敛.,3.改进的Euler方法,把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式),则得改进的Euler方法:,或,也称为预估-校正法.,有时为了方便,预估-校正格式也写成下面形式:,改进的欧拉公式比欧拉公式精度高的原因是：改进欧拉公式用梯形面积代替曲边梯形面积，而欧拉公式用矩形面积代替曲边梯形面积。数值积分的梯形公式比矩形公式的精度高。,二、单步法的局部截断误差及精度,Def 1:先假设,再估计误差这种误差称为单步迭代法在 xk+1处的局部截断误差.,Def 2:若某种数值方法的局部截断误差为,则称该数值方法的精度为P 阶的.,注:通常情况下,P 越大,h 越小,则截断误差越小,数值方法越精确.,所以Euler方法为一阶方法.,而,设,10.Euler方法是一阶方法.,20.梯形法是二阶方法.,Taylor展开,将 代入上式,得,而,代入上式得:,当h充分小时,若,则可选取 h,使得,故梯形法的精度为2.,同样可以证明改进的Euler法也是二阶方法.,梯形法的局部截断误差为:,从而,例1:,取步长 h=2/10,2/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公式可得:,(2).改进的Euler法:,(3).梯形法(计算过程略),Euler法误差:,改进的Euler法误差:,预-校方法,h=0.2时误差最大值:0.0123,欧拉方法,h=0.2时误差最大值:0.1059,解析解:,三、Runge-Kutta 方法,1、Taylor 级数法,设初值问题 有解 y(x),由Tayler公式得:,令,当 时,有.此时为 p 阶Taylor方法.p=1时即为Euler公式.,称之为Taylor级数法.其中,例2:取步长 h=0.1,用一阶、二阶和四阶Taylor方法求解下列初值问题,解:(1)一阶Taylor法,(2)二阶Taylor法,(3)四阶Taylor法,记,由,得,称为xk,xk+1上的平均斜率.故,2、Runge-Kutta方法,只要对K*提供不同的算法,就会得出不同的计算公式.如取,则得改进的Euler公式,它是利用xk,xk+1两点的斜率值K1,K2 的算术平均值作为K*,精度比Euler法高.,则得Euler公式;取,Runge-Kutta法的基本思想:,设法在xk,xk+1内多预报几个点的斜率,再将它们的加权平均值作为平均斜率K*,一般显式Runge-Kutta公式为:,其中 为待定参数,且.称为r 级Runge-Kutta方法计算公式.,即可得 p 个方程,从而确定出待定参数.代入表达式即可得到计算公式.如果要求两个表达式的前p+1项完全重合,即局部截断误差达到,则称式为 p 阶 r 级的Runge-Kutta方法.常用的是 r=2,3,4 级的R-K方法,且适当选取参数使得 p=r.,如要求:,注:式中待定参数的确定:先将式右端在(xk,yk)处展成h的幂级数(即将 yk+1 展成 h 的幂级数);再将 y(xk+1)作Taylor 级数展开;最后比较两式中hk(k=0,1,2,)的系数,以确定出所有待定参数.,Runge-Kutta方法的推导(以r=2为例):,当r=2 时,记,则,又,这是一个四个参数三个方程的非线性方程组.它有一个自由度.称满足上述方程组的一族公式为二级二阶Runge-Kutta方法.,为使局部截断误差为,比较上述两式右端同次幂系数,应取,(1),常用的二阶Runge-Kutta方法:,预估-校正算法,(2),中间点方法,注:二级Runge-Kutta方法的精度最高是二阶的,不可能达到三阶.要提高计算方法的阶,就必须增加预报点.,常用的三阶Runge-Kutta方法(r=3):,(1),Heun(休恩)方法,(3),三阶Kutta方法,(1),三阶Heun方法,标准(经典)四阶Runge-Kutta方法,(2),常用的四阶Runge-Kutta方法(r=4):,(2),称为Gill(吉尔)方法,注:从理论上讲,可以构造任意高阶的计算方法.但事实上,精度的阶数与预报点的个数之间并非等量关系.,一般情况下,四阶Runge-Kutta方法已可满足精度要求.,例3:用经典Runge-Kutta方法求解下列初值问题(取 h=0.1),解:,标准Runge-Kutta公式为:,计算结果见下表.为比较在相同计算量条件下近似解的精度,表中列出了Euler法(h=0.025)和改进的Euler法(h=0.05)在相应节点上的计算结果.,注:用表中每种方法计算 yi 都需要计算四次 f 的值,即它们的计算量基本相等.,四、单步法的进一步讨论收敛性、相容性与稳定性,注:由定义可知,数值方法的收敛性并不涉及计算过程的舍入误差,只与方法的截断误差有关.若格式收敛,则整体截断误差必趋于零.,1.收敛性,由于,且 关于 y 满足Lipschitz条件,得,则存在常数 c 0 使得,且单步法中函数 关于 y 满足Lipschitz条件,则,定理1:若初值问题的一个单步法的局部截断误差为,记,证:由局部截断误差的定义知,故,从而有,故,若 y(x0)=y0,则e0=0,由不等式,得,设单步法为,注:定理表明,数值方法的整体截断误差比局部截断误差低一阶.收敛的方法至少是一阶方法.在该定义条件下,Euler方法是一阶的,预估-校正方法是二阶.当f(x,y)关于 y 也满足Lipschitz条件,r 级Runge-Kutta方法中的 关于 y 也满足Lipschitz条件,故定理中的条件得到满足,解的收敛性得到保证.,由于R n,h0(h0),且 xn为任意点,故该式相当于用近似方程,当x=xn+1固定时,所以有,2.相容性,通过在 x=xn 处求解近似方程而获得原方程的近似解.因此,必须要求当h0 时,近似方程应逼近于原方程.,来代替,因此,要使 h0 时,近似方程的极限状态为原微分方程,需且只需下列极限成立:,由于,由于假设 是连续函数,故上式可表示为,相容的充要条件:,事实上:,Remark:可以证明若单步法的阶大于或等于1,则单步法与微分方程相容;反之,如果单步法与微分方程相容,且 关于 y 满足Lipschitz条件,则单步法至少为一阶方法.,(h0),(1)若单步法的阶大于或等于1,由,知,即单步法与微分方程相容.,故有,(2)如果单步法与微分方程相容,且 关于y 满足Lipschitz条件,则,关于 单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的,但在实际计算时通常会有舍入误差及其积累,数值求解微分方程的过程是一个递推公式,必须考,即与微分方程相容的单步法至少为一阶方法.,Remark:在定理条件下,Euler方法、预估-校正方法以及Runge-Kutta方法都与原微分方程相容.,中连续,且关于变量y 满足Lipschitz条件,则单步法收敛的充要条件为相容性条件成立.,Th1.设增量函数 在区域,3.稳定性,如果数值方法在计算过程中舍入误差的积累越来越大,得不到有效控制,则称其是不稳定的;反之如果计算结果对初始数据的误差及计算过程中的误差不敏感,即舍入误差不增长,则称相应的算法是稳定的.数值方法的稳定性有各种定义,这里仅考虑绝对稳定性概念.,虑误差积累能否得到控制.,Remark:从上面的定义可知,单步法是绝对稳定的,与模型方程,Def:若某种数值方法在计算 yn 时有扰动 但在计算后面的 ym(m n)由 引起的误差 满足 则称该数值方法是绝对稳定的.,设 f(x,y)关于 y 满足Lipschitz条件,下面仅对典型方程(模型方程)进行讨论.(其中为复常数,且),(1)Euler显式公式:,几个常用公式的稳定性,要求误差不增加,即 必须,中的复数以及所用步长 h 有关.若对复平面上的某个区域G,当 时单步法绝对收敛,则称G为单步法的绝对收敛域,G与实轴的交集称为绝对稳定区间.显然绝对稳定域越大,数值方法的绝对稳定性越好.,将Euler显式公式用于模型方程,则有,设 yn 有误差,参与运算的量为 由此引起的 yn+1 有误差,则实际得到近似 yn+1的量为,即,故有,故,当为实数时,得绝对稳定区间为(-2,0).,即,是保持绝对稳定性对步长 h 所加的限制.因此Euler法的绝对稳定域为(Euler法是条件稳定的):,(2)梯形公式:,将梯形公式用于模型方程,则有,由于,设,当x=Re(h)0时,上式右端总小于1.故梯形法的绝对稳定区域为:Re(h)0,即左半平面.,(3)四阶经典R-K方法:,四阶经典R-K方法用于模型方程,则有,