数值积分法

在计算数学中，不可避免地需要对函数在给定区间上求定积分。即计算

\[\int_{a}^{b}f(x)\mathrm{d}x.
\]

根据定积分的定义，定积分是积分上和与积分下和共同极限，因此只要函数可积，且区间被拆分得足够细小，那么积分就可以近似表示为

\[\int_{a}^{b} f(x)\mathrm{d}x=\lim_{n\to \infty}\sum_{i=1}^{n}\frac{b-a}{n}\cdot f(a+\frac{b-a}{n}i).
\]

因此，如果我们知道函数在足够多的\(n\)个点处的取值，就可以用积分的定义近似计算其数值积分。不过这样做的缺点是运算量大，且误差难以估计。本文要介绍的Newton-Cotes积分和Gauss积分法，是两种常用的数值积分方法。

一、Newton-Cotes公式积分

利用Newton-Cotes公式对函数进行积分，是基于之前我们说过的插值理论的。在插值理论中我们说过，如果已知函数在不同的\(n+1\)个点处的取值，就可以构造一个\(n\)次多项式作为原函数的近似。而多项式的积分是很容易计算的，这样就简化了定积分的计算。

为方便接下来的讨论，先回顾插值理论中的一些结果，设\(\omega(x)=\prod_{i=0}^{n}(x-x_i)\)。

首先，过\(n+1\)个不同的节点给出的插值多项式是唯一的，利用Lagrange插值基函数可以给出过\(x_0,x_1,\cdots,x_n\)点的插值函数为

对任何\(x\in [a,b]\)，用\(P_n(x)\)作为\(f(x)\)的近似的误差估计是

特别当\(n=1\)时和\(n=2\)时，有

简单情形：梯形积分公式与抛物线积分公式

梯形积分公式基于两个点：\((a,f(a))\)和\((b,f(b))\)，利用这两个点构造出的插值多项式最简单，是一条直线，可以很容易地写出其表达式：

对其进行积分，就得到

直观上来看，这就是以\((a,f(a))\)和\((b,f(b))\)到\(x\)轴的线段为平行边围成的梯形的面积，因而将这个求积公式叫做梯形求积公式，它是插值多项式积分的最简单情形。

如果我们需要稍复杂一些，增加一个点，最好的选择是增加\(a,b\)的中点。此时，插值多项式是一个二次多项式，即一条抛物线，可以写出其表达式为：

对其进行积分，得到的结果是

如果对比梯形求积公式，会发现，它们共同的特点是：都包含\((b-a)\)这一项，同时另一个因子都是函数值的加权平均。如果将插值多项式推广到更高的阶数，是否仍有这样的结论？下一Part的结论将给出肯定的答复。

Newton-Cotes公式

这一部分的内容，就是对梯形求积公式和抛物线求积公式的一个推广，它的要求是：在\([a,b]\)上选择的分点是等间距的，即

（下面的推导过程可以略看，与计算过程无关）根据前言部分的结论，我们有

于是

令

则有

下面要计算\(A_i\)，就要用到节点等间距的要求，用\(a+ih\)替代\(x_i\)，并令\(x=a+th\)，这样\(t\)就在\([0,n]\)中连续取值，并且

代入就得到

（略看部分结束）观察\(A_i\)式，会发现\((nh)\)就是区间长度\(b-a\)，再令剩余部分为\(c_{i}^{(n)}\)，就得到积分式为：

当然，\(c_i^{(n)}\)的计算有一些繁琐，但是注意到\(c_i^{(n)}\)的表达式中只有\(n,i\)而没有\(h\)，这就说明\(c_i^{(n)}\)的计算不依赖于区间，只依赖于节点个数与节点的排位，因此可以事先计算出来以供查阅。这样\(c_i^{(n)}\)被称为Newton-Cotes系数，有表如下：

可以发现，当\(n\le 7\)时所有的\(c_{i}^{(n)}\)都是正数，而当\(n\ge 8\)时Cotes系数出现了负数，这是我们不使用过多结点进行Newton-Cotes公式计算数值积分的原因，因为出现了负数会导致误差估计难以确定，这是我们下一Part将讨论的内容。

现在我们用Cotes公式计算\(\sin x^2\)在\([-1,1]\)上的积分，不妨取\(n=2\)和\(n=3\)，插值多项式的图示如下（黑色为原函数，红色为插值多项式）：

当\(n=2\)时，计算得

\[\int_{-1}^{1}\sin(x^2)\mathrm{d}x\approx 2\left[\frac{1}{6}\sin (-1)^2+\frac{4}{6}\sin(0)^2+\frac{1}{6}\sin(1^2) \right]=0.5609807,
\]

