数值分析讲义第四章 数值积分

1、第四章 数值积分与数值微分 /* Numerical Integration and differentiation*/,近似计算,1 引言, 对f()采用不同的近似计算方法，从而得到各种不同的求积公式。, 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有,求积节点,求积系数，与被积函数无关,像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。, 求积误差, 机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题： 精确度的度量标准； 如何构造具体的求积公式； 具体求积公式构造出来后，误差如何

2、估计？, 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有,求积节点,求积系数，与被积函数无关,像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。, 求积误差, 机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题： 精确度的度量标准； 如何构造具体的求积公式； 具体求积公式构造出来后，误差如何估计？, 以上三种方法都是用被积函数值的线性组合来表示积分值。推广，一般地有,求积节点,求积系数，与被积函数无关,像这样，将积分用若干节点上被积函数值的线性组合来表示的数值积分公式称为机械求积公式。

3、, 求积误差, 机械型求积公式的构造归结为，确定求积节点xk和求积系数Ak，使在某种意义下精确度较高。总之，要解决三个问题： 精确度的度量标准； 如何构造具体的求积公式； 具体求积公式构造出来后，误差如何估计？,代数精度与误差的关系：代数精度越高，求积误差越小。,问题1,由上面代数精度条件确定求积公式可分两种情形： 若事先给定求积节点xk(k=0,n),例如被积函数以表的形式给出时xk确定，可令m=n,由上式确定n+1个系数Ak即可-待定系数法和插值法。 若xk和Ak都可选择，令m=2n +1，确定xk和法Ak -Gauss法,要使求积公式具有m阶代数精度，则它对1,x,xm均准确成立，即,m

4、+1个方程，2n+2个未知数,问题2,Case 1-方法1,1 插值型求积 公式, 在a, b上取 a x0 x1 xn b，做 f 的 n 次插值多项式 ，即得到,节点,f (x),插值型积分公式 /*interpolatory quadrature*/,误差,Case 1-方法2,1 Newton-Cotes Formulae,梯形公式 /* trapezoidal rule*/,解：逐次检查公式是否精确成立,代入 P0 = 1：,=,代入 P1 = x ：,=,代入 P2 = x2 ：,代数精度 = 1,1 Newton-Cotes Formulae,Th1.形如 的求积公式至少有 n

5、次代数精度 该公式为插值型（即： ）, 当节点等距分布时：,令,Cotes系数,注：Cotes 系数仅取决于 n 和 i，可查表得到。与 f (x) 及区间a, b均无关。,2 Newton-Cotes 公式,NewtonCotes formula,1 Newton-Cotes Formulae,n = 1:,Trapezoidal Rule,/* 令 x = a+th, h = ba, 用中值定理 */,代数精度 = 1,n = 2:,Simpsons Rule,代数精度 = 3,n = 4: Cotes Rule, 代数精度 = 5,偶数阶N-C公式具有n+1阶代数精度,n = 3: Si

6、mpsons 3/8-Rule, 代数精度 = 3,对称节点的系数相同,Cotes公式是用不同节点的函数值（高度）的加权平均来近似区间的平均高度,注：当n8时，Cotes系数有负，造成公式不稳定，因此常用低阶Cotes公式。,证明：只需证明n为偶数时， N-C公式对f(x)=xn+1的余项R(f)=0即可。 因 f(n+1)(x)=(n+1)!, 由余项公式得,Th2. n为偶数时， N-C公式至少具有n+1阶代数精度。,注：当n 为偶数时，Cotes公式具有n+1阶精度，与n+1阶Cotes公式精度相同，但少计算一个节点上的函数值，因此一般常用偶数阶Cotes公式。,偶数阶N-C公式具有n+

7、1阶代数精度,N-C公式具有n阶代数精度余项R=o(h n+2),Hint：construct a interpolation polynomial of order 5, H(x), satisfying H(a)=f(a), H(b)=f(b), H (k)(a+b)/2) = f (k)(a+b)/2).,HW: p.151-152 #1-6,数值稳定性的一般概念,N-C的稳定性,3 复合求积 /* Composite Quadrature */,Havent we had enough formulae? Whats up now?,Oh come on, you dont serio

8、usly consider h=(ba)/2 acceptable, do you?,Why cant you simply refine the partition if you have to be so picky?,Dont you forget the oscillatory nature of high- degree polynomials!,Uh-oh,高次插值有Runge 现象，故采用分段低次插值 分段低次合成的 Newton-Cotes 复合求积公式。, 复合梯形公式：,在每个 上用梯形公式：,= Tn,/*中值定理*/,2 Composite Quadrature, 复化

