《第9章常微分方程初值问题数值解法幻灯片课件》由会员分享,可在线阅读,更多相关《第9章常微分方程初值问题数值解法幻灯片课件(91页珍藏版)》请在金锄头文库上搜索。
1、第9章 常微分方程初值问题数值解法,9.1 引言 9.2 简单的数值方法与基本概念 9.3 龙格-库塔方法 9.4 单步法的收敛性与稳定性 9.5 线性多步法 9.6 方程组和高阶方程,9.1 引 言,科学技术中常常需要求解常微分方程的定解问题. 这类问题最简单的形式,是本章将着重考察的一阶方程的初值问题,我们知道,只有f(x, y)适当光滑譬如关于y满足利普希茨(Lipschitz)条件,理论上就可以保证初值问题的解yf(x)存在并且唯一.,虽然求解常微分方程有各种各样的解析方法,但解析方法只能用来求解一些特殊类型的方程,实际问题中归结出来的微分方程主要靠数值解法.,所谓数值解法, 就是寻求
2、解y(x)在一系列离散节点,上的近似值 y1,y2,yn,yn+1,. 相邻两个节点的间距hn=xn+1-xn称为步长. 今后如不特别说明,总是假定 hi=h(i=1,2,)为定数, 这时节点为xn=x0+nh(i=0,1,2,) (等距节点).,初值问题的数值解法有个基本特点,他们都采取“步进式”,即求解过程顺着节点排列的次序一步一步地向前推进. 描述这类算法,只要给出用已知信息yn,yn-1,yn-2,计算yn+1的递推公式.,首先,要对微分方程离散化,建立求解数值解的递推公式. 一类是计算yn+1时只用到前一点的值yn,称为单步法. 另一类是用到yn+1前面 k 点的值yn,yn-1,
3、yn-k+1,称为k步法. 其次,要研究公式的局部截断误差和阶,数值解yn与精确解y(xn)的误差估计及收敛性,还有递推公式的计算稳定性等问题.,一般地,设已做出该折线的顶点Pn,过Pn(xn, yn)依方向场的方向再推进到Pn+1(xn+1, yn+1),显然两个顶点Pn,Pn+1的坐标有关系,这就是著名的(显式)欧拉(Euler)公式. 若初值y0已知,则依公式(2.1)可逐次逐步算出各点数值解.,即,例1 用欧拉公式求解初值问题,解 取步长h=0.1,欧拉公式的具体形式为,其中xn=nh=0.1n (n=0,1,10), 已知y0 =1, 由此式可得,依次计算下去,部分计算结果见下表.,
4、与准确解 相比,可看出欧拉公式的计算结果精度很差.,欧拉公式具有明显的几何意义, 就是用折线近似代替方程的解曲线,因而常称公式(2.1)为欧拉折线法.,还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上,那么,,按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线,可见欧拉方法是相当粗糙的.,为了分析计算公式的精度,通常可用泰勒展开将y(xn+1)在xn处展开,则有,在yn=y(xn)的前提下,f(xn,yn )=f(xn,y(xn)=y(xn).于是可得欧拉法(2.1)的公
5、式误差为,称为此方法的局部截断误差.,如果对方程(1.1)从xn到xn+1积分,得,右端积分用左矩形公式hf(xn,y(xn)近似,再以yn代替y(xn),yn+1代替y(xn+1)也得到欧拉公式(2.1),局部截断误差也是(2.3).,称为(隐式)后退的欧拉公式.,如果右端积分用右矩形公式hf(xn+1,y(xn+1)近似,则得到另一个公式,后退的欧拉公式与欧拉公式有着本质的区别, 后者是关于yn+1的一个直接计算公式,这类公式称作是显式的;前者公式的右端含有未知的yn+1,它实际上是关于yn+1的一个函数方程,这类方程称作是隐式的.,显式与隐式两类方法各有特点,考虑到数值稳定性等其他因素,
6、人们有时需要选用隐式方法,但使用显式算法远比隐式方便.,隐式方程通常用迭代法求解,而迭代过程的实质是逐步显式化.,设用欧拉公式,给出迭代初值 ,用它代入(2.5)式的右端,使之转化为显式,直接计算得,然后再用 代入(2.5)式,又有,如此反复进行,得,由于f(x, y)对y满足Lipschitz条件(1.3). 由(2.6)减(2.5)得,由此可知,只要hL1,迭代法(2.6)就收敛到解.关于后退欧拉方法的公式误差,从积分公式看到它与欧拉法是相似的.,9.2.2 梯形方法,为得到比欧拉法精度高的计算公式,在等式(2.4) 右端积分用梯形求积公式近似, 并用yn代替y(xn), yn+1代替y(
7、xn+1),则得,称为矩形方法.,矩形方法是隐式单步法,用迭代法求解,同后退的欧拉方法一样,仍用欧拉法提供迭代初值,则矩形迭代公式为,为了分析迭代过程的收敛性, 将(2.7)与(2.8)相减, 得,于是有,使得,则当k时有 , 这说明迭代过程(2.8)是收敛的.,9.2.3 单步法的局部截断误差与阶,初值问题(1.1),(1.2)的单步法可用一般形式表示为,其中多元函数与f(x, y )有关,当含有yn+1时,方法是隐式的,若不含yn+1则为显式方法,所以显式单步法可表示为,(x, y, h)称为增量函数,例如对欧拉法(2.1)有,它的局部截断误差已由(2.3)给出, 对一般显式单步法则可如下
8、定义.,定义1 设y(x)是初值问题(1.1),(1.2)的准确解, 称,为显式单步法(2.10)的局部截断误差.,Tn+1之所以称为局部的,是假设在xn前各步没有误差.当yn=y(xn)时,计算一步,则有,所以,局部截断误差可理解为用方法(2.10)计算一步的误差,也即公式(2.10)中用准确解y(x)代替数值解产生的公式误差. 根据定义, 显然欧拉法的局部截断误差为,即为(2.3)的结果. 这里 称为局部截断误差主项. 显然Tn+1=O(h2). 一般情形的定义如下,定义2 设y (x)是初值问题的准确解,若存在最大整数p使显式单步法(2.10)的局部截断误差满足,则称方法(2.10)具有
9、p阶精度.,若将(2.10)展开式写成,则 称为局部截断误差主项.,以上定义对隐式单步法(2.9)也是适用的.例如,对后退欧拉法(2.5)其局部截断误差为,这里p=1是1阶方法,局部截断误差主项为,同样对梯形法(2.7)有,所以梯形方法(2.7)是二阶的. 其局部截断误差主项为,9.2.4 改进的欧拉公式,我们看到,梯形方法虽然提高了精度,但其算法复杂,在应用迭代公式(2.9)进行实际计算时,每迭代一次,都要重新计算函数f(x, y )的值,而迭代又要反复进行若干次,计算量很大,而且往往难以预测. 为了控制计算量,通常只迭代一两次就转入下一步的计算,这就简化了算法.,具体地说,我们先用欧拉公式
10、求得一个初步的近似值 ,称之为预测值,此预测值 的精度可能很差,再用梯形公式(2.7)将它校正一次,即按(2.8)式迭代一次,这个结果称之为校正值.,这样建立的预测校正系统通常称为改进的欧拉公式:,或表为下列平均化形式,(2.13),预测,校正,例2 用改进的欧拉法解例1中的初值问题(2.2).,解 仍取步长h=0.1,改进的欧拉公式为,部分计算结果见下表,同例1中的欧拉法的计算结果比较,明显改善了精度.,例 (两种方法的精度比较),用欧拉公式和改进的欧拉公式解下述初值问题,并与其准确解y=e-x+x进行比较.,解 取步长h=0.1,xk=kh(k=0,1,6).用两种方法进行计算对应结果及绝
11、对误差见下表,9.3 龙格库塔方法,对许多实际问题来说,欧拉公式与改进欧拉公式精度还不能满足要求,为此从另一个角度来分析这两个公式的特点,从而探索一条构造高精度方法的途径.,9.3.1 显式龙格库塔法的一般形式,上节给出了显式单步法的表达式(2.10), 其局部截断误差为O(hp+1),对欧拉法Tn+1=O(h2),即方法为p=1阶,若用改进欧拉法(2.13),它可表为,此时增量函数为,它比欧拉法的(xn, yn, h)=f(xn, yn), 增加了计算一个右函数f 的值,可望 p=2.若要使得到的公式阶数p更大, 就必须包含更多的f 值. 实际上从方程(1.1)等价的积分形式(2.4) ,即
12、,若要使公式阶数提高,就必须使右端积分的数值求积公式精度提高,它必然要增加求积节点,为此可将(3.3)的右端用求积公式表示为,一般说来,点数r 越多,精度越高,上式右端相当于增量函数(xn, yn, h),为得到便于计算的显式方法,可类似于改进欧拉法(3.1),(3.2),将公式表示为,其中,这里ci, i, ij均为常数. (3.4)和(3.5)称为r级显式龙格-库塔(Runge-Kutta)法, 简称R-K方法.,当r=1, (xn, yn, h)=f(xn, yn)时,就是欧拉法,此时方法的阶为p=1. 当r=2时,改进欧拉法(3.1),(3.2)是其中一种,下面将证明其阶p=2. 要使
13、公式(3.4),(3.5)具有更高的阶p,就要增加点数r. 下面我们只就r=2推导R-K方法. 并给出 r=3,4 时的常用公式,其推导方法与r=2时类似,只是计算较复杂.,9.3.2 二阶显式R-K方法,对r=2的R-K方法,由(3.4),(3.5)式可得如下计算公式,这里 c1, c2, 2, 21 均为待定常数,我们希望适当选取这些系数,使公式阶数 p 尽量高. 根据局部截断误差定义,推导出(3.6)的局部截断误差为,其中,这里yn=y(xn), yn+1=y(xn+1). 为得到Tn+1的阶p,要将上式各项在(xn, yn)处做泰勒展开,由于f(x, y )是二元函数,故要用二元泰勒展
14、开,各项展开式为,将以上结果代入(3.7),则有,要使公式(3.6)具有p=2阶,必须使,即,(3.9)的解是不唯一的. 可令c2=a0,则得,这样得到的公式称为二阶R-K方法.,则由此可以看出在改进的欧拉公式中相当于取(xn,yn), (xn+1,yn+1)两点处斜率的平均值,近似代替平均斜率,其精度比欧拉公式提高了.,如取a=1/2,则c1= c2=1/2, 2=21=1. 这就是改进的欧拉公式(3.1).,称为中点公式(变形的欧拉公式),相当于数值积分的中矩形公式.也可以表示为,如取a=1,则c1=0, c2=1, 2=21=1/2. 得计算公式,对r=2的R-K公式(3.6)能否使局部
15、误差提高到O(h4)? 为此 需把K2多展开一项,从(3.8)的 看到展开式中的项 是不能通过选择参数削掉的,实际上要使 h3 的项为零,需增加3个方程,要确定4个参数c1, c2, 2, 21,这是不可能的. 故r2的显式R-K方法的阶只能是p=2,而不能得到三阶公式.,9.3.3 三阶与四阶显式R-K方法,要得到三阶显式R-K方法,必须r=3. 此时计算(3.4), (3.5)的公式表示为,其中c1, c2, c3及2, 21, 3, 31, 32均为待定常数,公式(3.11)的局部截断误差为,只要K1, K2将按二元泰勒展开,使Tn+1O(h4),可得待定参数满足方程,这是8个未知数6个
16、方程的方程组,解不是唯一的. 可以得到很多公式. 满足条件(3.12)的公式(3.11)统称为三阶R-K公式. 下面只给出其中一个常见的公式.,此公式称为库塔三阶方法.,继续上述过程,经过较复杂的数学演算,可以导出各种四阶R-K公式,下列经典公式是其中常用的一个:,四阶R-K方法的每一步需要计算四次函数值f,可以证明其局部截断误差为O(h5).(例3见书p349),然而值得指出的是,龙格-库塔方法的推导基于泰勒展开方法,因而它要求所求的解具有较好的光滑性质. 反之,如果解的光滑性差,那么,使用龙格-库塔方法求得的数值解,其精度可能反而不如改进的欧拉方法. 实际计算时,我们应当针对问题的具体特点选择合适的算法.,9.3.4 变步长的龙格-库塔方法,单从每一步看,步长越小,截断误差就越小,但随着步长的缩小,在一定求解范围内所要完成的步数就增加了. 步数的增加不但引起计算量的增大,而且可能导致舍入误差的严重积累. 因此同积分的数值计算一样,微分方程的数值解法也有个选择步长的问题.,在选择步长时,需要考虑两个问题: 1.