当\(n=3\)时，计算得

\[\int_{-1}^{1}\sin(x^2)\mathrm{d}x\approx 2\left[\frac{1}{8}\sin(-1)^2+\frac{3}{8}\sin(-\frac{1}{3})^2+\frac{3}{8}\sin(\frac{1}{3})^2+\frac{1}{8}\sin 1^2 \right]=0.5870594.
\]

而积分的更精确值是\(0.6205366\)，从图上也可以看出红色线要比黑色线普遍低一些，这代表对于此函数，Newton-Cotes公式积分的结果会比精确值低。此时，如果选择\(n=1\)，结果会怎样？

误差分析

首先，在使用Newton-Cotes公式计算积分时，常常不会选择太多的节点，因为\(n\ge 8\)时，Cotes系数可能会出现负数。为什么对于负数Cotes系数时一般不选择，这是考虑到计算\(f(x_i)\)时可能出现的偏误。一般说来，即便\(f\)已知，计算\(f(x_i)\)也会有一定的舍入误差、截断误差等，这一般会有一定的绝对误差限\(\varepsilon\)，不妨假设计算出来的有误差的值为\(\tilde f(x_i)\)。如果Cotes系数全是正数，那么结合\(\sum_i c_i^{(n)}=1\)，有

这使得误差处于可以控制的范围内，这个性质称为数值稳定性。如果Cotes系数存在负数，这个纯由计算误差导致的误差就可能变得难以控制。

在讨论误差估计之前，先回顾一下积分中值定理。

积分中值定理：若函数\(f(x)\)、\(g(x)\)在闭区间\([a,b]\)上连续，且\(g(x)\)在\([a,b]\)上不变号，则存在一点\(\varepsilon\in[a,b]\)，使得

\[\int_{a}^{b}f(x)g(x)\mathrm{d}x=f(\varepsilon)\int_{a}^{b}g(x)\mathrm{d}x.
\]

现在，我们讨论Newton-Cotes公式的误差估计，方便起见，我们只讨论比较简单的梯形公式与抛物线公式。首先对于梯形求积公式，有如下的定理：若\(f(x)\in C^{2}[a,b]\)，则梯形求积公式有误差估计

显然

\[\begin{aligned}
R(f)&=\int_{a}^{b}[f(x)-P_1(x)]\mathrm{d}x\\
&=\int_{a}^{b} R_1(x)\mathrm{d}x\\
&=\int_{a}^{b}\frac{f^{(2)}(\xi)}{2}(x-a)(x-b)\mathrm{d}x.
\end{aligned}
\]

要注意，这里\(\xi\)是\(x\)的函数\(\xi(x)\)，所以不能贸然将其提到积分限以外，而应运用积分中值定理，因\((x-a)(x-b)\)在\([a,b]\)上非正，有

\[\int_{a}^{b}\frac{f^{(2)}(\xi(x))}{2}(x-a)(x-b)\mathrm{d}x=\frac{f^{(2)}(\eta)}{2}\int_{a}^{b}(x-a)(x-b)\mathrm{d}x.
\]

对于下面的积分，我们常用以下变换：\(x=a+th\)，\(h=\frac{b-a}{n}\)。这里\(n=1\)，所以

\[\int_{a}^{b}(x-a)(x-b)\mathrm{d}x=h^3\int_{0}^{1}t(t-1)\mathrm{d}t=-\frac{(b-a)^3}{6},
\]

所以

\[R(f)=-\frac{(b-a)^3f^{(2)}(\eta)}{12}.
\]

对于抛物线求积公式，有与梯形求积公式类似的、但又有所不同的误差估计：若\(f(x)\in C^4[a,b]\)，则

这是因为此时\(n\)为偶数。这里的证明，需要在下一部分给出代数精度这一概念时才能给出。

复化求积与逐次分半

在插值理论中，我们知道存在某些函数，即使节点无限加密，也难以对原始函数很好的拟合（尤其是负次数函数），即Runge现象。为了克服Runge现象，往往使用分段插值以缩小误差。同理，如果拟合效果差，用多项式代替求积分的效果也自然不好，我们可能也会对分段插值得到的分段函数进行分段求积，这便是复化求积公式。

这里不对复化梯形公式与复化抛物线求积公式过多介绍，直接给出公式。对于梯形求积公式，把\([a,b]\)区间\(n\)等分后得到的节点\(x_k=a+kh\)，有

