【科学计算】数值积分

数值积分法

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

\[\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公式
    • 误差分析
    • 复化求积与逐次分半
  • 二、代数精度
    • 代数精度的定义
    • Newton-Cotes公式的代数精度
    • 证明抛物线公式的误差估计
  • 三、Gauss求积公式
    • Gauss求积的目标
    • Gauss-Legendre求积公式

一、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\)点的插值函数为

\[P_n(x)=\sum_{i=0}^{n}\left(\prod_{j\ne i}\frac{x-x_i}{x_j-x_i} \right)f(x_i)=\sum_{i=0}^{\infty}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) \]

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

\[R(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega(x). \]

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

\[R_1(x)=\frac{f^{(2)}(\xi)}{2}(x-x_0)(x-x_1),\\ R_2(x)=\frac{f^{(3)}(\xi)}{6}(x-x_0)(x-x_1)(x-x_2). \]

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

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

\[P_1(x)=f(a)+\frac{f(b)-f(a)}{b-a}(x-a), \]

对其进行积分，就得到

\[\int_{a}^{b}P_1(x)\mathrm{d}x=(b-a)f(a)+\frac{f(b)-f(a)}{b-a}\frac{(b-a)^2}{2}=(b-a)\left(\frac{1}{2}f(a)+\frac{1}{2}f(b) \right). \]

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

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

\[\begin{aligned} &\quad P_2(x) \\ &=\frac{(x-a)(x-\frac{a+b}{2})}{(b-a)(b-\frac{a+b}{2})}f(b)+\frac{(x-a)(x-b)}{(\frac{a+b}{2}-a)(\frac{a+b}{2}-b)}f\left(\frac{a+b}{2} \right)+\frac{(x-\frac{a+b}{2})(x-b)}{(a-\frac{a+b}{2})(a-b)}f(a)\\ &=\frac{2}{(b-a)^2}\bigg[(x-a)\left(x-\frac{a+b}{2} \right)f(b)+\\ & \qquad 2(x-a)(x-b)f\left(\frac{a+b}{2} \right)+\left(x-\frac{a+b}{2}(x-b)\right)f(a) \bigg] \end{aligned} \]

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

\[\int_{a}^{b}P_2(x)\mathrm{d}x=(b-a)\left[\frac{1}{6}f(a)+\frac{4}{6}f\left(\frac{a+b}{2} \right)+\frac{1}{6}f(b) \right]. \]

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

Newton-Cotes公式

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

\[x_i=a+ih,\quad i=0,1,2,\cdots,n,h=\frac{b-a}{n}. \]

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

\[P_n(x)=\sum_{i=0}^{n}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i), \]

于是

\[\int_{a}^{b}P_n(x)\mathrm{d}x=\int_{a}^{b}\left(\sum_{i=0}^{n}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}f(x_i) \right)\mathrm{d}x=\sum_{i=0}^{n}\left(\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)} \right)f(x_i). \]

令

\[A_i=\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}, \]

则有

\[\int_{a}^{b}P_n(x)\mathrm{d}x=\sum_{i=0}^{n}A_if(x_i). \]

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

