1、北京科技大学数理学院卫宏儒 W,科学与工程计算,高斯求积公式,常见的四个正交多项式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存在，
2、且y(x)=f(x,y(x) 则称y(x)是初值问题（1.1）、（1.2）在a,b上的解。,关于初值问题解的存在、唯一及对初始条件的连续依赖性，有下列定理：定理1: 设f(x,y)是在D=(x,y)| axb, cy d 上的连续函数，其中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）式与初值问题
3、（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),记
4、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)在x n点作Taylor展开:y (x n+1)=y(xn)+hy (xn)+(h2/2!)y”(xn)+O(h3)取其关于h
5、的线性部分，并用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 解：由欧拉公式y n+1=y n+hf(xn,y n)=y n - 2hxny n2 计算如下 y 1 = y 0 - 2hx0y 02 =1-20.1012 =1 y 2 = y 1 - 2hx1y 12 =1-20
6、.10.112 =0.98 y 3= y2 - 2hx2y 22 =0.98-20.10.20.982 =0.9416,用改进欧拉法的预-校公式计算如下：,计算如下 y 1 = 0.99; y 2 = 0.9614; y 3= 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)
7、在x n点作Taylor展开:,（1）基本思想,对照标准形式 y n+1=y n+h(xn,y n,h) 。若取(x,y,h)=y(x)+(h/2!)y(x)+(hp-1/p!)y(p)(x)并以y n代替y(xn)，则得到一个p阶近似公式y n+1=y n+h(xn ,y n ,h)(n=0,1,2,) (*),R-K方法不是直接使用Taylor级数,而是利用它的思想,即计算f(x,y)在不同结点的函数值,然后作这些函数值的线性组合,构造近似公式,式中有一些可供选择的参数。将近似公式与Taylor展开式相比较,使前面的若干项密合,从而使近似公式达到一定的精度。下面以二级二阶R-K方法为例说明
8、这一方法的基本思想。,在xn , xn+1 上,取f(x,y)在两个点的函数值作线性组合,即得到二级R-K方法:y n+1=y n+h(c1K1+c2K2)K1=f(xn,y n) （*） K2=f(xn+a2h,y n+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!
9、)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=1 c2a2=1/2 c2b21=1/2,当 =1时,c1=0,c2=1，得yn+1= yn+hK2 n=0,1,.N-1K1=f(xn,yn)K2=f(xn+h/2,yn+hK1/2) 这就是变形的欧拉方法或中点方法。,二级R-K方法是显式单步式，每前进一步需要计算两个函数值。由上面的讨论可知，
10、适当选择四个参数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方法的最高阶却不
11、一定是N阶。N表示R-K方法的级数表示公式中计算函数值f的次数。Butcher给出了R-K方法计算函数值f的次数与阶数之间的关系表，如下：计算f的次数 1 2 3 4 5 6 7 方法的最高阶数 1 2 3 4 4 5 6由表可见，四级以下R-K的方法其最高阶数与计算f的次数一致，对m阶R-K公式，当m4，虽然计算f的次数增加，但是方法阶数不一定增加。因此四级四阶R-K公式是应用最为广泛的公式。,4、绝对稳定性问题,(1)Euler方法的绝对稳定性,将Euler方法应用于实验方程得到：,这是一个齐次线性差分方程，其对应的特征方程为：,由上可知：,5、经典R-K法应用中步长的自动选取,注：一阶微分方程组与高阶方程的数值解法,（1）一阶微分方程组的解法前面介绍的单个方程的各种数值解法完全可以推广到一阶微分方程组的情形，只要将一阶微分方程组中的函数换为向量函数，得到方程组初值问题转化的类似于一阶微分方程的初值问题，按照前面的解法可给出欧拉格式或者经典的R-K方法。 （2）高阶方程的数值解法对于高阶微分方程，可通过引入变量代换，将其化为一阶方程组，再按（1）的方法求解。,例题：,因此得到：,作业10：,谢 谢！,
