北京科技大学计算方法课件10第十章高斯求积公式及常微分方程初值问题

1、北京科技大学数理学院北京科技大学数理学院 卫宏儒卫宏儒W科学与工程计算科学与工程计算高斯求积公式高斯求积公式 常见的四个正交多项式常见的四个正交多项式P213-217勒让德（勒让德（Legendre）多项式）多项式 切比雪夫（切比雪夫（Chebyshev）多项式）多项式 拉盖尔（拉盖尔（Laguerre）多项式）多项式 埃尔米特（埃尔米特（Hermite）多项式）多项式 常微分方程初值问题的数值解法常微分方程初值问题的数值解法 1 1、基本概念和定理：、基本概念和定理：一阶常微分方程初值问题是：一阶常微分方程初值问题是： y =f(x,y) (1.1) y(x0)=y0 (1.2) 其中其中f

2、是已知的是已知的xoy平面上某个区域平面上某个区域D上连续函数，上连续函数，式（式（1.1）是微分方程，有无穷多解，式（）是微分方程，有无穷多解，式（1.2）是）是确定解的初始条件。确定解的初始条件。 如一元函数如一元函数y(x)对一切对一切axb b 满足满足（1 1） （x,y(x)Dx,y(x)D（2 2） y(xy(x0 0)=y)=y0 0（3 3） y y 存在，且存在，且y y (x)=f(x,y(x)(x)=f(x,y(x)则称则称y(x)y(x)是初值问题（是初值问题（1.11.1）、（）、（1.21.2）在）在a,ba,b上上的解。的解。 关于初值问题解的存在、唯一及对初始

3、条关于初值问题解的存在、唯一及对初始条件的连续依赖性，有下列定理：件的连续依赖性，有下列定理： 定理定理1:1: 设设f(x,y)f(x,y)是在是在D=(x,y)|D=(x,y)| axb, b, cycy d 上的连续函数，其中上的连续函数，其中a a，b b为有限实为有限实数，而且数，而且f(x,y)f(x,y)满足对满足对y y的的lipschitzlipschitz条件，条件，则对则对(x(x0 0,y,y0 0) ) D,D,初值问题（初值问题（1.11.1），（），（1.21.2）在在a,ba,b的解存在且唯一。的解存在且唯一。 若若y(x)y(x)是式（是式（1.11.1），（

4、），（1.21.2）的解，从方程）的解，从方程（1.11.1） 两边积分，再利用式（两边积分，再利用式（1.21.2） 可得可得积分积分方程方程 反之，若反之，若y(x) y(x) 满足积分方程满足积分方程 （1.41.4），可验证），可验证它满足（它满足（1.11.1） 和（和（1.21.2），所以），所以 （1.41.4）式与）式与初值问题（初值问题（1.11.1），（），（1.21.2） 等价等价, ,这说明可用积这说明可用积分方程构造初值问题的数值解法。分方程构造初值问题的数值解法。 ).()(,()()(4100 xxdttytfxyxy定义定义3：若一种数值方法的：若一种数值方法的

5、局部截断误局部截断误差差O(hp+1)，则称相应数值方法是，则称相应数值方法是p阶方阶方法，其中法，其中p为正整数。为正整数。定义定义4：设：设y(x)是初始问题（是初始问题（1.1）的精）的精确解，确解，yn表示用某种数值方法算出的数表示用某种数值方法算出的数值解，值解， en= y(xn) - yn称为该方法在称为该方法在xn的的整体截断误差整体截断误差。 为了研究数值方法的绝对稳定为了研究数值方法的绝对稳定性，下面给出常系数线性差分方性，下面给出常系数线性差分方程的有关概念。程的有关概念。定义定义5 5： 方程方程111100n kkn knnyaya ya y 1110()0kknkr

6、ara ra r 1 12 2njnjnjnjk kyc rc rc r 0,1,jk 例例讨论线性多步法的绝对稳定性条件讨论线性多步法的绝对稳定性条件 得到相应的齐次线性差分方程：得到相应的齐次线性差分方程：00,0,0(,),.kkjkjjnjjjnnnjjyhfff xy 其中为常数不全为零将将应用于实验方程应用于实验方程yy 0()0kjjnjjhy 其对应的特征方程为：其对应的特征方程为：0()0kjjjjhr 0()0kjjnjjh 1knni iib r 2、数值解法的构造途径、数值解法的构造途径 （1）差商代替导数）差商代替导数 设初值问题设初值问题(1.1(1.1）的准确解）

7、的准确解y(x)y(x)在节在节点点x xn n之值为之值为y(xy(xn n),),记记y(xy(xn n) )的近似值为的近似值为y yn n ，又记，又记f fn n=f=f（x xn n ， y yn n ）, ,则初值问则初值问题题(1.1(1.1）离散化为）离散化为: : hxyhxyhxyhxyxykkkkhk)()()()(lim)( 0()()()(, ()(,)()()(,)0,1,2,.nnnnnnnnnnny xhy xy xf xy xf xyhy xhy xhf xyn即有得计算公式1(,)nnnnyyhf xy它称为（向前）它称为（向前）欧拉欧拉(Euler)公式

8、。（公式。（类类似地可以用向后差商、中心差商代替导似地可以用向后差商、中心差商代替导数产生相应的欧拉数产生相应的欧拉(Euler)公式公式） （2）数值积分法）数值积分法把把 y =f(x,y) 在在xn,xn+1积分，得积分，得对右端的定积分用数值积分方法做离对右端的定积分用数值积分方法做离散化，可得计算公式，如用散化，可得计算公式，如用矩形公式矩形公式可得欧拉公式，若用可得欧拉公式，若用梯形公式梯形公式可得改可得改进的欧拉公式，它也称为梯形公式：进的欧拉公式，它也称为梯形公式：11()()( , ( )nnxnnxy xy xf x y x dx111( (,)(,)2nnnnnnhyf

9、xyyf xy1111()()( , ( )( (, ()(, ()2nnxnnxnnnny xy xf x y x dxhf xy xf xy x （3）Taylor展开法 设设 f(x,y) f(x,y) 充分光滑充分光滑, ,将将y(xy(xn+1n+1) )在在x x n n点作点作TaylorTaylor展展开开: : y (x y (x n+1 n+1)=y(x)=y(xn n)+hy)+hy (x(xn n)+(h)+(h2 2/2!)y/2!)y”(x xn n)+O(h)+O(h3 3) )取其关于取其关于h h 的的线性部分线性部分，并用，并用y yn n 代替代替 y (

10、x y (xn n),),就得就得到到EulerEuler公式公式。易知易知EulerEuler公式的局部截断误差为公式的局部截断误差为 T T1 1= (h= (h2 2/2!)y/2!)y ”(x(xn n)+O(h)+O(h3 3) = O(h) = O(h2 2) )改进欧拉法的预改进欧拉法的预- -校公式校公式1111(,)(,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy11(,)(,)2nnnnnnnnyyhf xyf xyhf xyEuler公式的几何意义公式的几何意义abY=y(x)xy0例题：用例题：用EulerEuler公式和改进的公式和改进的Eule

11、rEuler公式分别公式分别求下列初值问题的数值解求下列初值问题的数值解( (取步长取步长h=0.1h=0.1计算计算到到y y3 3) )： y y =-2xy =-2xy2 2 y(0)=1y(0)=1解：由欧拉公式解：由欧拉公式 y y n+1 n+1=y=y n n+hf(x+hf(xn n,y,y n n)=y)=y n n - 2hx - 2hxn ny y n n2 2计算如下计算如下y y 1 1 = y = y 0 0 - 2hx - 2hx0 0y y 0 02 2 =1-20.101=1-20.1012 2 =1=1y y 2 2 = y = y 1 1 - 2hx -

12、2hx1 1y y 1 12 2 =1-20.10.11=1-20.10.112 2 =0.98=0.98y y 3 3= y= y2 2 - 2hx - 2hx2 2y y 2 22 2 =0.98-20.10.20.98=0.98-20.10.20.982 2 =0.9416=0.9416用改进欧拉法的预用改进欧拉法的预-校公式计算如下：校公式计算如下： 21221112(22)2nnnnnnnnnnyyhx yhyyx yxy计算如下计算如下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

13、)=0.9173 可见改进欧拉公式比欧拉公式精度高。可见改进欧拉公式比欧拉公式精度高。3、Runge-Kutta 方法 Runge-Kutta 方法是一种高精度的单步方法是一种高精度的单步法，简称法，简称R-K法。得到高精度方法的一个法。得到高精度方法的一个直接想法是利用直接想法是利用Taylor展开。展开。 假设式假设式 y =f(x,y) (axb) 中的中的 f(x,y) 充分光滑，将充分光滑，将y(xn+1)在在x n点作点作Taylor展开展开: （1）基本思想）基本思想对照标准形式对照标准形式 y n+1=y n+h(xn,y n,h) 。若若取取(x,y,h)=y(x)+(h/2

14、!)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-KR-K方法不是直接使用方法不是直接使用TaylorTaylor级数级数, ,而是利用它的思想而是利用它的思想, ,即计算即计算f(x,y)f(x,y)在不在不同结点的函数值同结点的函数值, ,然后作这些函数值的然后作这些函数值的线性组合线性组合, ,构造近似公式构造近似公式, ,式中有一些式中有一些可 供 选 择 的 参 数 。 将 近 似 公 式 与可 供 选 择 的 参 数 。

15、 将 近 似 公 式 与TaylorTaylor展开式相比较展开式相比较, ,使前面的若干项使前面的若干项密合密合, ,从而使近似公式达到一定的精度。从而使近似公式达到一定的精度。 下面以二级二阶下面以二级二阶R-KR-K方法为例说明方法为例说明这一方法的基本思想。这一方法的基本思想。 在在 xn , , xn+1 上上, ,取取f(x,y)f(x,y)在两个点的在两个点的函数值作线性组合函数值作线性组合, ,即得到二级即得到二级R-KR-K方法方法: : y y n+1=y=y n+h(c+h(c1 1K K1 1+c+c2 2K K2 2) ) K K1 1=f(=f(xn,y,y n)

16、) （* * *）K K2 2=f(=f(xn+a+a2 2h,yh,y n+b+b2121hKhK1 1) )其中其中c c1 1,c,c2 2,a,a2 2,b,b2121为待定参数。对照式为待定参数。对照式( (* *) )有：有： (x,y,h)=c1 1f(x,y)+c2 2f(x+a2 2h,y+b2121hf(x,y)（2）二级二阶）二级二阶R-K方法方法(xn,y(xn);h)=(c1 1+c2 2)y(xn)+c2 2(a2 2hfx+b2 21 1hfyf)+O(h2 2)因为因为y(xn+1)在在xn处的处的Taylor展开为展开为 y(xn+1)=y(xn)+hy(xn

17、)+(h2/2!)y(xn)+O(h3)由由显式单步法在显式单步法在xn+1的局部截断误差的局部截断误差定义有：定义有： Tn+1=y(y(xn+1 )-y( )-y(xn )-h )-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-1 K1=f(xn,yn) K2=f(xn+h/2,yn+hK1/2)这就

18、是变形的这就是变形的欧拉方法或中点方法欧拉方法或中点方法。 二级二级R-KR-K方法是显式单步式，每前进一步方法是显式单步式，每前进一步需要需要计算两个函数值计算两个函数值。由上面的讨论可知，由上面的讨论可知，适当选择四个参数适当选择四个参数c c1 1,c,c2 2,a,a2 2,b,b2121, ,可使每步计可使每步计算两次函数值的二阶算两次函数值的二阶R-KR-K方法达到二阶精度。方法达到二阶精度。能否在计算函数值次数不变的情况下能否在计算函数值次数不变的情况下, ,通过通过选择四个参数选择四个参数, ,使得二阶使得二阶R-KR-K方法的精度再提方法的精度再提高呢高呢? ? 答案是否定的

19、答案是否定的。无论四个参数怎样选择，。无论四个参数怎样选择，都不能使公式都不能使公式( (* * *) )提高到三阶。提高到三阶。 这说明每一步计算两个函数值的二阶这说明每一步计算两个函数值的二阶R-KR-K方法最高阶为二阶。若要获得更高阶得数值方法最高阶为二阶。若要获得更高阶得数值方法，就必须增加计算函数值的次数。方法，就必须增加计算函数值的次数。 仿照二级仿照二级R-KR-K方法，在方法，在 xn , , xn+1 上上, ,取取f f在在m m个个点的函数值做线性组合，即得到点的函数值做线性组合，即得到m m级级R-KR-K方法：方法：（3） m级显式级显式Runge-Kutta 方法方

20、法11111,2,3,mnnrrrnnrrnrnrsssyyhc KKfxyKfxhayhb Krm112341234123226,11,2211,22,nnnnnnnnnnhyyKKKKKfxyKfxh yhKfxh yhKfxh yKKhK 前面已经看到前面已经看到, ,二级、四级二级、四级R-KR-K方法可分别达到方法可分别达到最高阶数二阶、四阶，但是最高阶数二阶、四阶，但是N N级级R-KR-K方法的最高阶却方法的最高阶却不一定是不一定是N N阶阶。N N表示表示R-KR-K方法的级数表示公式中计方法的级数表示公式中计算函数值算函数值f f的次数。的次数。ButcherButcher给

21、出了给出了R-KR-K方法计算函方法计算函数值数值f f的次数与阶数之间的关系表，如下：的次数与阶数之间的关系表，如下：计算计算f f的次数的次数 1 2 3 1 2 3 4 4 5 6 7 5 6 7方法的最高阶数方法的最高阶数 1 21 2 3 3 4 44 54 5 6 6由表可见，四级以下由表可见，四级以下R-KR-K的方法其最高阶数与计算的方法其最高阶数与计算f f的次数一致，对的次数一致，对m m阶阶R-KR-K公式，当公式，当m4m4，虽然计算，虽然计算f f的的次数增加，但是方法阶数不一定增加。因此次数增加，但是方法阶数不一定增加。因此四级四四级四阶阶R-KR-K公式是应用最为广泛的公式。公式是应用最为广泛的公式。4 4、绝对稳定性问题、绝对稳定性问题(1)Euler方法的绝对稳定性 将将Euler方法应用于实验方程得到方法应用于实验方程得到：1(,)(1)nnnnnnnyyhf x yyhyh y 这是一个齐次线性差分方程，这是一个齐次线性差分方程，其对应的特征方程为：其对应的特征方程为： 10rh 由定理由定理2 2可知，当特征根可知，当特征根 1rh 满足满足 11rh 时时, , Euler方法是绝对