第6章--常微分方程数值解法-计算方法-《代码优化》课件

,第二级,第三级,第四级,第五级,第6章 常微分方程数值解法,计算方法,第6章 常微分方程数值解法,1 引言,2 欧拉法和改进的欧拉法,3 龙格库塔法,4 阿达姆斯方法,5 二阶线性常微分方程边值问题的数值解,1 引言,在常微分方程中，我们已经掌握了一些典型方程的解法。但许多形式的方程只能用数值方法求近似解，也就是求在某些点上满足一定精度的近似解。现以求一阶常微分方程初值问题,(61),在区间a，b上的解为例，介绍数值方法的基本思想。,设f(x，y)在带形区域,R：axb，-y+,上为x，y的连续函数，且对任意的y满足李普希茨(Libusize)条件,f(x，y,1,)-f(x，y,2,)Ly,1,-y,2, (62),其中(x，y,1,)、(x，y,2,)R，L为正常数。在求初值问题(61) 的数值解时，我们通常采用离散化方法(数值微分、数值积分、泰勒展式等)，求在自变量x的离散点,a=x,0,x,1,x,2,x,n,=b,图 6.1,2 欧拉法和改进的欧拉法,2.1 欧拉法(折线法),若将函数y(x)在点x,i,处的导数y(x,i,)用差商来表示，即,再用y,i,近似地代替y(x,i,)，则初值问题(61) 就化为,(63),式(63)就是所求的欧拉公式。,欧拉公式有很明显的几何意义。我们知道初值问题(61)中的微分方程的解是xoy平面上的一簇积分曲线，这簇积分曲线上任意点(x，y)?的斜率为f(x，y)，而初值问题(61)的解是过点(x,0,，y,0,)的一条特定的积分曲线。,例1 用欧拉法求初值问题,的数值解(取h=0.1)。,解 因为,故由欧拉计算公式(63)得,(64),图 6.2,图 6.3,代入(63)式得,(65),这样得到的点列仍为一折线，只是用平均斜率,来代替原来一点处的斜率。式(65)称为改进的欧拉,公式。,直到满足,(为预给精度),取,再转到下一步计算。,这里必须特别说明，因为初值问题(61)满足李普,希茨条件,当h足够小时，可使得,于是有,当k时，有q,k,0，故公式(66)收敛。,虽然式(67)仅迭代一次，但因进行了预先估计，故精度却有较大的提高。,在实际计算时，还常常将式(67)写成下列形式：,(68),2.4误差估计,初值问题(61)的等价积分方程为,(69),若对式(69)右端的积分采用各种不同的近似计算,方法，就可以得到初值问题(61)的各种不同的数值,解法。,例 如积分采用左矩形公式,图 6.4,用y,i,、y,i+1,分别代替y(x,i,)、y(x,i+1,)便得到欧拉公式(63)。,若积分采用梯形公式,在进行误差分析时，我们假设y,i,=y(x,i,)，考虑用,y,i+1,代替y(x,i+1,)而产生慕囟衔蟛睿康氖俏伺卸吓拉公式和改进的欧拉公式的精确度。,设初值问题(61)的准确解为y=y(x)，则利用泰勒公式,1. 欧拉公式的截断误差,由式(63)知,(611),对(69)式右端的积分采用梯形公式并根据梯形公式的误差可得到,(614),其中(x,i,，x,i+1,)，比较式(613)和(614)得,(615),因此,所以，改进的欧拉公式的截断误差为O(h3)，也即改进的欧拉法为二阶的。,可以验证，预估校正公式(67)与改进的欧拉公式的截断误差相同，均为O(h3)。这里略去证明。,例 2求解初值问题,表 62,3 龙格库塔法,3.1 泰勒级数展开法,我们还是假设y,i,=y(x,i,)，利用泰勒级数展开求y(x,i+1,)。式(610)就是y(x,i+1,)的泰勒展开式，若取右端前有限项作为y(x,i+1,)的近似值，就可得到计算y(x,i+1,)的各种不同截断误差的数值公式。,例 如，取前两项可得到,这里,y(x,i,)=f(x,i,，y(x,i,),y(x,i,)=fx(x,i,，y(x,i,)+fy(x,i,，y(x,i,)y(x,i,),=fx(x,i,，y(x,i,)+f(x,i,，y(x,i,)fy(x,i,，y(x,i,),类似地，若取前k项作为y(x,i+1,)的近似值，便得到截断误差为O(h,k,)的数值计算公式。这些公式的计算必须依赖于求y(x,i,)的k阶导数，除非f(x，y)足够简单，否则直接用泰勒展开法求解较为复杂。但是泰勒级数展开法的基本思想是许多数值方法的基础。,3.2 龙格库塔法,前面已经知道，初值问题(61)等价于,龙格库塔法的基本思想是：用f(x，y)在几个不同点的,数值加权平均来代替f(x,i,+,h,，y(x,i,+h)的值，而使截断,误差的阶尽可能高。,1. 二阶龙格库塔公式,将预估校正公式(68)改写成更一般的形式,(616),适当选取%、,1,、,2,%的值，使截断误差,y(x,i+1,)-y,i+1,的阶数尽可能高。这里仍假定y,i,=y(x,i,)，显然,2. 四阶龙格库塔公式,二阶龙格库塔公式是由使用在两个不同点上的函数值的线性组合而得到的。同样，我们用四个不同点上的函数值的线性组合就可得到四阶龙格库塔公式。设,y,i+1,=yi+h(,1,k,1,+,2,k,2,+,3,k,3,+,4,k,4,) (620),这里k,1,、k,2,、k,3,、k,4,为四个不同点上的函数值，分别设其为,k,1,=f(x,i,，y,i,),k,2,=f(x,i,+,1,h，y,i,+,11,k,1,h),k,3,=f(x,i,+,2,h，y,i,+,21,k,1,h+,22,k,2,h),k,4,=f(x,i,+,3,h，y,i,+,31,k,1,h+,32,k,2,h+,33,k,3,h),其中,1,、,2,、,3,、,4,、,1,、,2,、,3,、,11,、,21,、,22,、,31,、,32,、,33,均为待定系数。,类似于前面的讨论，把k,2,、k,3,、k,4,分别在x,i,点展成h的幂级数，代入线性组合式(620)中，将得到的公式与y(x,i+1,)在xi点上的泰勒展开式比较，使其两式右端直到h,4,的系数相等，经过较复杂的运算便可得到关于,i,，,i,，,ij,的一组特解,1,=,2,=,11,=,22,=1/2,21,=,31,=,32,=0,3,=,33,=1,1,=,4,=1/6,2,=,3,=1/3 (622),从而得到常用的标准四阶龙格库塔公式：,(623),图 6.5,表 63,4 阿达姆斯方法,我们已经知道，初值问题(61)等价于积分方程(69)，即,对积分式分别采用矩形公式和梯形公式可得到欧拉公式和改进的欧拉公式，截断误差分别为O(h,2,)和O(h,3,)。为此，我们自然可以想到，若用更高次的插值多项式来代替f(x，y)，则所得公式的精度会更高。这就是线性多步法的起源思想。,本章前面介绍的方法称为单步法，因为在计算y,i+1,时，只用到前面yi的值。而对于线性多步法是要利用前面已经算出的若干个值y,i-k,，y,i-1,，y,i,来求y,i+1,。,现用k次多项式P,k,(x)来代替f(x，y(x),(624),舍去余项,并设y,i,=y(x,i,)，而y,i+1,为y(x,i+1,)的近似值，于是可得到线性多步法的计算公式,4.1 阿达姆斯(Adams)显式,取q+1个基点x,i,，x,i-1,，x,i-q,，并作牛顿后差插值多项式见式(438)，则,其中,将式(627)代入式(626)得,(628),这里,(629),式(628)称为阿达姆斯显式。,对于余项,亦即,显然当q=3时,(630),m,是多项式积分，易算出结果如下：,m,0,1,2,3,4,m,1,1/2,5/12,3/8,251/270,例如,为了易于在电子计算机上实现，常将式(628)中的 用各点的已知函数值表示。特别，当q=2时，有,当q=3时，有,(631),(632),4.2 阿达姆斯隐式,类似于4.1，取q+1个基点x,i+1,，x,i,，x,i-q+1,，并作牛顿后差插值多项式，则,(633),其中,将式(633)代入式(626)得,(634),其中,(635),式(634)称为阿达姆斯隐式。,类似于阿达姆斯显式余项的求法，可得到阿达姆斯隐式的余项为,(636),例当q=3时,m,的计算结果如下：,m,0,1,2,3,4,m,1,-1/2,-1/12,-1/24,-19/720,若将式(634)中各阶差分mfi+1用各点的已知函数值表示，则可得到便于在电子计算机上实现的数值公式。,例如，当q=2时,(637),(638),当q=3时,4.3 阿达姆斯预估校正公式,我们常把阿达姆斯显式及隐式联立使用，即构造所谓阿达姆斯预估校正公式。现以q=2为例构造预估校正公式,(639),并取,与同阶的龙格库塔方法相比较，阿达姆斯方法计算量小，公式简单，程序易于实现。但它的主要缺点是不能自动开始，开始的前几个值要依赖于其它方法获得。这里介绍两种计算开始值的方法。,(1) 用单步法中的数值方法求出开始值。,(2) 使用y(x)的泰勒展开式,例4 用阿达姆斯方法求初值问题,(640),的数值解。,解 首先用泰勒展式求其三个点的值，因为,表 64,设常微分方程组的初值问题为,(641),这里,初值问题(641)与(61)式形式上完全相似。因此，对于(61)适用的数值计算公式，只要将其中的y,0,，y，f，都改写成相应的向量形式 s，y，f ，就能写出求解(641)的数值公式。,例如，初值问题(641)的标准四阶龙格库塔公式为,5二阶线性常微分方程边值问题的数值解,设二阶线性常微分方程的边值问题为,y+p(x)y+q(x)y=f(x),y(a)=，y(b)=，axb (642),其中p(x)，q(x)，f(x)为区间崐a，b上足够光滑的已知函数，且q(x)0，、为已知常数。,在上述条件下，边值问题(642)式存在连续可微的解，且是唯一的。若采用差分方法来解边值问题，其基本步骤是：,(1)将区间a，b“离散化”，即给a，b一个分划，此分划常考虑等距;,(2)对每一个基点，将各阶导数用差商来近似表示，将微分方程转化为差分方程，进而转化为线性代数方程组;,(3)解线性代数方程组，求得各基点上的近似解。,现具体给出求解边值问题(642)的方法步骤。,首先将区间a，b进行等距分划，即令,x,i,=a+ih，i=0，1，2，n,其中,一般称,为边界点，称x,1,，x,2,，x,n-1,为内,其次，在各基点x,i,上，将y，y用差商来近似表示。这里要求有相同阶数的截断误差，以保证精度协调。我们知道，由(463)式可得到,(643),又对式(458)再求一次导数，注意到,有,令t=1，可得到,取上式右边第一项作为pn(x,1,)的近似有,由差分与导数的关系可得余项为,所以，有,(644),略去(643)和(644)式的截断误差O(h,2,)，并用,y,i,代替y(x,i,)可有,(645),(646),将(645)、(646)代入(642)得到近似差分方程为,(647),其中,p,i,=p(x,i,)，q,i,=q(x,i,)，f,i,=f(x,i,),将式(447)整理后，得,其中,(649),这个方程组的系数矩阵是三对角的，可以利用追赶法求解。,根据(645)、(646)式知，对于y(x,i,)-y,i,的误差方程只要把fi取成O(h,2,)即可。现剩下两个问题：其一是线性方程组(648)即差分方程(647) 解的存在唯一性;其二是(648)的解的收敛性，也即当h0时，解是否收敛于微分方程(642)的准确解。,下面对差分方程(647)给出结论，不加以证明。,定理若a,i,0，c,i,0，-b,i,a,i,+c,i,，i=1，2，n-1，则二阶差分方程,a,i,y,i-1,+b,i,y,i,+c,i,y,i+1,=d,i,y,0,=，y,n,=，i=1，2，n-1,例5 用差分法解边值问题,y-y=x,y(0)=0，y(1)=1，0x1，h=01,解 这里步长h=1/10，则基点,表 65,