对于抛物线求积公式，把\([a,b]\)区间\(2n\)等分后得到的节点为\(x_k=a+kh\)，这里\(h=(b-a)/2n\)，于是

其误差估计为：

对于梯形求积公式，有

\[\begin{aligned}
R(f,T_n)&=\sum_{k=0}^{n-1}\left(-\frac{h^3}{12}f^{(2)}(\eta_k) \right),\quad \eta_k\in(x_{k},x_{k+1})\\
&=-\frac{h^3}{12}\sum_{k=0}^{n-1}f^{(2)}(\eta_k)\\
&=-\frac{h^3}{12}nf^{(2)}(\eta),\quad \eta\in(a,b)\\
&=-\frac{b-a}{12}h^2f^{(2)}(\eta).
\end{aligned}
\]

这里第三个等号是因为\(f^{(2)}(x)\)作为连续函数的性质。

对于抛物线求积公式，有

\[\begin{aligned}
R(f,S_n)&=\sum_{k=0}^{n-1}\left(-\frac{(2h)^5}{2880}f^{(4)}(\eta_k) \right),\quad \eta_k\in(x_k,x_{k+1})\\
&=-\frac{(2h)^5}{2880}\sum_{k=0}^{n-1}f^{(4)}(\eta_k)\\
&=-\frac{2h^5}{180}nf^{(4)}(\eta),\quad \eta\in(a,b)\\
&=-\frac{b-a}{180}h^{4}f^{(4)}(\eta).
\end{aligned}
\]

接下来介绍逐次分半法，它是一种计算复化求积的方式。虽然有了梯形求积公式的误差估计，但是这是一种事前估计，即为了达到我们预期的误差限，往往会选择超过实际所需非常多的节点，这是一种计算量上的浪费。而逐次分半法，每次对节点加密一倍，从而获得更精确的积分结果，这样将得到一个积分数列\(T_1,T_2,T_4,\cdots\)，当相邻两个积分结果小于一定值时，就不继续计算，将最终结果作为返回值。

逐次分半法还有一个显著优点是，重复的节点不需要重复计算。不妨设将\([a,b]\)四等分，节点是\(a,x_1,x_2,x_3,b\)，看下面的计算过程：

在这样的计算过程下，重复的结点没有加入计算，所以每次只需要计算新加入的节点即可。

一般将停步的定值（即最终相邻两个积分值的差）设定为\(3\varepsilon\)，\(\varepsilon\)是误差限。

二、代数精度

在Newton-Cotes公式中，最终将求积公式的形式表现为

这也是我们要介绍的求积公式的一般形式，注意各个\(A_i\)是不依赖于被积函数\(f\)的。这里介绍代数精度概念，一方面是作为评判的标准，另一方面是为下面Gauss求积公式做准备。

代数精度的定义

前面所叙述的求积，都是以某种方式将被积函数\(f(x)\)用\(n\)次多项式\(P_n(x)\)代替，对多项式求积。特别在Newton-Cotes公式求积分时，使用的是插值多项式，而我们知道如果\(f(x)\)是不高于\(n\)次的多项式，则\(P_n(x)\equiv f(x)\)，这时候积分的结果是精确的；但如果\(f(x)\)是高于\(n\)次的多项式，那么\(f^{(n+1)}(x)\)不会恒等于\(0\)，用\(P_n(x)\)替代\(f(x)\)求积分就会存在误差。代数精度这个概念，就是用来表征某种积分法是相当于用几次多项式来替代求积的。

现给出代数精度的定义。对一个一般求积公式

其中\(A_k\)是不依赖于函数\(f(x)\)的常数。如果求积公式中的\(f(x)\)是任意一个次数不高于\(m\)次多项式时，求积公式的约等号可改为精确成立的等号，但存在\(m+1\)次多项式\(f(x)\)使求积公式不能精确成立，就称求积公式具有\(m\)次代数精度。

Newton-Cotes公式的代数精度

在引入Gauss求积公式之前，先对Newton-Cotes公式的代数精度进行讨论。

由于Newton-Cotes公式使用\(n\)次插值多项式进行替代求值，所以对于有\(n+1\)个节点\(x_0,x_1,\cdots,x_n\)的Cotes公式，Newton-Cotes求积至少应当具有\(n\)次代数精度，因为当\(f(x)\in\mathcal P_n(\mathbb{R})\)时\(f(x)\equiv P_n(x)\)。比如，梯形求积公式对于一次函数的积分是精确的，而抛物线求积公式对于二次函数的积分也是精确的。

