第九章 微分方程初值问题的数值解法,内容提纲,引言 Euler法及其改进 Runge-Kutta方法 线性多步法 误差分析 数值解法的收敛性、相容性和稳定性 边值问题数值解法简介,引言,初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进式”,即求解过程顺着节点排列的次序一步步向前推进.,基本知识:,(1),定理1: 如果函数 f (x , y)在区域 上连续,且关于 y 满足Lipschitz条件,常微分方程初值问题:,求未知函数 y= y (x) .,此时Lipschitz条件显然成立. 故常用 在D上连续有界来代替 f (x , y)关于 y 满足Lipschitz条件.,注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f (x , y)对变量 y 可微的情形下, 若偏导数 连续有界, 则可取L为,除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即(1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.,(其中L 称为Lipschitz常数),则对任何 ,初值问题(1)在a ,b上存在唯一连续可微解 y = y (x).,定理2: 如果函数 f (x , y)在区域 上关于 y 满足Lipschitz条件, 则(1)是稳定的.,单步迭代: 计算 yn+1时仅用 yn ;,初值问题(1)与下列积分方程的解等价:,初值问题的数值解就是求一系列节点 上函数 y= y (x)的近似值 . 称为步长. 一般取等步长 h .,多步迭代: 计算 yn+1时除用 yn 外, 还要用到 yn-1, yn-2 ,; k 步迭代要用到 yn-1, yn-2 , yn-k+1 .,显式单步迭代:,隐式单步迭代:,(2),一、Euler方法及其改进,将a , b n 等分, 记,微分法:,积分法:,积分项利用矩形公式计算,1. 显式Euler方法,(),Taylor公式推导:,Euler公式几何意义:,2. 梯形法,称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法:,取,先用Euler法求出初值 ,即 ,将其代入梯形公式的右端,使之转化为显式公式,即,注: 当 f (x , y)关于y满足Lipschitz条件且步长h 满足,直至满足:,若采用梯形公式计算()中的积分项,则有,类似地,可得,(),时,迭代格式 () 收敛 .,3. 改进的Euler方法,把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式 ),则得改进的Euler方法:,或,也称为预估-校正法.,有时为了方便,预估-校正格式也写成下面形式:,二、单步法的局部截断误差及精度,Def 1: 先假设 ,再估计误差 这种误差称为单步迭代法在 xk+1处的局部截断误差.,Def 2: 若某种数值方法的局部截断误差为 ,则称该数值方法的精度为P 阶的.,注: 通常情况下,P 越大, h 越小,则截断误差越小,数值方法越精确.,设,10 .Euler方法是一阶方法.,所以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方法计算公式.,注: 式中待定参数的确定: 先将式右端在(xk , yk ) 处展成h的幂级数(即将 yk+1 展成 h 的幂级数); 再将 y(xk+1) 作Taylor 级数展开;最后比较两式中hk ( k=0,1,2,)的系数, 以确定出所有待定参数.,即可得 p 个方程, 从而确定出待定参数. 代入表达式即可得到计算公式. 如果要求两个表达式的前p+1项完全重合, 即局部截断误差达到 , 则称式为 p 阶 r 级的Runge-Kutta方法. 常用的是 r =2,3,4 级的R-K方法, 且适当选取参数使得 p = r .,如要求:,Runge-Kutta方法的推导(以r =2为例):,当r =2 时,则,记,又,(1),常用的二阶Runge-Kutta方法:,预估-校正算法,(2),这是一个四个参数三个方程的非线性方程组. 它有一个自由度. 称满足上述方程组的一族公式为二级二阶Runge-Kutta方法.,为使局部截断误差为 ,比较上述两式右端同次幂系数,应 取,注: 二级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, 当