微分方程初值问题的数值解法

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

第二页，共四十一页，编辑于2023年，星期二此时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)是稳定的.第三页，共四十一页，编辑于2023年，星期二单步迭代:计算yn+1时仅用yn;初值问题(1)与下列积分方程的解等价:

初值问题的数值解就是求一系列节点上函数y=y(x)的近似值.

称为步长.一般取等步长h.多步迭代:计算yn+1时除用yn

外,还要用到yn-1,yn-2,…;k步迭代要用到yn-1,yn-2,…,yn-k+1.

显式单步迭代:

隐式单步迭代:

(2)第四页，共四十一页，编辑于2023年，星期二一、Euler方法及其改进

将[a,b]n等分,记

微分法:

积分法:

积分项利用矩形公式计算

1.显式Euler方法

(★)第五页，共四十一页，编辑于2023年，星期二Taylor公式推导:

Euler公式几何意义:

P1P2Pk也称折线法

P0xy第六页，共四十一页，编辑于2023年，星期二2.梯形法

称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法:

取

先用Euler法求出初值,即,将其代入梯形公式的右端,使之转化为显式公式,即

注:

当f(x,y)关于y满足Lipschitz条件且步长h满足

直至满足:

若采用梯形公式计算(★)中的积分项,则有类似地,可得(☆)

第七页，共四十一页，编辑于2023年，星期二时,迭代格式

(☆)收敛

.

3.改进的Euler方法

把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式

),则得改进的Euler方法:或也称为预估-校正法.有时为了方便,预估-校正格式也写成下面形式:第八页，共四十一页，编辑于2023年，星期二二、单步法的局部截断误差及精度

Def1:先假设,再估计误差这种误差称为单步迭代法在xk+1处的局部截断误差.Def2:若某种数值方法的局部截断误差为,则称该数值方法的精度为P阶的.注:通常情况下,P越大,h越小,则截断误差越小,数值方法越精确.设10.Euler方法是一阶方法.第九页，共四十一页，编辑于2023年，星期二所以Euler方法为一阶方法.而20.梯形法是二阶方法.Taylor展开

第十页，共四十一页，编辑于2023年，星期二将代入上式,得而代入上式得:当h充分小时,若,则可选取h,使得第十一页，共四十一页，编辑于2023年，星期二故梯形法的精度为2.同样可以证明改进的Euler法也是二阶方法.梯形法的局部截断误差为:从而第十二页，共四十一页，编辑于2023年，星期二例1:

取步长h=2/10,2/20,2/30,2/40,分别用欧拉法、改进的欧拉法和梯形法求解.解:记f(x,y)=y－xy2,xk=kh(k=0,1,2,···,n)(1).Euler法:yk+1=yk+h(yk－xkyk2)(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第十三页，共四十一页，编辑于2023年，星期二(2).改进的Euler法:

k01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3).梯形法(计算过程略)

第十四页，共四十一页，编辑于2023年，星期二n10203040h0.20.10.06670.05误差0.10590.05210.03420.0256Euler法误差:改进的Euler法误差:n10203040h0.20.10.06670.05误差0.01230.00260.00115.9612e-004第十五页，共四十一页，编辑于2023年，星期二预-校方法,h=0.2时误差最大值:0.0123欧拉方法,h=0.2时误差最大值:0.1059解析解:第十六页，共四十一页，编辑于2023年，星期二三、Runge-Kutta方法1、Taylor级数法设初值问题有解y(x),由Tayler公式得:令当时,有.此时①为p阶Taylor方法.p=1时即为Euler公式.称之为Taylor级数法.其中例2:取步长h=0.1,用一阶、二阶和四阶Taylor方法求解下列初值问题①第十七页，共四十一页，编辑于2023年，星期二解:(1)一阶Taylor法k01234yk+11.11.2211.370081.557791.80046(2)二阶Taylor法k01234yk+11.111.246891.421751.652631.97088第十八页，共四十一页，编辑于2023年，星期二(3)四阶Taylor法k01234yk+11.11111.249961.428481.666441.99942第十九页，共四十一页，编辑于2023年，星期二记由得称为[xk,xk+1]上的平均斜率.故2、Runge-Kutta方法只要对K*提供不同的算法,就会得出不同的计算公式.如取则得改进的Euler公式,它是利用xk,xk+1两点的斜率值K1,K2的算术平均值作为K*,精度比Euler法高.则得Euler公式;取第二十页，共四十一页，编辑于2023年，星期二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,…)的系数,以确定出所有待定参数.第二十一页，共四十一页，编辑于2023年，星期二即可得p个方程,从而确定出待定参数.代入表达式即可得到计算公式.如果要求两个表达式的前p+1项完全重合,即局部截断误差达到,则称②式为p阶r级的Runge-Kutta方法.常用的是r=2,3,4级的R-K方法,且适当选取参数使得p=r

.如要求:Runge-Kutta方法的推导(以r=2为例):当r=2时第二十二页，共四十一页，编辑于2023年，星期二则记又第二十三页，共四十一页，编辑于2023年，星期二(1)常用的二阶Runge-Kutta方法:预估-校正算法(2)这是一个四个参数三个方程的非线性方程组.它有一个自由度.称满足上述方程组的一族公式为二级二阶Runge-Kutta方法.为使局部截断误差为,比较上述两式右端同次幂系数,应取第二十四页，共四十一页，编辑于2023年，星期二注:二级Runge-Kutta方法的精度最高是二阶的,不可能达到三阶.要提高计算方法的阶,就必须增加预报点.常用的三阶Runge-Kutta方法(r=3):