\[\omega(x)=\prod_{i=0}^{n}(x-x_i)=\prod_{i=0}^{n}[(a+ih)-(a+th)]=h^{n+1}\cdot t(t-1)\cdots(t-n),\\ \omega'(x_i)=\prod_{j\ne i}^{n}(x_i-x_j)=h^n (-1)^{n-i}[i!(n-1)!]. \]

代入就得到

\[\begin{aligned} A_i&=\int_{a}^{b}\frac{\omega(x)}{(x-x_i)\omega'(x_i)}\mathrm{d}x\\ &\xlongequal{x=a+th}\int_{0}^{n}\frac{h^{n+1}t(t-1)\cdots(t-n)}{(-1)^{n-i}h^{n}[i!(n-1)!]h(t-i)}h\mathrm{d}t \\ &=\frac{(-1)^{n-i}h}{i!(n-1)!}\int_{0}^{n}\frac{t(t-1)\cdots(t-n)}{t-i}\mathrm{d}t\\ &=(nh)\frac{(-1)^{n-i}}{n\cdot i!(n-i)!}\int_{0}^{n}\frac{t(t-1)\cdots(t-n)}{t-i}\mathrm{d}t. \end{aligned} \]

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

\[\int_{a}^{b}P_n(x)\mathrm{d}x=(b-a)\sum_{i=0}^{n}c_i^{(n)}f(x_i). \]

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

\[\begin{array}{c|cc} \hline & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline 1 & \frac{1}{2} & \frac{1}{2} \\ \hline 2 & \frac{1}{6} & \frac{4}{6} & \frac{1}{6} \\ \hline 3 & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} & \frac{1}{8} \\ \hline 4 & \frac{7}{90} & \frac{32}{90} & \frac{12}{90} & \frac{32}{90} & \frac{7}{90} \\ \hline 5 & \frac{19}{288} & \frac{75}{288} & \frac{50}{288} & \frac{50}{288} & \frac{75}{288} & \frac{19}{288} \\ \hline 6 & \frac{41}{840} & \frac{216}{840} & \frac{27}{840} & \frac{572}{840} & \frac{27}{840} & \frac{216}{840} & \frac{41}{840} \\ \hline 7 & \frac{751}{17280} & \frac{3577}{17280} & \frac{1323}{17280} & \frac{2989}{17280} & \frac{2989}{17280} & \frac{1323}{17280} & \frac{3577}{17280} & \frac{751}{17280} \\ \hline 8 & \frac{989}{28350} & \frac{5888}{28350} & \frac{-928}{28350} & \frac{10496}{28350} & \frac{-4540}{28350} & \frac{10496}{28350} & \frac{-925}{28350} & \frac{5888}{28350} & \frac{989}{28350}\\ \hline \end{array} \]

可以发现，当\(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\)，有

\[\begin{aligned} &\quad \left|(b-a)\sum_{i=1}^{n}c_i^{(n)}f(x_i)-(b-a)\sum_{i=1}^{n}c_i^{(n)}\tilde f(x_i) \right| \\ &=(b-a)\left| \sum_{i=1}^{n}c_i^{(n)}[f(x_i)-\tilde f(x_i)] \right|\\ &\le (b-a)\sum_{i=1}^{n}|c_i^{(n)}||f(x_i)-\tilde f(x_i)|\\ &\le (b-a)\varepsilon \sum_{i=1}^{n}|c_i^{(n)}|\\ &=(b-a)\varepsilon. \end{aligned} \]

这使得误差处于可以控制的范围内，这个性质称为数值稳定性。如果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]\)，则梯形求积公式有误差估计

\[R(f)=\int_{a}^{b}f(x)\mathrm{d}x-\frac{b-a}{2}(f(a)+f(b))=-\frac{(b-a)^3}{12}f^{(2)}(\eta),\quad \eta\in(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]\)，则

\[R(f)=-\frac{(b-a)^5}{2880}f^{(4)}(\eta). \]

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

复化求积与逐次分半

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

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

\[\begin{aligned} I&\approx \sum_{k=0}^{n-1}\frac{x_{k+1}-x_k}{2}(f(x_{k+1})+f(x_{k}))\\ & =\frac{h}{2}\left[f(a)+f(b)+2\sum_{k=1}^{n-1}f(a+kh)\right]\\ &\xlongequal{def}T_n. \end{aligned} \]

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

\[\begin{aligned} I&\approx \sum_{k=1}^{n}\frac{h}{3}[f(x_{2k-2})+4f(x_{2k-1})+f(x_{2k})] \\ &=\frac{h}{3}\left[ f(a)+f(b)+4\sum_{k=1}^{n}f(x_{2k-1})+2\sum_{k=1}^{n-1}f(x_{2k}) \right]\\ &\xlongequal{def} S_n. \end{aligned} \]

其误差估计为：

\[R(f,T_n)=-\frac{b-a}{12}h^2f^{(2)}(\eta),\quad \eta\in(a,b),\\ R(f,S_n)=-\frac{b-a}{180}h^4f^{(4)}(\eta),\quad \eta\in(a,b). \]

对于梯形求积公式，有

\[\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\)，看下面的计算过程：

\[T_1=(b-a)\left(\frac{f(a)}{2}+\frac{f(b)}{2} \right),\\ T_2=\frac{(b-a)}{2}\left(\frac{f(a)}{2}+\frac{f(b)}{2}+f(x_2) \right)=\frac{T_1}{2}+\frac{b-a}{2}f(x_2),\\ T_4=\frac{b-a}{4}\left(\frac{f(a)}{2}+\frac{f(b)}{2}+f(x_1)+f(x_2)+f(x_3) \right)=\frac{T_2}{2}+\frac{b-a}{4}(f(x_1)+f(x_3)),\\ \cdots \]

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

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

二、代数精度

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

\[\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{i=0}^{n}A_if(x_i).\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{i=0}^{n}A_if(x_i). \]

这也是我们要介绍的求积公式的一般形式，注意各个\(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)\)求积分就会存在误差。代数精度这个概念，就是用来表征某种积分法是相当于用几次多项式来替代求积的。

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

\[\int_{a}^{b}f(x)\mathrm{d}x\approx \sum_{k=0}^{n}A_kf(x_k), \]