9、 Simpson 公式：,= Sn,注：为方便编程，可采用另一记法：令 n = 2n 为偶数， 这时 ，有,2 Composite Quadrature, 收敛速度与误差估计：,定义,若一个复化积分公式的误差满足 且C 0，则称该公式是 p 阶收敛的。,/*中值定理*/,类似的，可得,2阶收敛,4阶收敛,6阶收敛,例1：计算,解：,其中,= 3.138988494,其中,= 3.141592502,Q: 给定精度 ，如何取 n ?,例如：要求 ，如何判断 n = ？,上例中若要求 ，则,即：取 n = 409,2 Composite Quadrature,4 龙贝格积分 /* Romberg

10、Integration */,复化梯形公式算法简单，但精度较差，收敛速度（2阶收敛）较慢，如何提高收敛速度？,注：按上面规律，可以构造线性组合系数为 的新的积分公式，但当m4时，前一个系数接近于1，后一个系数接近于0，这样构造出的新公式与前一个公式结果差别不大，反而增加计算量，因此实际上常做到Romberg公式为止。,例：计算,已知对于 = 106 须将区间对分 9 次，得到 T512 = 3.14159202,由 来计算 I 效果是否好些？,考察,= 3.141592502,= S4,一般有：,Romberg 序列, Romberg 算法：, ?, ?, ?, ,3 Romberg Inte

11、gration, 理查德森外推法 /* Richardsons extrapolation */,利用低阶公式产生高精度的结果。,设对于某一 h 0，有公式 T0(h) 近似计算某一未知值 I。由Taylor展开得到： T0(h) I = 1 h + 2 h2 + 3 h3 + ,i 与 h 无关,Q：如何将公式精度由 O(h) 提高到 O(h2) ？,即：,5 高斯型积分 /* Gaussian Quadrature */,构造具有2n+1次代数精度的插值型求积公式,5 高斯型积分 /* Gaussian Quadrature */,构造具有2n+1次代数精度的插值型求积公式,令 f (x)

12、 = 1, x, x2, , x2n+1 代入可求解，得到的公式具有2n+1 次代数精度。,例：求 的 2 点 Gauss 公式。,代入 f (x) = 1, x, x2, x3,不是线性方程组，不易求解。怎么办？先求出Gauss点，则问题转化为前面的问题，或者将方程组变为关于求积系数的线性方程组,将节点 x0 xn 以及系数 A0 An 都作为待定系数。 这样的节点称为Gauss 点，公式称为Gauss 型求积公式。,4 Gaussian Quadrature,证明： “”,x0 xn 为 Gauss 点, 则公式 至少有 2n+1 次代数精度。,对任意次数不大于n 的多项式 Pm(x)，

13、Pm(x) w(x)的次数不大于2n+1，则代入公式应精确成立：,= 0,“” 要证明 x0 xn 为 Gauss 点，即要证公式对任意次数不大于2n+1 的多项式 Pm(x) 精确成立，即证明：,设,求 Gauss 点 求w(x),4 Gaussian Quadrature, 正交多项式族 0, 1, , n, 有性质：任意次数不大于n 的多项式 P(x) 必与n+1 正交。,Step 1：构造正交多项式2,即：,4 Gaussian Quadrature,Step 2：求2 = 0 的 2 个根，即为 Gauss 点 x0 ，x1,Step 3：代入 f (x) = 1, x 以求解 A0

14、 ，A1,解线性方程组，简单。,结果与前一方法相同：, 利用此公式计算 的值,注：构造正交多项式也可以利用 Schmidt正交化方法进行（P.57）。,4 Gaussian Quadrature, 特殊正交多项式族：, Legendre 多项式族：,满足：,由 有递推,以 Pn+1 的根为节点的求积公式称为Gauss-Legendre 公式。, Chebyshev 多项式族：,注意到积分端点 1 可能是积分的奇点，用普通Newton-Cotes公式在端点会出问题。而Gauss公式可能避免此问题的发生。,其它公式见教材p.66,4 Gaussian Quadrature, Gauss 公式的余项

15、：,/* 设P为f 的过x0 xn的插值多项式 */,/*只要P 的阶数不大于2n+1，则下一步等式成立*/,插值多项式的余项,Q：什么样的插值多项式在 x0 xn 上有 2n+1 阶？,A：Hermite 多项式！,满足, Gauss 公式的稳定性：,求积系数的另一种计算方法,结论：Gauss型积分公式是数值稳定的。（稳定性分析类似于n 7的N-C公式）,区间-1,1上权函数(x)=1的Gauss型求积公式,称为Gauss-Legendre求积公式,其Gauss点为Legendre多项式的零点.,几种Gauss型求积公式,(1) Gauss-Legendre求积公式,公式的Gauss点和求积

16、系数可在数学用表中查到 .,用3点Gauss公式计算积分,解 查表得x1=-0.7745966692,x2=0,x3=0.7745966692, A1=A3=0.5555555556,A2=0.8888888889, 所以有,Gauss-Legendre求积公式的余项为,例13,误差为,实际上,I=2sin1=1.68294197, 误差为|R|=6.15810-5 .,用Simpson公式,则有I1.69353487, 误差为|R|=1.0610-2 .,由于,因此,a,b上权函数(x)=1的Gauss型求积公式为,用3点Gauss公式计算积分,结果远比Simpson公式的结果精确.,例14,解 这里Gauss点和积分系数与上例相同,所以,求积误差可表示为,区间0,)上权函数(x)=e-x的Gauss型求积公式,称为Gauss-Laguerre求积公式,其Gauss点为Laguerre多项式的零点.,(2) Gauss-Laguerre求积公式,公式的Gauss点和求积系数可在数学用表中查到 .,Gauss-Laguerre求积公式为,求积公式的误差为,由于,所以,对0, +)上权函数(x)=