《第6章常微分方程数值解法 (2)讲义教材》由会员分享，可在线阅读，更多相关《第6章常微分方程数值解法 (2)讲义教材（70页珍藏版）》请在金锄头文库上搜索。
1、第6章 常微分方程数值解法,1 引言 2 欧拉法和改进的欧拉法 3 龙格-库塔法 4 阿当姆斯方法,1 引言,在高等数学里我们已经接触过常微分方程，对于一些典型的常微分方程，有求解析解的基本方法，但多数情况下，遇到的问题比较复杂，此时，只能利用近似方法求解，一般有两种近似方法 。,近似解析方法,数值方法,实际求解的常微分方程，大多是定解问题满足指定条件的特解,初值问题,边值问题,本章讨论常微分方程，数值解的最简单问题 一阶方程初值问题 ，即函数f(x)满足下列微分方程和初值条件： 在几何问题是(6-1)表现为一簇曲线，称(6-1)的积分曲线，初值问题(6-1) (6-2)就是要求一条过(x0
2、,y0)的积分曲线,(61),(62),定理1 对初值问题(6-1)(6-2)，若f(x,y)在区域 G= axb , |y| 内连续，且关于y满足李普希兹条件，即存在常数L，使 |f(x , y1)-f(x , y2)| L|y1-y2|（6-3） 对G中任意两个y1,y2均成立，其中L是与x,y无关的常数， 则初值问题(6-1)(6-2)在(a,b)内存在唯一解，且解是连续 可微的。,设f(x，y)在带形区域R：axb，-y+上为x，y的连续函数，且对任意的y满足李普希茨(Libusize)条件 f(x ，y1)-f(x ，y2)Ly1-y2 其中( x ，y1)、( x ，y2)R，L为
3、正常数。,在求初值问题(6-1)(6-2)的数值解时，我们通常采用离散化方法，求在自变量x的离散点a=x0 x1x2xn=b,上的准确解y(x)的近似值 y0，y1，y2，yn,常取离散点x0，x1，x2，xn为等距，即 x i+1-xi =h，i=0，1，2，n-1 h称为步长。图6.1表示为初值问题(61) (62)在n+1个离散点上的准确解y(x)的近似值。,图 6.1,数值解法的重点不在于求准确解（即解析解），而是直接求一系列点上的近似解。,求解过程顺着节点排列的顺序一步步向前推进，也即按递推公式由y0，y1.yi推yi+1，下面各种方法的实质是建立递推公式。,初值问题（6.1）(6.
4、2)的数值解法的基本特点是：,2 欧拉法和改进的欧拉法,一、欧拉方法 1. 基本思想 区间a，b上给定n+1个点x0，x1，x2，xn,再用yi近似地代替y(xi)，则初值问题(6-1)(6-2) 就化为,从x0出发根据初值问题，y(x0)=y0 再利用上式得 y(x1)y1=y(0)+hf(x0，y0)， 再以y1作为y(x1)的近似值，代入上式求y2.yn y(x2) y2=y1+hf(x1,y1) . y(xi+1) yi+1=yi+hf(xi,yi) i=0,1,(6.4) 称为解初值问题的欧拉方法,2几何意义,欧拉公式有很明显的几何意义。我们知道初值问题(6.1) 中的微分方程的解是
5、xoy平面上的一簇积分曲线 这簇积分曲线上任意点(x，y)的斜率为f(x，y)， 而初值问题(6.1) (6.2)的解是过点(x0，y0)的一条特定的积分曲线。,O,x0,x1,x2,xn,(x0,y0),(x1,y1),(x2,y2),(xn,yn),y(x)过点P0(x0,y0)，从P0出发以f(x0,y0)为斜率做一直线与直线x=x1交于点p1(x1,y1),显然有: y1=y0+hf(x0,y0), 再从p1出发,以f(x1,y1)为斜率做一直线推进到x=x2上一点p2(x2,y2),依此类推,这样得到解曲线的一条近似曲线,它就是折线 p0p1p2 所以欧拉方法又叫欧拉折线法,欧拉法是
6、用yi通过 yi+1=y i + hf(x i ,y i) i=0,1, 求yi+1,这样利用 y0y1y2 计算yi+1用前一步的y i单步法 计算yi+1用前几步的yi-n多步法,例1：用欧拉法求解方程,0 x1.2 h=0.2,解: 欧拉法的具体形式为: yi+1=yi+hf(xi, yi)=(1-0.4xi)yi 所以:,y1 = y0 + h f(x0,y0),y2 = y1 + h f(x1,y1) =(1-0.4x1)y1 =0.920000 ,=(1-0.4x0)y0,=1,可见欧拉法的精度是很差的,所求值用下表列出,并与精确值对比,二、欧拉方法的误差分析,定义1(p146)
7、对于初值问题，当假设yi是准确的时，用某种方法求yi+1时所产生的截断误差称为该方法的局部截断误差 。,我们来看在第i+1步使用欧拉方法所得yi+1的局部截断误差y(xi+1)-yi+1,假定yi是准确的,即yi=y(xi),由y(xi+1)=y(xi+h)，应用泰勒展开 y(xi+1)=y(xi+h) =y(xi)+hy(xi)+y()/2*h2 而由欧拉公式算出yi+1=yi+hf(xi,yi)= y(xi)+hf(xi,y(xi)=y(xi)+hy(xi),y(xi+1)=y(xi+h)= y(xi) + hy(xi) + y()/2*h2 yi+1=yi+hf(xi,yi)=y(xi)
8、 + hf(xi,y(xi) = y(xi) + hy(xi),两式相减得,y(xi+1)-yi+1=(h2/2)* y()=O(h2) 即欧拉方法所得yi+1的局部截断误差为O(h2),注意:只计算了一步，事实上每一步都有可能产生误差， 有时误差会原来越大，有时又会得到很好的控制，因此还 要考虑整体截断误差。,称为局部截断误差的主项,定义2 (p147),设yi是用某种方法计算初值问题(6-1)(6-2)在xi点的近似解，而y(xi)是它的精确解，则称 为该方法的整体截断误差，也称为该方法的精度。 补：若某方法的局部截断误差为O(hp+1)，则该方法的精度为p阶的。,欧拉方法的精度为O(h)
9、，一阶的,三、改进的欧拉法 欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当y=y(x)的曲线曲率较大时，欧拉法的效果更差。为了构造较高精度的数值解法，对初值问题 再做分析,对y=f（x，y）等式两边在（xi，xi+1）上取积分,xi,xi+1,(xi,yi),(xi+1,yi+1),在上图也可利用,矩形公式精度不高，再次说明欧拉精度低,这样就利用数值积分公式计算y(xi+1)的近似值。如果用左矩阵计算右面的积分：,为了得到更精确的方法我们可使用梯形公式,(6-9),改进的欧拉公式，又称为梯形公式,xi,xi+1,(xi,yi),(xi+1,yi+1),这样得到的点列仍为一折线，只是
10、用平均斜率来代替原来一点处的斜率。,不难发现: 欧拉公式 yi+1=yi +hf(x i,yi) 是关于yi+1的显式，只要已知yi，经一次计算可立即得到yi+1的值； 而改进欧拉公式,中的yi+1以隐式给出，且yi+1含在函数f(xi+1，yi+1)中， 所以梯形法是隐式单步法，,一般来说，这是一个非线性方程（除非f对y是线性的），可用我们前面讲过的非线性方程的各种方法求解，比如用迭代法,可证明其收敛，局部截断误差O(h3),k=0,1,2,具体做法是：先用欧拉公式(64)求出一个y(0)i+1作为初始近似值，然后再用改进的欧拉公式(69)进行迭代，即,直到满足,取,再转到下一步计算,例2:
11、 证明解常微分方程初值问题的梯形方法精度是二阶的 证明:,T=y(xi+1)-Th=y(x i+h)-T h,所以:,局部截断误差为O(h3 ) 所以精度为2阶的,四、预报-校正方法 我们看到梯形法虽然提高了精度，但其算法复杂，每算一点，都需进行反复迭代，为了控制计算量，通常只迭代一两次就转入下一步计算，以简化算法。 具体地说，我们先用欧拉公式求一个初步的近似值,称为预测值，其精度可能很差，再用梯形公式将其校正一次,预报-校正公式,(6-11),在实际计算时，还常常将式(611)写成下列形式：,(6-12),例3 用预报-校正公式求解初值问题,解：由预报-校正公式有,h=0.1, x0=0 ,
12、 x1=0.1 , x2=0.2,., x9=0.9 , x10=1,i=0,i=1,i=2,最后,得下表,其解析解为,表 62,预报-校正公式的截断误差为O(h3) （证明，p150) 预报-校正公式的整体截断误差为O(h2) 预报-校正公式是单步显式,练习1: 用Euler法解,取 h=0.2 计算到x=0.6,练习2: 用改进Euler法（梯形公式）解,取 h=0.2 小数点后至少保留5位,练习3 :用预报校正方法求,取步长 h=0.2 计算y(1.2)及y(1.4)的近似值，小数点后保留5位,3 龙格库塔法,一、基本思想 对初值问题,(axb),目的是求y(x)在一些给定点上的值y(x
13、1)，y(x2)，y(xn),研究差商,微分中值定理,得到一种微分方程的数值计算公式， 关键在求,称为在区间xi,xi+1上的平均斜率， 记为,K*,若用这种方法研究欧拉公式，可以发现欧拉公式仅取一个点(xi,yi)的斜率f(xi,yi)代替k*，梯形公式却是利用xi与xi+1两个点的斜率值,然后取值,由此也可知，梯形公式精度要高些。,它也启发我们若设法在xi,xi+1上多找几个点的斜率值，然后将它们加权平均作为K*的近似值，则有可能构造出精度更高的计算公式，这就是Runge-Kutta方法的基本思想。,即,于是得到一般公式,其中j，j，j为待定系数,(6-14),二、Runge-Kutta公
14、式,我们以两点xi，xi+p=xi+h为例，说明(6-14)式中系数的求法,其中1，2，为待定参数,(6-15),由二元函数Taylor展开式有,代入式(6-15)，得到yi+1的展开式,(6-16),利用,得到,把它与式(6-16)比较，使它们前三项相等， 得,则y(xi+1)-yi+1=O(h3),(6-17),满足(6-17)式的解有无穷多，对应于(6-15)的公式均具有二阶精度，统称为二阶龙格-库塔公式。,特别地，取 则(6-15)式正好是预报-校正公式。 如取,(6-18),则有,三阶龙格-库塔公式为,局部截断误差为O(h4),(6-19),四阶龙格-库塔公式为,(6-20),此公式
15、又称为经典的龙格-库塔公式，局部截断误差为O(h5),例 4 用经典的龙格-库塔法计算 取步长h=0.2,解：由x0=0，y0=1，h=0.2，利用公式(6-20)可计算出,从而,表6-3,补充-单步法的收敛性和稳定性 计算过程中舍入误差总是存在的，以Euler法为例。 假设，由于舍入误差的影响，实际得到的是.,两式相减,可得:,由微分中值定理得:,由此可知:只有当,误差不增,算法稳定,1. 单步法的收敛性 补充定义1: 设y(x)是初值问题(6.2)的精确解,对单步法,若有:,则称单步法收敛,补充定义2：用单步法解模型方程 得到解满足稳定性方程 若 ,就称此方法是绝对稳定的 在复平面上,所有
16、满足 所围成的区域称为方法的绝对稳定区域.,例5：求Euler法的绝对稳定区间,6.4 阿达姆斯方法,我们已经知道，初值问题(61)等价于积分方程，即,对积分式分别采用矩形和梯形面积公式 可得到欧拉公式和改进欧拉公式(梯形法)， 截断误差分别为O(h2)和O(h3)。 若追溯到数值积分原理，梯形法则 实际上是用节点 x i和x i+1 的线性插值函数代替f(x，y) 而得到的。 为此，我们自然可以想到，若用更高次的插值多项式来代替f(x，y)，则所得公式的精度会更高。 -这就是线性多步法的起源思想。,(624),对于线性多步法 是要利用前面已经算出的若干个值 yi-k，yi-1，yi来求yi+1。 现用k次多项式Pk(x)来代替f(x，y) 有: f (x ,y) =P k (x) + R k (x),舍去余项