但值得注意的是，Newton-Cotes公式对于\(n\)为偶数的情形，可以具有\(n+1\)次代数精度。

记任意\(n+1\)次多项式为

\[q_{n+1}(x)=\sum_{k=0}^{n+1}b_{k}x^{k},\quad q_{n+1}^{(n+1)}(x)=(n+1)!b_{n+1},
\]

用\(P_n(x)\)代表其\(n\)次插值多项式，于是其误差为

\[R_n(x)=\frac{q_{n+1}^{(n+1)}(\xi)}{(n+1)!}\omega_n(x)=b_{n+1}\omega_n(x),
\]

在前面，我们已经得到了这样的结果：

\[\int_{a}^{b}\omega_n(x)\mathrm{d}x=h^{n+2}\int_{0}^{n}t(t-1)\cdots(t-n)\mathrm{d}t,
\]

接下来的目标是证明上面的积分为\(0\)，这样就代表用\(P_n(x)\)替代\(q_{n+1}(x)\)的积分是精确的。令\(n=2k\)并引进变换\(u=t-k\)，有

\[\begin{aligned}
&\quad \int_{0}^{n}t(t-1)\cdots(t-n)\mathrm{d}t\\
&=\int_{0}^{2k}t(t-1)\cdots(t-k)(t-k-1)\cdots(t-2k+1)(t-2k)\mathrm{d}t\\
&=\int_{k}^{k}(u+k)(u+k-1)\cdots u(u-1)\cdots(u-k-1)(u-k)\mathrm{d}u,
\end{aligned}
\]

令积分号里面的部分为\(H(u)\)，即\(H(u)=(u+k)\cdots u(u-1)\cdots(u-k)\)，有

\[\begin{aligned}
H(-u)&=(-u+k)(-u+k+1)\cdots(-u)(-u-1)\cdots(-u-k) \\
&=(-1)^{2k+1}H(u)\\
&=-H(u),
\end{aligned}
\]

从而\(H(u)\)是一个奇函数，这说明积分结果是\(0\)，所以

\[\int_{a}^{b}\omega_n(x)\mathrm{d}x=0,
\]

即\(P_n(x)\)代替\(q_{n+1}(x)\)使用Newton-Cotes求积公式是精确的。

证明抛物线公式的误差估计

之前，我们给出了梯形求积公式与抛物线求积公式的误差估计，并证明了梯形求积公式的误差估计，但没有给出抛物线公式的误差估计。这是因为抛物线求积公式具有三次代数精度，这使得它可以用更高次的插值多项式逼近。

再回顾抛物线公式的误差估计：

由于可以构造三次多项式，故除了\(f(a),f(\frac{a+b}{2}),f(b)\)以外，还可以增加一个条件，这里选择的是\(f\'(\frac{a+b}{2})\)。利用Newton插值法，可知\(f\'(\frac{a+b}{2})\)实际上是在\((\frac{a+b}{2})\)处选择了重节点。

构造一个三次插值多项式\(P_3(x)\)满足

\[P_3(a)=f(a),\quad P_3(b)=f(b),\\
P_3\left(\frac{a+b}{2} \right)=f\left(\frac{a+b}{2} \right),\\
P_3\’\left(\frac{a+b}{2} \right)=f\’\left(\frac{a+b}{2} \right).
\]

从而插值函数的误差是

\[f(x)-P_3(x)=\frac{f^{(4)}(\xi(x))}{4!}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b),
\]

易知后面的部分在\([a,b]\)上不变号。积分，由积分中值定理，并引入变换\(x=a+t(b-a)\)，就有

\[\begin{aligned}
&\quad \int_{a}^{b}\frac{f^{(4)}(\xi(x))}{4!}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b)\mathrm{d}x\\
&=\frac{f^{(4)}(\eta)}{4!}\int_{a}^{b}(x-a)\left(x-\frac{a+b}{2} \right)^2(x-b)\mathrm{d}x\\
&=\frac{(b-a)^5f^{(4)}(\eta)}{4!}\int_0^1 t\left(t-\frac{1}{2}\right)^2 (t-1)\mathrm{d}t\\
&=-\frac{(b-a)^5f^{(4)}(\eta)}{4!\cdot 120}.
\end{aligned}
\]

容易计算最后一个积分值为\(-1/120\)，且\(4!\cdot 120=2880\)，故结论得证。

三、Gauss求积公式

Gauss型求积公式是对Newton-Cotes求积公式的改进。Newton-Cotes求积公式基于插值多项式推出，选择了最方便的等距节点，并且根据\(A_i\)中与区间无关的部分编制了Cotes系数表，从而使计算变得方便。