(1)Heun(休恩)方法

中间点方法

(3)三阶Kutta方法

第二十五页，共四十一页，编辑于2023年，星期二(1)三阶Heun方法

标准(经典)四阶Runge-Kutta方法

(2)常用的四阶Runge-Kutta方法(r=4):

第二十六页，共四十一页，编辑于2023年，星期二(2)称为Gill(吉尔)方法

注:从理论上讲,可以构造任意高阶的计算方法.但事实上,精度的阶数与预报点的个数之间并非等量关系.预报点的个数r123456789r≥10精度的阶数123445667≤r-2一般情况下,四阶Runge-Kutta方法已可满足精度要求.第二十七页，共四十一页，编辑于2023年，星期二例3:用经典Runge-Kutta方法求解下列初值问题(取h=0.1)解:标准Runge-Kutta公式为:计算结果见下表.为比较在相同计算量条件下近似解的精度,表中列出了Euler法(h=0.025)和改进的Euler法(h=0.05)在相应节点上的计算结果.第二十八页，共四十一页，编辑于2023年，星期二xiEuler法h=0.025改进Euler法h=0.05经典R-K法h=0.1准确解0.11.1114391.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.5788043.5788091.04.0551924.1515734.1548394.154845注:用表中每种方法计算yi都需要计算四次f的值,即它们的计算量基本相等.第二十九页，共四十一页，编辑于2023年，星期二四、单步法的进一步讨论—收敛性、相容性与稳定性注:由定义可知,数值方法的收敛性并不涉及计算过程的舍入误差,只与方法的截断误差有关.若格式收敛,则整体截断误差必趋于零.Def:

(整体截断误差)称为某一数值方法在点xk处的整体截断误差.它不仅与xk有关,也与xk-1,xk-2,…,x1,x0有关.则称该单步法收敛.Def:对满足解存在唯一性条件的初值问题(1),如果一个显式单步法(3)产生的近似解对于任一固定的,均有1.收敛性第三十页，共四十一页，编辑于2023年，星期二由于,且关于y满足Lipschitz条件,得则存在常数c>0使得且单步法中函数关于y满足Lipschitz条件,则定理1:若初值问题的一个单步法的局部截断误差为记证:由局部截断误差的定义知第三十一页，共四十一页，编辑于2023年，星期二故从而有故若y(x0)=y0,则e0=0,由不等式得第三十二页，共四十一页，编辑于2023年，星期二设单步法为注:定理表明,数值方法的整体截断误差比局部截断误差低一阶.收敛的方法至少是一阶方法.在该定义条件下,Euler方法是一阶的,预估-校正方法是二阶.当f(x,y)关于y也满足Lipschitz条件,r级Runge-Kutta方法中的φ关于y也满足Lipschitz条件,故定理中的条件得到满足,解的收敛性得到保证.由于Rn,h→0(h→0),且xn为任意点,故该式相当于用近似方程当x=xn+1固定时,,所以有2.相容性第三十三页，共四十一页，编辑于2023年，星期二通过在x=xn处求解近似方程而获得原方程的近似解.因此,必须要求当h→0时,近似方程应逼近于原方程.来代替因此,要使h→0时,近似方程的极限状态为原微分方程,需且只需下列极限成立:由于由于假设是连续函数,故上式可表示为Def:如果当h→0时,近似方程逼近微分方程,则称数值公式与原微分方程相容.相容的充要条件:

第三十四页，共四十一页，编辑于2023年，星期二事实上:Remark:可以证明若单步法的阶大于或等于1,则单步法与微分方程相容;反之,如果单步法与微分方程相容,且关于y满足Lipschitz条件,则单步法至少为一阶方法.(h→0)(1)若单步法的阶大于或等于1,由知即单步法与微分方程相容.故有(2)如果单步法与微分方程相容,且关于y满足Lipschitz条件,则第三十五页，共四十一页，编辑于2023年，星期二关于单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的,但在实际计算时通常会有舍入误差及其积累,数值求解微分方程的过程是一个递推公式,必须考即与微分方程相容的单步法至少为一阶方法.Remark:在定理条件下,Euler方法、预估-校正方法以及Runge-Kutta方法都与原微分方程相容.中连续,且关于变量y满足Lipschitz条件,则单步法收敛的充要条件为相容性条件成立.Th1.设增量函数在区域3.稳定性第三十六页，共四十一页，编辑于2023年，星期二如果数值方法在计算过程中舍入误差的积累越来越大,得不到有效控制,则称其是不稳定的;反之如果计算结果对初始数据的误差及计算过程中的误差不敏感,即舍入误差不增长,则称相应的算法是稳定的.数值方法的稳定性有各种定义,这里仅考虑绝对稳定性概念.虑误差积累能否得到控制.Remark:从上面的定义可知,单步法是绝对稳定的,与模型方程设某数值方法在节点xn处对初值问题的数值解为yn,实际计算得到的近似解为,称为第n步数值解得扰动.=-Def:若某种数值方法在计算yn

时有扰动但在计算后面的ym(m>n)由引起的误差满足