关于数值积分的几点教学注记

Several Teaching Notes on Numerical Integration

摘要: 本文对常用的数值积分公式给出了相关教学注记。首先，我们引入Newton-Cotes公式，并利用该公式证明了例1中定积分的数值估计。此例说明数值积分公式在数学理论上的有趣应用，为此，我们给出教学注记。然后，我们引入Gauss求积公式，通过例2比较了Newton-Cotes公式，Gauss求积公式和Clenshaw-Curtis公式的数值模拟过程，并且基于数值积分公式的代数精度，数值稳定性和收敛性给出了相关详细的教学注记。

Abstract: This paper gives relevant teaching notes on commonly used numerical integration formulas. First, we introduce the Newton-Cotes formula and use it to prove the numerical estimation of the definite integral in Example 1. This example illustrates the interesting application of the numerical integration formula in mathematical theory. To this end, we give teaching notes. Then, we introduce the Gauss quadrature formula and compare the numerical simulation processes of the Newton-Cotes formula, the Gauss quadrature formula and the Clenshaw-Curtis formula through Example 2. Relevant detailed teaching notes are given based on the algebraic precision, numerical stability and convergence of the numerical integration formula.

1. 引言

本文主要讨论如下定积分

$I\left(f\right)={\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\left(b>a\right)$ ，

由微积分知识可得，若可求得$f\left(x\right)$ 的原函数$F\left(x\right)$ ，则$I\left(f\right)=F\left(b\right)-F\left(a\right)$ ，该公式称为Newton-Leibniz公式。然而在实际应用中，$f\left(x\right)$ 的原函数可能不是初等函数，从而该公式无法直接应用，例如${\displaystyle {\int }_{-1}^1{\text{e}}^{-x^2}\text{d}x}$ 就无法精确计算。数值积分公式为近似逼近$I\left(f\right)$ 提供了一个较好的计算手段。数值积分的主要思想是通过如下机械求积公式进行近似计算[1]：

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}\approx {\displaystyle \sum _{k=0}^n{w}_{k}f\left(x_k\right)}$ ，

其中$x_k$ 为求积节点，$w_k$ 为求积系数。

2. Newton-Cotes公式

Newton-Cotes公式为如下形式[1]：

