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

关于微分方程初值问题的数值解法第1页，共52页，2023年，2月20日，星期四

包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中,自变量的个数只有一个,称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数8.1

引言8.1.1

微分方程知识回顾都是一次的,则称它是线性的,否则称为非线性的。

第2页，共52页，2023年，2月20日，星期四

在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法：

一阶：可分离变量法、齐次方程、一阶线性方程以及伯努利方程

高阶：可降阶方程、常系数齐次线性方程的解法、常系数非齐次线性方程的解法但能求解的常微分方程仍然是有限的，大多数的常微分方程是求不出解析解的。

这个一阶微分方程就不能用初等函数及其积分来表达它的解。8.1.1

微分方程知识回顾

例如第3页，共52页，2023年，2月20日，星期四

从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题

(8.1)

在区间

可以证明,如果函数在带形区域R={a≤x≤b,-∞＜y＜∞}内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使

对R内任意两个都成立,则方程(8.1)的解在a,b上存在且唯一。

上的数值解法。第4页，共52页，2023年，2月20日，星期四

常微分方程初值问题(8.1)式的数值解法，首先要算出精确解y(x)在区间a,b上的一系列离散节点处的函数值的近似值:

相邻两个节点的间距称为步长，本章总是假定h为定数，称为定步长，这时节点可表示为8.1.2

数值方法的基本思想1、数值解法需要把连续性的问题加以离散化，从而求出离散节点处的数值解。第5页，共52页，2023年，2月20日，星期四

描述这类算法，要求给出用已知信息计算的递推公式。建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数进行不同的离散化处理。2、数值解法的基本特点是采用“步进式”：即求解过程按照递推公式顺着节点排列的次序一步一步地向前推进。8.1.2

数值方法的基本思想第6页，共52页，2023年，2月20日，星期四

递推公式通常有两类，一类是计算yi+1时只用到xi+1,xi和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法，其代表是龙格—库塔法。另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到，即前面k步的值，此类方法称为多步法，其代表是亚当斯法。8.1.2

数值方法的基本思想第7页，共52页，2023年，2月20日，星期四一、Euler方法及其改进

将[a,b]n等分,记

微分法:

积分法:

积分项利用矩形公式计算

1.显式Euler方法

(★)第8页，共52页，2023年，2月20日，星期四Taylor公式推导:

第9页，共52页，2023年，2月20日，星期四Oyxy=y(x)(x1,y1)p1p0x0x1x2xixi+1xn-1xnpipi+1pn-1pnp’1p’2p’ip’i+1p’n-1p’n切线p0p1的斜率为f(x0,y0)p2(x2,y2)

欧拉公式的几何意义：第10页，共52页，2023年，2月20日，星期四Euler法的求解过程是:

从初始点P0(即点(x0,y0))出发,作积分曲线y=y(x)在P0点上切线(其斜率为),与直线x=x1相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值）这样就获得了P1点的坐标。当时,得

重复以上过程,就可获得一系列的点:p1,p2,…,pn，相应的可求出y1,y2,…,yn,取第11页，共52页，2023年，2月20日，星期四从图形上看,就获得了一条近似于曲线y=y(x)的折线通常取(常数),则Euler法的计算格式

i=0,1,…,n(7.2)第12页，共52页，2023年，2月20日，星期四2.梯形法

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

取

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

注:

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

直至满足:

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

第13页，共52页，2023年，2月20日，星期四时,迭代格式

(☆)收敛

.

3.改进的Euler方法

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

),则得改进的Euler方法:或也称为预估-校正法.第14页，共52页，2023年，2月20日，星期四有时为了方便,预估-校正格式也写成下面形式:第15页，共52页，2023年，2月20日，星期四

改进的欧拉公式比欧拉公式精度高的原因是：改进欧拉公式用梯形面积代替曲边梯形面积，而欧拉公式用矩形面积代替曲边梯形面积。数值积分的梯形公式比矩形公式的精度高。第16页，共52页，2023年，2月20日，星期四二、单步法的局部截断误差及精度

Def1:

先假设,再估计误差这种误差称为单步迭代法在xk+1处的局部截断误差.Def2:

若某种数值方法的局部截断误差为,则称该数值方法的精度为P阶的.注:

通常情况下,P越大,h越小,则截断误差越小,数值方法越精确.第17页，共52页，2023年，2月20日，星期四所以Euler方法为一阶方法.而设10.Euler方法是一阶方法.第18页，共52页，2023年，2月20日，星期四20.梯形法是二阶方法.Taylor展开

第19页，共52页，2023年，2月20日，星期四将代入上式,得而代入上式得:当h充分小时,若,则可选取h,使得第20页，共52页，2023年，2月20日，星期四故梯形法的精度为2.同样可以证明改进的Euler法也是二阶方法.梯形法的局部截断误差为:从而第21页，共52页，2023年，2月20日，星期四例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第22页，共52页，2023年，2月20日，星期四(2).改进的Euler法:

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

第23页，共52页，2023年，2月20日，星期四n10203040h0.20.10.06670.05误差

0.10590.05210.03420.0256Euler法误差:改进的Euler法误差:n10203040h0.20.10.06670.05误差

0.01230.00260.00115.9612e-004第24页，共52页，2023年，2月20日，星期四预-校方法,h=0.2时误差最大值:0.0123欧拉方法,h=0.2时误差最大值:0.1059解析解:第25页，共52页，2023年，2月20日，星期四三、Runge-Kutta方法1、Taylor级数法

设初值问题有解y(x),由Tayler公式得:令当时,有.此时①为p阶Taylor方法.p=1时即为Euler公式.称之为Taylor级数法.其中例2:

取步长h=0.1,用一阶、二阶和四阶Taylor方法求解下列初值问题①第26页，共52页，2023年，2月20日，星期四解:(1)一阶Taylor法k01234yk+11.11.2211.370081.557791.80046(2)二阶Taylor法k01234yk+11.111.246891.421751.652631.97088第27页，共52页，2023年，2月20日，星期四(3)四阶Taylor法k01234yk+11.11111.249961.428481.666441.99942第28页，共52页，2023年，2月20日，星期四记由得称为[xk,xk+1]上的平均斜率.故2、Runge-Kutta方法只要对K*提供不同的算法,就会得出不同的计算公式.如取则得改进的Euler公式,它是利用xk,xk+1两点的斜率值K1,K2的算术平均值作为K*,精度比Euler法高.则得Euler公式;取第29页，共52页，2023年，2月20日，星期四Runge-Kutta法的基本思想:

设法在[xk,xk+1]内多预报几个点的斜率,再将它们的加权平均值作为平均斜率K*一般显式Runge-Kutta公式为:其中为待定参数,且.称为r级Runge-Kutta方法计算公式.②第30页，共52页，2023年，2月20日，星期四即可得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,…)的系数,以确定出所有待定参数.第31页，共52页，2023年，2月20日，星期四Runge-Kutta方法的推导(以r=2为例):当r=2时记第32页，共52页，2023年，2月20日，星期四则又第33页，共52页，2023年，2月20日，星期四这是一个四个参数三个方程的非线性方程组.它有一个自由度.称满足上述方程组的一族公式为二级二阶Runge-Kutta方法.为使局部截断误差为,比较上述两式右端同次幂系数,应取第34页，共52页，2023年，2月20日，星期四(1)常用的二阶Runge-Kutta方法:预估-校正算法(2)中间点方法

第35页，共52页，2023年，2月20日，星期四注:

二级Runge-Kutta方法的精度最高是二阶的,不可能达到三阶.要提高计算方法的阶,就必须增加预报点.常用的三阶Runge-Kutta方法(r=3):

(1)Heun(休恩)方法

(3)三阶Kutta方法

第36页，共52页，2023年，2月20日，星期四(1)三阶Heun方法

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

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

第37页，共52页，2023年，2月20日，星期四(2)称为Gill(吉尔)方法

注:

从理论上讲,可以构造任意高阶的计算方法.但事实上,精度的阶数与预报点的个数之间并非等量关系.预报点的个数r123456789r≥10精度的阶数123445667≤r-2一般情况下,四阶Runge-Kutta方法已可满足精度要求.第38页，共52页，2023年，2月20日，星期四例3:

用经典Runge-Kutta方法求解下列初值问题(取h=0.1)解:标准Runge-Kutta公式为:计算结果见下表.为比较在相同计算量条件下近似解的精度,表中列出了Euler法(h=0.025)和改进的Euler法(h=0.05)在相应节点上的计算结果.第39页，共52页，2023年，2月20日，星期四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的值,即它们的计算量基本相等.第40页，共52页，2023年，2月20日，星期四四、单步法的进一步讨论—收敛性、相容性与稳定性注:

由定义可知,数值方法的收敛性并不涉及计算过程的舍入误差,只与方法的截断误差有关.若格式收敛,则整体截断误差必趋于零.Def:

(整体截断误差)

称为某一数值方法在点xk处的整体截断误差.它不仅与xk有关,也与xk-1,xk-2,…,x1,x0有关.则称该单步法收敛.Def:

对满足解存在唯一性条件的初值问题(1),如果一个显式单步法(3)产生的近似解对于任一固定的,均有1.收敛性第41页，共52页，2023年，2月20日，星期四由于,且关于y满足Lipschitz条件,得则存在常数c>0使得且单步法中函数关于y满足Lipschitz条件,则定理1:

若初值问题的一个单步法的局部截断误差为记证:由局部截断误差的定义知第42页，共52页，2023年，2月20日，星期四故从而有故若y(x0)=y0,则e0=0,由不等式得第43页，共52页，2023年，2月20日，星期四设单步法为注:定理表明,数值方法的整体截断误差比局部截断误差低一阶.收敛的方法至少是一阶方法.在该定义条件下,Euler方法是一阶的,预估-校正方法是二阶.当f(x,y)关于y也满足Lipschitz条件,r级Runge-Kutta方法中的φ

关于y也满足Lipschitz条件,故定理中的条件得到满足,解的收敛性得到保证.由于Rn,h→0(h→0),且xn为任意点,故该式相当于用近似方程当x=xn+1固定时,,所以有2.相容性第44页，共52页，2023年，2月20日，星期四通过在x=xn处求解近似方程而获得原方程的近似解.因此,必须要求当h→0

时,近似方程应逼近于原方程.来代替因此,要使h→0

时,近似方程的极限状态为原微分方程,需且只需下列极限成立:由于由于假设是连续函数,故上式可表示为Def:

如果当h→0时,

近似方程逼近微分方程,则称数值公式与原微分方程相容.相容的充要条件:

第45页，共52页，2023年，2月20日，星期四事实上:Remark:

可以证明若单步法的阶大于或等于1,则单步法与微分方程相容;反之,如果单步法与微分方程相容,且关于y满足Lipschitz条件,则单步法至少为一阶方法.(h→0)(1)若单步法的阶大于或等于1,由知即单步法与微分方程相容.故有(2)如果单步法与微分方程相容,且关于y满足Lipschitz条件,则第46页，共52页，2023年，2月20日，星期四关于单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的,但在实际计算时通常会有舍入误差及其积累,数值求解微分方程的过程是一个递推公式,必须考

即与微分方程相容的单步法至少为一阶方法.Remark:

在定理条件下,Euler方法、预估-校正方法以及Runge-Kutta方法都与原微分方程相容.中连续,且关于变量y满足Lipschitz条件,则单步法收敛的充要条件为相容性条件成立.Th1.

设增量函数在区域3.稳定性第47页，共52页，2023年，2月20日，星期四

如果数值方法在计算过程中舍入误差的积累越来越大,得不到有效控制,则称其是不稳定的;反之如果计算结果对初始数据的误差及计算过程中的误差不敏感,即舍入误差不增长,则称相应的算法是稳定的.数值方法的稳定性有各种定义,这里仅考虑绝对稳定性概念.虑误差积累能否得到控制.Remark:从上面的定义可知,单步法是绝对稳定的,与模型方程

设某数值方法在节点xn处对初值问题的数值解为yn,实际计算得到的近似解为,