其中\(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求积公式是精确的。

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

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

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

\[R(f)=-\frac{(b-a)^5}{2880}f^{(4)}(\eta),\quad \eta\in(a,b). \]

由于可以构造三次多项式，故除了\(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\)给定时，最大化代数精度。要最大化代数精度，首先就要知道代数精度最高能达到多少，这就成为一个规划问题：

\[\max m,\\ \mathrm{s.t.}\forall q\in\mathcal P_{m}(\mathbb{R}),\quad \sum_{k=1}^{n}A_kq(x_k)=\int_{a}^{b}q(x)\mathrm{d}x. \]

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

\[q_m(x)=\sum_{i=0}^{m}a_ix^{i}, \]

于是

\[\int_{a}^{b}q_m(x)\mathrm{d}x=\sum_{i=0}^{m}\int_{a}^{b}a_ix^{i}\mathrm{d}x=\sum_{i=0}^{m}a_i\int_{a}^{b}x^{i}\mathrm{d}x\xlongequal{def}\sum_{i=0}^{m}a_i\mu_i, \]

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

\[\sum_{i=0}^{m}a_i\mu_i=\sum_{k=1}^{n}A_kq(x_k)=\sum_{k=1}^{n}\sum_{i=0}^{m}A_ka_ix_k^{i}=\sum_{i=0}^{m}a_i\sum_{k=1}^{n}A_kx_k^{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\)，即

\[\begin{array}{rl} A_1+A_2+\cdots+A_n&=\mu_0,\\ A_1x_1+A_2x_2+\cdots+A_nx_n&= \mu_1,\\ A_1x_1^2+A_2x_2^2+\cdots+A_nx_n^2&=\mu_2,\\ \cdots & \cdots \\ A_1x_1^{m}+A_1x_2^{m}+\cdots+A_nx_n^{m}&=\mu_m. \end{array} \]

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

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

\[\int_{a}^{b}q(x)\mathrm{d}x\xlongequal{x=\frac{b+a}{2}+\frac{b-a}{2}t}\frac{b-a}{2}\int_{-1}^{1}q\left(\frac{b+a}{2}+\frac{b-a}{2}t \right)\mathrm{d}t. \]

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

\[q(x)=p(x)\omega(x)+r(x), \]

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

\[\int_{-1}^{1}q(x)\mathrm{d}x=\int_{-1}^{1}p(x)\omega(x)\mathrm{d}x+\int_{-1}^{1}r(x)\mathrm{d}x. \]

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

\[\int_{-1}^{1}p(x)\omega(x)\mathrm{d}x=0,\quad (1), \]

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

\[\int_{-1}^{1}(x-x_1)\cdots(x-x_n)=0,\\ \int_{-1}^{1}x(x-x_1)\cdots(x-x_n)=0,\\ \cdots \\ \int_{-1}^{1}x^{n-1}(x-x_1)\cdots(x-x_n)=0. \]

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

\[\int_{-1}^{1}p(x)\mathrm{d}x=\int_{-1}^{1}r(x)\mathrm{d}x, \]

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

\[A_k=\int_{0}^{1}\frac{\omega(x)}{(x-x_k)\omega'(x_k)}\mathrm{d}x. \]

这相当于把一个\(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多项式的性质都有

\[\int_{-1}^{1}x^{i}\frac{\mathrm{d}^n}{\mathrm{d}x^n}(x^2-1)^{n}\mathrm{d}x\xlongequal{def}\int_{-1}^{1}x^{i}P(x)\mathrm{d}x=0, \]

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

\[\frac{\mathrm{d}^n}{\mathrm{d}x^{n}}(x^2-1)^{n}=(x-x_1)(x-x_2)\cdots(x-x_n). \]

并且

\[A_k=\frac{2(1-x_{k}^2)}{[nP_{n-1}(x_k)]^2},k=1,2,\cdots,n. \]

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

\[\begin{array}{|c|c|c:c|c|c|} \hline n & x_k & A_k & n & x_k & A_k \\ \hline 2 & \pm 0.5773503 & 1.0000000 & 8 & \pm 0.9602899 & 0.1012285 \\ 4 & \pm 0.8611363 & 0.3478548 & & \pm 0.7966665 & 0.2223810 \\ & \pm 0.3399810 & 0.6521452 & & \pm 0.5255324 & 0.3137066 \\ 6 & \pm 0.9324695 & 0.1713245 & & \pm 0.1834346 & 0.3626838 \\ & \pm 0.6612094 & 0.3607616 \\ & \pm 0.2386192 & 0.4679139 \\ \hline \end{array} \]

在实际应用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)\)的值不能重复利用。