$I_N^n\left(f\right)={\displaystyle \sum _{k=0}^n{w}_{k}f\left(x_k\right)}\approx {\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ ，

其中$w_k=\left(b-a\right)C_k^{\left(n\right)}$ ，$C_k^{\left(n\right)}=\frac{{\left(-1\right)}^{n-k}}{nk!\left(n-k\right)!}{\displaystyle {\int }_0^n{\displaystyle \prod _{\begin{array}{l}j=0\\ j\ne k\end{array}}^n\left(t-j\right)\text{d}t}}$ ，$x_k=a+kh$ ，$h=\frac{b-a}{n}$ 。当$n=1$ 时，即为梯形公式

$I_N^1\left(f\right)=\frac{b-a}{2}\left[f\left(a\right)+f\left(b\right)\right]$ ， (1)

相应，我们有如下误差式

$I\left(f\right)-I_N^1\left(f\right)=\frac{-{\left(b-a\right)}^3}{12}{f}^{″}\left(\theta \right),\theta \in \left(a,b\right)$ 。 (2)

此外，中矩形公式为

$I_D\left(f\right)=\left(b-a\right)f\left(\frac{a+b}{2}\right)\approx {\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ ， (3)

相应的误差式为

$I\left(f\right)-I_D\left(f\right)=\frac{{\left(b-a\right)}^3}{24}{f}^{″}\left(\theta \right),\theta \in \left(a,b\right)$ 。 (4)

下面应用中矩形和梯形公式证明如下积分的估计，该例题取自于文献[2] 268页习题5。文献[2]只提供了该习题，没有解答，本文利用式(1)~(4)的结论给出如下详细解答。

例1 证明：$8.9<{\displaystyle {\int }_3^5\sqrt{x^2+4}\text{d}x}<9$ 。

证明：令$f\left(x\right)=\sqrt{x^2+4}$ ，记$I\left(f\right)={\displaystyle {\int }_3^5\sqrt{x^2+4}}\text{d}x$ ，由(3)式可得

$I_D\left(f\right)=\left(5-3\right)f\left(\frac{5+3}{2}\right)=2\sqrt{4^2+4}=2\sqrt{20}=4\sqrt{5}$ ，

根据误差式(4)，有

$I\left(f\right)-I_D\left(f\right)=\frac{{\left(5-3\right)}^3}{24}{f}^{″}\left(\theta \right)=\frac{1}{3}\times \frac{4}{{\left({\theta }^2+4\right)}^{\frac{3}{2}}}=\frac{4}{3{\left({\theta }^2+4\right)}^{\frac{3}{2}}},\theta \in \left(3,5\right)$ ，

由于$\theta \in \left(3,5\right)$ ，从而$\frac{4}{3{\left({\theta }^2+4\right)}^{\frac{3}{2}}}>0$ ，因此有

${\displaystyle {\int }_3^5\sqrt{x^2+4}}>4\sqrt{5}>8.9$ ，

从而证得题中左边不等式。

另外，由(1)式可得

$I_N^1\left(f\right)=\frac{5-3}{2}\left[f\left(3\right)+f\left(5\right)\right]=\sqrt{3^2+4}+\sqrt{5^2+4}=\sqrt{13}+\sqrt{29}$ ，

根据误差式(2)，有

$I\left(f\right)-I_N^1\left(f\right)=\frac{-{\left(5-3\right)}^3}{12}{f}^{″}\left(\theta \right)=\frac{-8}{3{\left({\theta }^2+4\right)}^{\frac{3}{2}}},\theta \in \left(3,5\right)$ ，

由于$\theta \in \left(3,5\right)$ ，从而$\frac{-8}{3{\left({\theta }^2+4\right)}^{\frac{3}{2}}}<0$ ，因此有

${\displaystyle {\int }_3^5\sqrt{x^2+4}}\text{d}x<\sqrt{13}+\sqrt{29}<9$ ，

从而证得题中右边不等式，即题中结论成立。

教学注记1：从例1可以看出数值积分在数学理论上也是很有用的工具，例1通过两类简单的数值积分公式的误差估计式给出了积分的估计，该方法的好处在于无需计算精确积分就给出了定积分的估计，且该估计是非常有效的。由微积分的知识，我们可得

${\displaystyle {\int }_3^5\sqrt{x^2+4}}\text{d}x={\frac{x\sqrt{x^2+4}}{2}+2\ln \left(x+\sqrt{x^2+4}\right)|}_3^5\approx \text{8}\text{.959520962808250}$ ，

尽管上式给出了精确值，然而由于涉及对数函数，要得到最后的近似值也是需要利用其它逼近手段的。

教学注记2 (表1)：Newton-Cotes公式采用等距节点作为求积节点，形式简单，易于实现。这本质上属于插值型求积公式，利用已知的等距节点$x_k\left(k=0,1,2,\cdots ,n\right)$ 结合插值多项式很容易计算出相应的求积系数$w_k\left(k=0,1,2,\cdots ,n\right)$ 。然而该公式当$n\ge 8$ 时，求积系数一般会出现负数，将导致数值不稳定性，从而在实际计算中失效[1] [3]。为此，我们模拟如下数值实验，用Newton-Cotes公式计算定积分${\displaystyle {\int }_0^1\sin x\text{d}x}$ ，计算结果如下表所示，从表中可以观察到：$n$ 从1到8，绝对误差逐步变小，但是$n$ 从9到20，绝对误差开始变大，主要原因就是Newton-Cotes公式$n$ 较大时的数值不稳定性质引起的。因此在实际计算时，往往采用低阶格式构造复合求积公式，如复合梯形公式($n=1$ )，复合Simpson公式($n=2$ )等。

Table 1. Calculation of definite integral by Newton-Cotes formula

表1. 用Newton-Cotes公式计算定积分${\displaystyle {\int }_0^1\sin x\text{d}x}$

n	精确值	数值积分	绝对误差
1	0.459697694131860	0.420735492403948	3.896220e−02
2	0.459697694131860	0.459862189870785	1.644957e−04
3	0.459697694131860	0.459770560550695	7.286642e−05
4	0.459697694131860	0.459697448597744	2.455341e−07
5	0.459697694131860	0.459697555984819	1.381470e−07
6	0.459697694131860	0.459697694435672	3.038114e−10
7	0.459697694131860	0.459697694318177	1.863164e−10
8	0.459697694131860	0.459697694130196	1.663947e−12
9	0.459697694131860	0.459697694086807	4.505318e−11
10	0.459697694131860	0.459697694038827	9.303291e−11
11	0.459697694131860	0.459697693778442	3.534181e−10
12	0.459697694131860	0.459697697130789	2.998929e−09
13	0.459697694131860	0.459697692241116	1.890745e−09
14	0.459697694131860	0.459697683183847	1.094801e−08
15	0.459697694131860	0.459698735199662	1.041068e−06
16	0.459697694131860	0.459685961017148	1.173311e−05
17	0.459697694131860	0.459642516462767	5.517767e−05
18	0.459697694131860	0.458674838166899	1.022856e−03
19	0.459697694131860	0.462647491563952	2.949797e−03
20	0.459697694131860	0.433828486370864	2.586921e−02

3. Gauss求积公式

Gauss求积公式具有如下形式[4]：

$I_G^n\left(f\right)={\displaystyle \sum _{g=0}^n{A}_{g}f\left(x_g\right)}\approx {\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ ，

其中$x_g\left(g=0,1,2,\cdots ,n\right)$ 为区间$\left[a,b\right]$ 上$n+1$ 次正交多项式的零点，$A_g={\displaystyle {\int }_a^b{l}_g\left(x\right)\text{d}x}$ ，其中$l_g\left(x\right)={\displaystyle \prod _{\begin{array}{l}p=0\\ p\ne g\end{array}}^n\frac{x-x_p}{x_g-x_p}}$ 。其相应的误差式为：

$I\left(f\right)-I_G^n\left(f\right)=\frac{f^{\left(2n+1\right)}\left(\theta \right)}{\left(2n+2\right)!}{\displaystyle {\int }_a^b{w}_{n+1}^2}\left(x\right)\text{d}x,\theta \in \left(a,b\right)$ ，

其中$w_{n+1}\left(x\right)={\displaystyle \prod _{g=0}^n\left(x-x_g\right)}$ 。

当$n=1$ ，$\left[a,b\right]=\left[-1,1\right]$ 时，Gauss求积公式为

$I_G^1\left(f\right)=f\left(\frac{-\sqrt{3}}{3}\right)+f\left(\frac{\sqrt{3}}{3}\right)$ ，

当$n=2$ ，$\left[a,b\right]=\left[-1,1\right]$ 时，Gauss求积公式为

$I_G^2\left(f\right)=\frac{5}{9}f\left(\frac{-\sqrt{15}}{5}\right)+\frac{8}{9}f\left(0\right)+\frac{5}{9}f\left(\frac{\sqrt{15}}{5}\right)$ 。

例2 利用Newton-Cotes和Gauss求积公式近似计算定积分${\displaystyle {\int }_{-1}^1{\text{e}}^{-x^2}\text{d}x}$ ，并进行比较。

解：令$f\left(x\right)={\text{e}}^{-x^2}$ ，记$I\left(f\right)={\displaystyle {\int }_{-1}^1{\text{e}}^{-x^2}}\text{d}x$ 。

用梯形公式($n=1$ )计算得

$I_N^1\left(f\right)=\frac{1+1}{2}\left({\text{e}}^{-1}+{\text{e}}^{-1}\right)\approx \text{0}\text{.735758882342885}$ ，

用Simpson公式($n=2$ )计算得

$I_N^2\left(f\right)=\frac{1+1}{6}\left({\text{e}}^{-1}+4{\text{e}}^0+{\text{e}}^{-1}\right)\approx \text{1}\text{.578586294114295}$ ，

用四点Newton-Cotes公式($n=3$ )计算得

$I_N^3\left(f\right)=\frac{1+1}{8}\left({\text{e}}^{-1}+3{\text{e}}^{\frac{-1}{9}}+3{\text{e}}^{\frac{-1}{9}}+{\text{e}}^{-1}\right)\approx \text{1}\text{.526198695807276}$ ，

用Cotes公式($n=4$ )计算得

$I_N^4\left(f\right)=\frac{1+1}{90}\left(7{\text{e}}^{-1}+32{\text{e}}^{\frac{-1}{4}}+12{\text{e}}^0+32{\text{e}}^{\frac{-1}{4}}+7{\text{e}}^{-1}\right)\approx \text{1}\text{.488745828732669}$ ，

用两点Gauss求积公式($n=1$ )计算得

$I_G^2\left(f\right)={\text{e}}^{\frac{-1}{3}}+{\text{e}}^{\frac{-1}{3}}\approx \text{1}\text{.433062621147578}$ ，

用三点Gauss求积公式($n=2$ )计算得

$I_G^3\left(f\right)=\frac{1}{9}\left(5{\text{e}}^{\frac{-3}{5}}+8{\text{e}}^0+5{\text{e}}^{\frac{-3}{5}}\right)\approx \text{1}\text{.498679595660029}$ ，

用二十点Gauss ($n=19$ )求积公式计算得

$I_G^{20}\left(f\right)\approx 1.493648265624854$ 。

把$I_G^{20}\left(f\right)$ 作为积分$I\left(f\right)$ 的参考值，经计算

$\left|I_N^1\left(f\right)-I_G^{20}\left(f\right)\right|=\text{0}\text{.757889383281969}$ ，

$\left|I_N^2\left(f\right)-I_G^{20}\left(f\right)\right|=\text{0}\text{.084938028489441}$ ，

$\left|I_G^1\left(f\right)-I_G^{20}\left(f\right)\right|=\text{0}\text{.060585644477276}$ ，

因此，两点Gauss求积公式比梯形公式(两点Newton-Cotes公式)和Simpson公式(三点Newton-Cotes公式)计算结果都要好。经计算，

$\left|I_G^2\left(f\right)-I_G^{20}\left(f\right)\right|=0.005031330035175$ ，

$\left|I_N^4\left(f\right)-I_G^{20}\left(f\right)\right|=\text{0}\text{.004902436892185}$ ，

因此三点Gauss求积公式和Cotes公式(五点Newton-Cotes公式)从误差角度来看是差不多的，也是上面计算结果中最好的。从如上结果及其下面图中显示的收敛性结果，我们可看出在求积节点个数相同情况下，对于本例，Gauss求积公式比Newton-Cotes公式更有效，而Clenshaw-Curtis的收敛速度在两者之间(图1)。

Figure 1. Convergence comparison of Newton-Cotes, Gauss quadrature and Clenshaw-Curtis formulas for calculating definite integral ${\displaystyle {\int }_0^1{\text{e}}^{-x^2}\text{d}x}$

图1. Newton-Cotes公式、Gauss求积公式与Clenshaw-Curtis公式计算定积分${\displaystyle {\int }_0^1{\text{e}}^{-x^2}\text{d}x}$ 的收敛性比较

教学注记3：在例2中，我们把$I_G^{20}\left(f\right)$ 作为积分$I\left(f\right)$ 的参考值，这是因为对于Gauss求积公式，我们有如下收敛性结果：设$f\left(x\right)\in C\left[a,b\right]$ ，则$\underset{n\to \infty }{\lim }{I}_G^n\left(f\right)=\underset{n\to \infty }{\lim }{\displaystyle \sum _{g=0}^n{A}_{g}f\left(x_g\right)}={\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ 。对于例2中的$f\left(x\right)={\text{e}}^{-x^2}\in C\left[-1,1\right]$ 满足如上假设条件，从而$n$ 越大，计算效果越好。事实上，对此例，在Matlab软件运行下，$I_G^{20}\left(f\right)\approx 1.493648265624854$ ，$I_G^{21}\left(f\right)\approx 1.493648265624854$ ，两者在计算机可显示的范围是一致的，因此我们把$I_G^{20}\left(f\right)$ 作为积分$I\left(f\right)$ 的参考值。对于任意节点Gauss求积公式的Matlab程序实现，可见参考文献[5]。

教学注记4：在例2中，我们看到相同求积节点的情况下，Gauss求积公式比Newton-Cotes的计算结果更好，其原因之一为在求积节点个数相同情形下，Gauss求积公式的代数精度要高于Newton-Cotes公式。事实上，Gauss求积公式$I_G^n\left(f\right)$ 的代数精度为$2n+1$ ，而对于Newton-Cotes公式，一般来说，当$n$ 为奇数时，$I_N^n\left(f\right)$ 的代数精度为$n$ ，当$n$ 为偶数时，$I_N^n\left(f\right)$ 的代数精度为$n+1$ 。此外，Gauss求积公式$I_G^n\left(f\right)$ 的求积节点为相应的正交多项式的$n+1$ 零点，可以证明$n+1$ 个求积节点可以生成出最多$2n+1$ 阶代数精度的求积公式–即高斯公式。事实上，可构造$2n+2$ 次多项式$p_{2n+2}\left(x\right)=w_{n+1}^2\left(x\right)={\left[{\displaystyle \prod _{g=0}^n\left(x-x_g\right)}\right]}^2$ ，该函数在$\left[-1,1\right]$ 是非负函数，且只在有限个高斯点$x_g\left(g=0,1,2,\cdots ,n\right)$ 等于零，由数学分析知识可得：${\displaystyle {\int }_{-1}^1{p}_{2n+2}\left(x\right)\text{d}x}={\displaystyle {\int }_{-1}^1{\left[{\displaystyle \prod _{g=0}^n\left(x-x_g\right)}\right]}^2\text{d}x}>0$ ，而用高斯积分计算可得$I_G^n\left(p_{2n+2}\right)=0$ ，从而说明最高代数精度的公式严格小于$2n+2$ ，即最多为$2n+1$ 。正交多项式零点即为高斯求积节点这个性质在数学理论上是非常优美的，利用此性质Gauss求积公式还可推广至带权函数积分的计算，例如${\displaystyle {\int }_{-1}^1\frac{1}{\sqrt{1-x^2}}f\left(x\right)\text{d}x}$ 和${\displaystyle {\int }_0^{+\infty }{\text{e}}^{-x}f\left(x\right)\text{d}x}$ 。

教学注记5：在教学注记2中，我们提到，Newton-Cotes公式当$n\ge 8$ 时，会出现负的求积系数从而使得格式数值不稳定。而对于Gauss求积公式，可证明，其求积系数是正的，从而是数值稳定的。事实上，由于Gauss求积公式$I_G^n\left(f\right)={\displaystyle \sum _{g=0}^n{A}_{g}f\left(x_g\right)}\approx {\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}$ 具有$2n+1$ 次代数精度，从而有，当$f\left(x\right)$ 为次数不超过$2n+1$ 的多项式时，该公式是精确的，特别地，取$2n$ 次多项式$f\left(x\right)=l_k^2\left(x\right)$ ，其中$l_k\left(x\right)={\displaystyle \prod _{\begin{array}{l}p=0\\ p\ne k\end{array}}^n\frac{x-x_p}{x_k-x_p}}$ ，则有$0<{\displaystyle {\int }_a^b{l}_k^2\left(x\right)\text{d}x}={\displaystyle \sum _{g=0}^n{A}_g{l}_k^2}\left(x_g\right)=A_k\left(k=0,1,2,\cdots ,n\right)$ 。需要注意的是，当被积函数带有奇异性的时候，直接用Gauss公式表现并不理想，此时复合梯形公式表现更稳健。为此，我们讨论用如上两种方法来计算定积分${\displaystyle {\int }_{0.0001}^1\frac{1}{x}\sin \left(\frac{1}{x}\right)\text{d}x}$ ，其计算结果放在图2，从图2中我们可以观察到，Gauss积分对于该定积分并不收敛，而复合梯形公式随着子区间细分一开始误差变小，后面也开始变大，这主要原因是被积函数$\frac{1}{x}\sin \left(\frac{1}{x}\right)$ 在$x=0.0001$ 附近会出现非常大的函数值，且是高震荡的。为进一步研究，我们把通过变量替换$t=\frac{1}{x}$ ，把原来的定积分转化成等价形式${\displaystyle {\int }_{0.0001}^1\frac{1}{x}\sin \left(\frac{1}{x}\right)\text{d}x}={\displaystyle {\int }_1^{10000}\frac{\sin t}{t}\text{d}t}$ ，我们再次用Gauss公式和复合梯形公式计算该等价定积分${\displaystyle {\int }_1^{10000}\frac{\sin t}{t}\text{d}t}$ ，结果放在图3，此时我们观察到，高斯积分仍然不收敛，而复合梯形公式是收敛的。从该例可知，没有普适性最好的数值积分公式，在实际问题中应该选取合适的数值积分方法(图2、图3)。

教学注记6：除了Newton-Cotes公式和Gauss求积公式，Clenshaw-Curtis求积公式也是常见的数值积分格式。Clenshaw-Curtis公式用Chebyshev点作为求积节点，$n+1$ 个求积节点对应的求积公式的代数精度为$n$ ，此外，若被积函数连续，Clenshaw-Curtis公式也具有收敛性。对于Gauss和Clenshaw-Curtis求积公式的详细且深刻的比较可见文献[5]。

Figure 2. Numerical performance of Gauss quadrature formula and composite trapezoidal formula for calculating singular definite integral

图2. Gauss求积公式与复合梯形公式计算奇异定积分的数值表现

Figure 3. Convergence comparison of Gauss quadrature formula and composite trapezoidal formula after variable substitution

图3. 变量替换后Gauss求积公式与复合梯形公式的收敛性比较

4. 结语

数值积分是逼近定积分的常见数值格式。本文首先讨论了Newton-Cotes公式，该公式利用等矩节点作为求积节点，格式简单，易于实现。在例1定积分的估计中，我们借助梯形公式和中矩形公式的误差估计式给出了相关证明，从此例可看出数值积分公式在数学理论上的作用。Newton-Cotes公式当$n\ge 8$ 时，通常会出现负的求积系数，从而导致数值不稳定性，因此在实际工程问题中往往结合低价格式构造复合求积公式，详细见教学注记1~2。然后，本文讨论了Gauss求积公式，该公式利用正交多项式零点作为求积节点，代数精度高，且对于连续函数是收敛的，从而是求解实际问题的有效数值格式。并且通过例2，我们比较了Newton-Cotes公式，Gauss求积公式和Clenshaw-Curtis公式的数值模拟过程，从而给出教学注记3~6。

基金项目

2023年度广东省本科高校教学质量与教学改革工程建设项目“大数据背景下的数值分析课程的教学研究和实践”(粤教高函[2024] 9号)；2024年度嘉应学院校级一流本科课程“数值分析”(嘉院教[2024] 9号)，2022年度嘉应学院教学质量与教学改革工程建设项目“一流本科课程建设背景下的《数值分析》课程的教学研究和实践”。