高斯求积公式及常微分方程初值问题

北京科技大学数理学院卫宏儒Weihr168,科学与工程计算,高斯求积公式,常见的四个正交多项式P213-217,勒让德（Legendre）多项式,切比雪夫（Chebyshev）多项式,拉盖尔（Laguerre）多项式,埃尔米特（Hermite）多项式,常微分方程初值问题的数值解法,1、基本概念和定理：一阶常微分方程初值问题是：y=f(x,y)(1.1)y(x0)=y0(1.2)其中f是已知的xoy平面上某个区域D上连续函数，式（1.1）是微分方程，有无穷多解，式（1.2）是确定解的初始条件。如一元函数y(x)对一切axb满足（1）（x,y(x)D（2）y(x0)=y0（3）y存在，且y(x)=f(x,y(x)则称y(x)是初值问题（1.1）、（1.2）在a,b上的解。,关于初值问题解的存在、唯一及对初始条件的连续依赖性，有下列定理：定理1:设f(x,y)是在D=(x,y)|axb,cyd上的连续函数，其中a，b为有限实数，而且f(x,y)满足对y的lipschitz条件，则对(x0,y0)D,初值问题（1.1），（1.2）在a,b的解存在且唯一。,若y(x)是式（1.1），（1.2）的解，从方程（1.1）两边积分，再利用式（1.2）可得积分方程反之，若y(x)满足积分方程（1.4），可验证它满足（1.1）和（1.2），所以（1.4）式与初值问题（1.1），（1.2）等价,这说明可用积分方程构造初值问题的数值解法。,定义3：若一种数值方法的局部截断误差O(hp+1)，则称相应数值方法是p阶方法，其中p为正整数。定义4：设y(x)是初始问题（1.1）的精确解，yn表示用某种数值方法算出的数值解，en=y(xn)-yn称为该方法在xn的整体截断误差。,为了研究数值方法的绝对稳定性，下面给出常系数线性差分方程的有关概念。,定义5：,方程,例,讨论线性多步法的绝对稳定性条件,得到相应的齐次线性差分方程：,其对应的特征方程为：,2、数值解法的构造途径,（1）差商代替导数设初值问题(1.1）的准确解y(x)在节点xn之值为y(xn),记y(xn)的近似值为yn，又记fn=f（xn，yn）,则初值问题(1.1）离散化为:,它称为（向前）欧拉(Euler)公式。（类似地可以用向后差商、中心差商代替导数产生相应的欧拉(Euler)公式）,（2）数值积分法把y=f(x,y)在xn,xn+1积分，得对右端的定积分用数值积分方法做离散化，可得计算公式，如用矩形公式可得欧拉公式，若用梯形公式可得改进的欧拉公式，它也称为梯形公式：,（3）Taylor展开法设f(x,y)充分光滑,将y(xn+1)在xn点作Taylor展开:y(xn+1)=y(xn)+hy(xn)+(h2/2!)y”(xn)+O(h3)取其关于h的线性部分，并用yn代替y(xn),就得到Euler公式。易知Euler公式的局部截断误差为T1=(h2/2!)y”(xn)+O(h3)=O(h2)改进欧拉法的预-校公式,Euler公式的几何意义,a,b,Y=y(x),x,y,0,例题：用Euler公式和改进的Euler公式分别求下列初值问题的数值解(取步长h=0.1计算到y3)：y=-2xy2y(0)=1解：由欧拉公式yn+1=yn+hf(xn,yn)=yn-2hxnyn2计算如下y1=y0-2hx0y02=1-20.1012=1y2=y1-2hx1y12=1-20.10.112=0.98y3=y2-2hx2y22=0.98-20.10.20.982=0.9416,用改进欧拉法的预-校公式计算如下：,计算如下y1=0.99;y2=0.9614;y3=0.9173精确解y(0.1)=0.99,y(0.2)=0.9614;y(0.3)=0.9173可见改进欧拉公式比欧拉公式精度高。,3、Runge-Kutta方法,Runge-Kutta方法是一种高精度的单步法，简称R-K法。得到高精度方法的一个直接想法是利用Taylor展开。假设式y=f(x,y)(axb)中的f(x,y)充分光滑，将y(xn+1)在xn点作Taylor展开:,（1）基本思想,对照标准形式yn+1=yn+h(xn,yn,h)。若取(x,y,h)=y(x)+(h/2!)y(x)+.+(hp-1/p!)y(p)(x)并以yn代替y(xn)，则得到一个p阶近似公式yn+1=yn+h(xn,yn,h)(n=0,1,2,.)(*),R-K方法不是直接使用Taylor级数,而是利用它的思想,即计算f(x,y)在不同结点的函数值,然后作这些函数值的线性组合,构造近似公式,式中有一些可供选择的参数。将近似公式与Taylor展开式相比较,使前面的若干项密合,从而使近似公式达到一定的精度。下面以二级二阶R-K方法为例说明这一方法的基本思想。,在xn,xn+1上,取f(x,y)在两个点的函数值作线性组合,即得到二级R-K方法:yn+1=yn+h(c1K1+c2K2)K1=f(xn,yn)（*）K2=f(xn+a2h,yn+b21hK1)其中c1,c2,a2,b21为待定参数。对照式(*)有：(x,y,h)=c1f(x,y)+c2f(x+a2h,y+b21hf(x,y),（2）二级二阶R-K方法,(xn,y(xn);h)=(c1+c2)y(xn)+c2(a2hfx+b21hfyf)+O(h2)因为y(xn+1)在xn处的Taylor展开为y(xn+1)=y(xn)+hy(xn)+(h2/2!)y(xn)+O(h3)由显式单步法在xn+1的局部截断误差定义有：Tn+1=y(xn+1)-y(xn)-h(xn,y(xn),h)=h(1-c1-c2)y(xn)+h2(1/2-a2c2)fx+(1/2-c2b21)fyf+O(h3)显然，若要求Tn+1=O(h3),则应有c1+c2=1c2a2=1/2c2b21=1/2,当=1时,c1=0,c2=1，得yn+1=yn+hK2n=0,1,.N-1K1=f(xn,yn)K2=f(xn+h/2,yn+hK1/2)这就是变形的欧拉方法或中点方法。,二级R-K方法是显式单步式，每前进一步需要计算两个函数值。由上面的讨论可知，适当选择四个参数c1,c2,a2,b21,可使每步计算两次函数值的二阶R-K方法达到二阶精度。能否在计算函数值次数不变的情况下,通过选择四个参数,使得二阶R-K方法的精度再提高呢?答案是否定的。无论四个参数怎样选择，都不能使公式(*)提高到三阶。这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶。若要获得更高阶得数值方法，就必须增加计算函数值的次数。,仿照二级R-K方法，在xn,xn+1上,取f在m个点的函数值做线性组合，即得到m级R-K方法：,（3）m级显式Runge-Kutta方法,前面已经看到,二级、四级R-K方法可分别达到最高阶数二阶、四阶，但是N级R-K方法的最高阶却不一定是N阶。N表示R-K方法的级数表示公式中计算函数值f的次数。Butcher给出了R-K方法计算函数值f的次数与阶数之间的关系表，如下：计算f的次数1234567方法的最高阶数1234456由表可见，四级以下R-K的方法其最高阶数与计算f的次数一致，对m阶R-K公式，当m4，虽然计算f的次数增加，但是方法阶数不一定增加。因此四级四阶R-K公式是应用最为广泛的公式。,4、绝对稳定性问题,(1)Euler方法的绝对稳定性,将Euler方法应用于实验方程得到：,这是一个齐次线性差分方程，其对应的特征方程为：,由上可知：,5、经典R-K法应用中步长的自动选取,注：一阶微分方程组与高阶方程的数值解法,（1）一阶微分方程组的解法前面介绍的单个方程的各种数