Gauss型求积公式则是为了更精确的积分结果，选择更复杂的节点与对应的系数\(A_i\)。如何获得更精确的积分结果，就是我们接下来需要考虑的问题。

Gauss求积的目标

先明确Gauss型求积公式的目标：在节点数\(n\)给定时，最大化代数精度。要最大化代数精度，首先就要知道代数精度最高能达到多少，这就成为一个规划问题：

写出\(m\)次多项式的一般形式：

于是

这里\(\int_{a}^{b}x_i\mathrm{d}x\)是定值，故使用\(\mu_i\)替代，从而我们得到一个等式：

注意这个等式中的几个量：

• \(a_i\)是任意实数。因为一组\(a_i\)唯一决定一个\(q_m(x)\)，所以\(a_i\)任意取值等式都应该成立。
• \(\mu_i\)是定值。因为\(\mu_i=\int_{a}^{b}x^{i}\mathrm{d}x\)。
• \(A_k\)和\(x_k\)是待取的定值，即取一组\(A_k,x_k\)使得对任何\(a_i\)都有等式成立。

由于\(a_i\)是任意实数，所以可以将等式展开成方程组：\(\mu_i=\sum_{k=1}^{n}A_kx_{k}^{i},i=0,1,\cdots,m\)，即

因未知数个数为\(2n\)，为使这个方程组成立，必有\(m\le 2n-1\)，且如果存在一组\(A_k,x_k\)使得\(m=2n-1\)，则\(2n-1\)就是\(n\)个互异节点能达到的最大代数精度。事实上，这样的\(A_k,x_k\)必定是可以取到的。

接下来讨论如何选择节点，即解上面的非线性方程组。为方便起见，我们将定积分的上下限控制为\([-1,1]\)，这样有

因为我们假定\(q(x)\)是\(2n-1\)次多项式，所以用一个含\(x_1,\cdots,x_n\)的\(n\)次多项式，这里就取为\(\omega(x)=\prod_{i=1}^{n}(x-x_i)\)去除\(q(x)\)，就有

这里多项式除法的存在唯一性保证了\(p(x),r(x)\)都是唯一的\(n-1\)次多项式，且

现在，最关键的一步在于，我们希望对任何不超过\(n-1\)次的多项式\(p(x)\)，都有

这等价于以下\(n\)个算式：

一般从中解出\(x_1,x_2,\cdots,x_n\)的值，可以达到降维的效果，求解难度也大幅下降；并且，从这里可以看出，\(x_1,\cdots,x_n\)的值是不依赖于函数\(f\)的。这时候

再利用求积公式对\(f(x)=x^{i}\)，\(i=0,1,\cdots,n-1\)精确成立，代入已经求得的\(x_1,\cdots,x_n\)的值，就可以解出\(A_1,\cdots,A_n\)。事实上，可以解得

这相当于把一个\(2n\)个未知数的非线性方程，转化成一个含\(n\)个未知数的非线性方程和一个\(n\)个未知数的线性方程组求解。

Gauss-Legendre求积公式

在之前的讨论中，我们使用\(\omega(x)\)，利用\(x^{i}\omega(x)\)的积分在\([-1,1]\)上为\(0\)来求解\(x_1,\cdots,x_n\)。而Gauss-Legendre求积公式，则将\(\omega(x)\)与Legendre多项式建立联系来作为除函数，这样做的好处是，对任何\(i=0,1,\cdots,n-1\)，因Legendre多项式的性质都有

由Legendre多项式的性质，如果取它的\(n\)个零点作为\(x_1,\cdots,x_n\)，则

并且

这些系数被制成表，需要使用时可以直接查阅。

在实际应用Gauss-Legendre公式求积时，一般步骤是先将其作积分限变换转化为\([-1,1]\)上的积分，再根据公式代入求解。

现使用Gauss-Legendre公式计算\(\sin(x^2)\)在\([-1,1]\)上的积分，之前已经计算过其积分的准确值大约是\(0.6205366\)。当\(n=4\)时，有

\[\int_{-1}^{1}\sin (x^2)\mathrm{d}x\approx \sum_{k=1}^{4}A_k\sin(x_k^2)=0.6203309.
\]

相比Newton-Cotes公式得到的结果，Gauss-Legendre积分显然更为精确。

尽管Gauss-Legendre求积有着更高的精度，但是随着\(n\)的增大，节点\(x_k\)的变化并没有什么规律，所以\(f(x_k)\)的值不能重复利用。