《数值分析 常微分方程初值问题数值解法.pptx》由会员分享，可在线阅读，更多相关《数值分析 常微分方程初值问题数值解法.pptx（90页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、9.1 引引 言言 科学技术中常常需要求解常微分方程的定解问题.这类问题最简单的形式，是本章将着重考察的一阶方程的初值问题 我们知道，只有f(x,y)适当光滑譬如关于y满足利普希茨(Lipschitz)条件理论上就可以保证初值问题的解yf(x)存在并且唯一.第1页/共90页 虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.所谓数值解法,就是寻求解y(x)在一系列离散节点上的近似值 y1,y2,yn,yn+1,.相邻两个节点的间距hn=xn+1-xn称为步长.今后如不特别说明，总是假定 hi=h(i=1,2,)为定数,
2、这时节点为xn=x0+nh(i=0,1,2,)(等距节点).第2页/共90页 初值问题的数值解法有个基本特点，他们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.描述这类算法，只要给出用已知信息yn,yn-1,yn-2,计算yn+1的递推公式.首先，要对微分方程离散化，建立求解数值解的递推公式.一类是计算yn+1时只用到前一点的值yn，称为单步法.另一类是用到yn+1前面 k 点的值yn,yn-1,yn-k+1，称为k步法.其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性，还有递推公式的计算稳定性等问题.第3页/共90页9.2 简单的数值方法
3、与基本概念简单的数值方法与基本概念欧拉法与后退欧拉法 我们知道，在xy平面上，微分方程(1.1)式的解y=f(x)称作它的积分曲线，积分曲线上一点(x,y)的切线斜率等于函数f(x,y)的值.如果按f(x,y)在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，我们从初始点P0(x0,y0)出发,先依方向场在该点的方向推进到x=x1上一点P1，然后再从P1点依方向场在该点的方向推进到 x=x2 上一点P2,循环前进做出一条折线P0 P1 P2.第4页/共90页 一般地，设已做出该折线的顶点Pn,过Pn(xn,yn)依方向场的方向再推进到P
4、n+1(xn+1,yn+1)，显然两个顶点Pn,Pn+1的坐标有关系这就是著名的(显式)欧拉(Euler)公式.若初值y0已知，则依公式(2.1)可逐次逐步算出各点数值解.即第5页/共90页 例1 用欧拉公式求解初值问题 解 取步长h=0.1，欧拉公式的具体形式为其中xn=nh=0.1n(n=0,1,10),已知y0=1,由此式可得第6页/共90页依次计算下去，部分计算结果见下表.与准确解 相比，可看出欧拉公式的计算结果精度很差.xn 欧拉公式数值解欧拉公式数值解yn准确解准确解y(xn)误差误差 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1
5、.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719第7页/共90页 欧拉公式具有明显的几何意义,就是用折线近似代替方程的解曲线，因而常称公式(2.1)为欧拉折线法.还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.第8页/共90页
6、 为了分析计算公式的精度，通常可用泰勒展开将y(xn+1)在xn处展开，则有在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn)=y(xn).于是可得欧拉法(2.1)的公式误差为称为此方法的局部截断误差.第9页/共90页 如果对方程(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)近似，则得到另一个公式第10页/共90页 后退的欧拉公式与欧拉公式有着本质的区别,后者
7、是关于yn+1的一个直接计算公式，这类公式称作是显式的；前者公式的右端含有未知的yn+1，它实际上是关于yn+1的一个函数方程,这类方程称作是隐式的.显式与隐式两类方法各有特点，考了到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便.隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化.第11页/共90页 设用欧拉公式给出迭代初值 ，用它代入(2.5)式的右端，使之转化为显式，直接计算得然后再用 代入(2.5)式，又有如此反复进行，得第12页/共90页由于f(x,y)对y满足Lipschitz条件(1.3).由(2.6)减(2.5)得由此可知，只要hL,我们反复将步
8、长折半计算,直至为止,这时取最终得到的 作为结果；2.如果为止，这时再将步长折半计算一次，就得到所要的结果.第46页/共90页9.4 单步法的收敛性与稳定性单步法的收敛性与稳定性收敛性与相容性 数值解法的基本思想是，通过某种离散化手段将微分方程转化为差分方程，如单步法 它在点xn处的解为yn，而初值问题在点xn处的精确解为y(xn)，记en=y(xn)-yn称为整体截断误差.收敛性就是讨论当 x=xn 固定且 时en0的问题.第47页/共90页 定义3 若一种数值方法对于固定的xn=x0+nh,当h0时有yny(xn)，其中y(x)是(1.1),(1.2)的准确解，则称该方法是收敛的.显然数值
9、方法收敛是指en=y(xn)-yn0，对单步法(4.1)有下述收敛性定理：定理1 假设单步法(4.1)具有p阶精度，且增量函数(x,y,h)关于y满足利普希次条件又设初值y0是准确的,即y0=f(x0),则其整体截断误差第48页/共90页 证明 设以yn+1表示取yn=y(xn)用公式(4.1)求得的结果，即 则y(xn)-yn+1为局部截断误差，由于所给方法具有p阶精度，按定义2，存在定数C，使又由式(4.4)与(4.1)，得第49页/共90页利用利普希次条件(4.2)，有从而有即对整体截断误差en=y(xn)-yn成立下列递推关系式据此不等式反复递推，可得第50页/共90页由此可以断定，判
10、断初值是准确的，即e0=0，则(4.3)式成立.定理证毕.再注意到当x=x0+nhT时最终得下列估计式依据这一定理，判断单步法(4.1)的收敛性，归结为验证增量函数能否满足利普希次条件(4.2).对于欧拉方法，由于其增量函数 就是f(x,y),故当f(x,y)关于y满足利普希次条件时它是收敛的.第51页/共90页再考察改进的欧拉方法，其增量函数已由(3.2)式给出，这时有设限定hh0(h0为定数)，上式表明关于y的利普希次常数为因此改进的欧拉方法也是收敛的.类似地,不难验证其它龙格-库塔方法的收敛性.第52页/共90页 定理1表明p1时单步法收敛,并且当y(x)是初值问题(1.1),(1.2)
11、的解,(4.1)具有p阶精度时,则有展开式所以p1的充分必要条件是 ，而 ，于是可给出如下定义：第53页/共90页 定义4 若单步法(4.1)的增量函数满足 以上讨论表明p阶方法(4.1)当p1时与(1.1),(1.2)相容，反之相容方法至少是1阶的.于是由定理1可知方法(4.1)收敛的充分必要条件是此方法是相容的.则称单步法(4.1)与初值问题(1.1),(1.2)相容.第54页/共90页绝对稳定性与绝对稳定域绝对稳定性与绝对稳定域 前面关于收敛性的讨论有个前提，必须假定数值方法本身的计算是准确的.实际情形并不是这样，差分方程的求解还会有计算误差.譬如由于数字舍入而引起的小扰动.这类小扰动在
12、传播过程中会不会恶性增长，以至于“淹没”了差分方程的“真解”呢？这就是差分方程的稳定性问题.在实际计算时，我们希望某一步产生的扰动值，在后面的计算中能够被控制，甚至是逐步衰减的.第55页/共90页 定义5 若一种数值方法在节点值yn上大小为的扰动，于以后各节点值ym(mn)上产生的偏差均不超过，则称该方法是稳定的.下面以欧拉法为例考察计算稳定性.例4 用欧拉公式求解初值问题 解 用欧拉法解方程y=-100y 得其准确解 是一个按指数曲线衰减很快的函数.第56页/共90页若取步长h=0.025，则欧拉公式的具体形式为节点节点xn欧拉方法欧拉方法yn后退欧拉方法后退欧拉方法yn0.0250.050
13、0.0750.100-1.5 2.25 -3.375 5.06250.28570.08160.02330.0067计算结果见表,明显计算过程不稳定,但取h=0.005,yn+1=-1.5yn,则计算过程稳定.对后退的欧拉公式，取h=0.025时，则计算公式为yn+1=-(1/3.5)yn.计算结果见表,这时计算过程是稳定的.第57页/共90页 例题表明稳定性不但与方法有关，也与步长h有关，当然与方程中的f(x,y)有关.为了只考察数值方法本身，通常只检验数值方法用于解模型方程的稳定性，模型方程为其中为复数，这个方程分析较简单，对一般方程可以通过局部线性化化为这种形式，例如在(x,y)的邻域，可
14、展开为第58页/共90页略去高阶项，再做变换即可得到的形式 u=u.对于m个方程的方程组,可线性化为y=Ay,这里A为mm雅可比矩阵(fi/yj)，若A有m个特征值1,2,m，其中i可能是复数，所以，为了使模型方程结果能推广到方程组，方程(4.8)中为复数.为保证微分方程本身的稳定性，还应假定Re()0.下面先研究欧拉方法的稳定性.模型方程y=y的欧拉公式为第59页/共90页设在节点yn上有一扰动值n，它的传播使节点值yn+1产生大小为的扰动值n+1，假设用yn*=yn+n，按欧拉公式得出yn+1*=yn+1+n+1的计算过程不再有新的误差，则扰动值满足可见扰动值满足原来的差分方程(4.9).
15、这样，如果差分方程的解是不增长的，即有则它就是稳定的.这一论断对于下面将要研究的其它方法同样适用.第60页/共90页显然，为要保证差分方程(4.9)的解是不增长的，只要选取h充分小，使在=h的复平面上，这是以(-1,0)为圆心，1为半径的单位圆.称为欧拉法的绝对稳定域，一般情形可由下面定义.定义6 单步法(4.1)用于解模型方程y=y，若得到的解yn+1=E(h)yn，满足|E(h)|1，则称方法(4.1)是绝对稳定的.在=h的平面上,使|E(h)|1的变量围成的区域，称为绝对稳定区域，它与实轴的交称为绝对稳定区间.第61页/共90页对欧拉法E(h)=1+h，其绝对稳定域为|1+h|1,绝对稳
16、定区间为-20,在例5中=-100,-2-100h0,即0h2/100=0.02为稳定区间，在例4中取h=0.025，故它是不稳定的，当取h=0.005时它是稳定的.对二阶R-K方法，解模型方程(4.1)可得到故绝对稳定域由|E(h)|1得到，于是可得绝对稳定区间为-2h0，即0h2/.第62页/共90页 类似可得三阶及四阶R-K方法的E(h)分别为由|1+h|1可得到相应的绝对稳定域.当为实数时，则得绝对稳定区间，它们分别为 三阶显式R-K方法：四阶显式R-K方法：从以上讨论可知显式R-K方法的绝对稳定域均为有限域，都对步长h有限制.如果h不在所给的绝对稳定区间内，方法就不稳定.第63页/共
17、90页 例4 分别取h=0.1及h=0.2，用经典的四阶R-K方法(3.1)计算初值问题 解 本例=-20,h分别为-2及-4，前者在绝对稳定区间内，后者则不在，用四阶R-K方法计算其误差见下表xn0.20.40.60.81.0h=0.1h=0.20.9310-14.980.1210-125.00.1410-2125.00.1510-3625.00.1710-43125.0从以上结果看到，如果步长h不满足绝对稳定条件，误差增长很快.第64页/共90页对隐式单步法,可以同样讨论方法的绝对稳定性,例如对后退欧拉法，用他解模型方程可得故由|E(h)|1，这是以(1,0)为圆心，1为半径的单位圆外部.
18、故方法的绝对稳定区间为-h0.当0时，则0h，即对任何步长均为稳定的.第65页/共90页 对隐式梯形法，它用于解模型方程(4.8)得故对Re()0有|E(h)|1，故绝对稳定域为=h的左半平面，绝对稳定区间为-h0，即0h时隐式梯形法均是稳定的.第66页/共90页9.5 线性多步法线性多步法 在逐步推进的求解过程中，计算yn+1之前事实上已经求出了一系列的近似值y0,y1,yn，如果充分利用前面多步的信息来预测yn+1，则可以期望会获得较高的精度.这就是构造所得线性多步法的基本思想.构造多步法的主要途径基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1.1)两端积分后利用插值求积公式得到
19、.本节主要介绍基于泰勒展开的构造方法.第67页/共90页线性多步法的一般公式线性多步法的一般公式 如果计算yn+k时，除用yn+k-1的值，还要用到yn+i(i=0,1,k-2)的值，则称此方法为线性多步法.一般的线性多步法公式可表示为 其中yn+1为y(xn+1)的近似，fn+i=f(xn+i,yn+i),这里xn+i=xn+ih，i,i为常数，0及0不全为零，则称(5.1)为线性k步法，计算时需先给出前面k个近似值y0,y1,yk-1，再由(5.1)逐次求出yk,yk+1,.第68页/共90页 如果k=0，则(5.1)称为显式k步法，这时yn+k可直接由(5.1)算出；如果k0,则(5.1
20、)称为隐式k步法,求解时与梯形法(2.7)相同,要用迭代法方可算出yn+k.(5.1)中系数i及i可根据方法的局部截断误差及阶确定，其定义为 定义7 设y(x)是初值问题(1.1),(1.2)的准确解，线性多步法(5.1)在xn+k上局部截断误差为若Tn+k=O(hp+1)，则称方法(5.1)是p阶的，p1则称方法(5.1)与方程(1.1)是相容的.第69页/共90页由定义7，对Tn+k在xn处泰勒展开，由于代入(5.2)得其中第70页/共90页若在公式(5.1)中选择系数i及i，使它满足由定义可知此时所构造的多步法是p阶的，且称右端第一项为局部截断误差主项,cp+1称为误差常数.第71页/共
21、90页根据相容性定义,p1，即c0=c1=0，由(5.4)得故线性多步法(5.1)与微分方程(1.1)相容的充分必要条件是(5.6)成立.显然，当k=1时，若1=0，则由(5.6)可求得0=1，0=1.此时公式(5.1)为第72页/共90页即为欧拉法.从(5.4)可求得c2=1/20，故方法为1阶精度，且局部截断误差为这和第2节给出的定义及结果是一致的.对k=1，若10，此时方法为隐式公式，为了确定系数0,0,1，可由c0=c1=c2=0解得0=1,0=1=1/2.于是得到公式即为梯形公式.第73页/共90页由(5.4)可求得c2=-1/12，故p=2，所以梯形法是二阶方法，其局部截断误差主项
22、是这与第2节中讨论是一致的.对k2的多步法公式都可利用(5.4)确定系数i,i,并由(5.5)给出局部截断误差，下面只就若干常用的多步法导出具体公式.第74页/共90页阿当姆斯显式与隐式公式阿当姆斯显式与隐式公式p362自学.米尔尼方法与辛普森方法p366自学.汉明方法p367自学.预测-校正方法p368自学.构造多步法公式的注记和例p371自学.第75页/共90页9.6 方程组和高阶方程方程组和高阶方程一阶方程组 前面我们研究了单个方程y=f 的数值解，只要把y 和f 理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形.考察一阶方程组的初值问题，初始条件给为第76页/共90页若
23、采用向量的记号，记(向量)则上述方程组的初值问题可表示为第77页/共90页求解这一初值问题的四阶龙格-库塔公式为(向量)式中(向量)第78页/共90页或表示为(分量)其中(分量)这里yin是第i个因变量yi(x)在节点xn=x0+nh的近似值.第79页/共90页为了帮助理解这一公式的计算过程，我们再考察两个方程的特殊情形这时四阶龙格-库塔公式具有形式第80页/共90页其中第81页/共90页这是一步法，利用节点xn上的值yn,zn，由(6.3)式顺序计算K1,L1,K2,L2,K3,L3,K4,L4，然后代入(6.2)式即可求得节点xn+1上的yn+1,zn+1.第82页/共90页化高阶方程为一
24、阶方程组化高阶方程为一阶方程组关于高阶微分方程(或方程组)的初值问题，原则上总可以归结为一阶方程组来求解.例如，考察下列m阶微分方程初始条件为只要引进新的变量即可将m阶方程(6.4)化为如下的一阶方程组：第83页/共90页初始条件(6.5)则相应地化为不难证明初始条件(6.4),(6.5)和(6.6),(6.7)是彼此等价的.第84页/共90页特别地，对于下列二阶方程的问题引进新变量z=y,即可化为下列一阶方程组的初值问题:第85页/共90页针对这个问题应用四阶龙格-库塔公式(6.2)，有由(6.3)式可得第86页/共90页如果消去K1,K2,K3,K4，则上述格式可表示为这里第87页/共90页刚性方程组刚性方程组p378自学.第88页/共90页课程结束!希望同学们好好复习!争取考个好成绩!再见!第89页/共90页感谢您的观看！第90页/共90页