数值积分

数值积分的定义

张瑞
中国科学技术大学数学科学学院

数值积分

定义 1.
若$\omega_i$, $i=0,1,\cdots,n$与函数$f(x)$无关，则称

\[I_n(f)=\sum_{i=0}^n \omega_i f(x_i) \]

为数值积分公式，$\omega_i$称为积分系数。

例 1. 如下公式，均为数值积分公式

  • $\displaystyle\int_a^b f(x)dx\approx (b-a)f(a)$
  • $\displaystyle\int_a^b f(x)dx\approx (b-a)f(\frac{a+b}2)$
  • $\displaystyle\int_a^b f(x)dx\approx \frac{b-a}2f(a)$

问题. 这么多公式，如何判定公式的优劣？

定义 2.
对于数值积分公式$I_n(f)$，若有

\[I_n(x^j)=\int_a^b x^jdx, j=0,1,2,\cdots,k \]

而$I_n(x^{k+1})\neq\int_a^b x^{k+1}dx$，则称数值积分公式具有$k$阶代数精度。

  • 显然，具有$k$阶代数精度的数值积分公式，对于任意次数不高于$k$次的多项式，积分计算没有误差。
  • 如果一个数值积分公式，对于任意次数不高于$k$次的多项式，误差为0，则该数值积分公式具有至少$k$阶代数精度
  • 如果一个数值积分公式，对于一个$k$次多项式，误差不为0,则该数值积分公式至多为$k-1$阶代数精度

例 2. 计算如下数值积分公式的代数精度

  1. $\int_a^b f(x)dx\approx (b-a)f(a)$
  2. $\int_a^b f(x)dx\approx (b-a)f(\frac{a+b}2)$
  3. $\int_a^b f(x)dx\approx \frac{b-a}2f(a)$

解. 按照代数精度的定义，逐步验证$x^i$的积分。

$x^j$ 积分值 公式1 公式2 公式3
$x^0$ $(b-a)$ $(b-a)$ $(b-a)$ $\frac{b-a}2$
$x^1$ $\frac{b^2-a^2}2$ $(b-a)a$ $(b-a)\frac{b+a}2$
$x^2$ $\frac{b^3-a^3}3$ $(b-a)(\frac{b+a}2)^2$

注. 代数精度可以给出评判数值积分公式优劣的一个标准。

例 3. 确定系数$\omega_i$，使得求积公式

\[\int_a^bf(x)dx\approx \omega_0 f(a)+\omega_1 f(\frac{a+b}2)+\omega_2f(b) \]

有尽可能高的代数精度，并求这个代数精度。

解. 令$f(x)$为$1$, $x$, $x^2$，代入数值积分公式使之精确成立，则可得方程组

\[\begin{cases} \omega_0+\omega_1+\omega_2&=b-a \\ a \omega_0 +\frac{a+b}2 \omega_1+b \omega_2&=\frac{b^2-a^2}2 \\ a^2 \omega_0+(\frac{a+b}2)^2 \omega_1+b^2 \omega_2&=\frac{b^3-a^3}3 \\ \end{cases} \]

得到了一个关于$\omega_i$的线性方程组。

\[\begin{pmatrix} 1 & 1 & 1 \\ a & \frac{a+b}2 & b \\ a^2 & (\frac{a+b}2)^2 & b^2 \end{pmatrix} \begin{pmatrix} \omega_0 \\ \omega_1 \\ \omega_2 \end{pmatrix} =\begin{pmatrix} b-a \\ \frac{b^2-a^2}2 \\ \frac{b^3-a^3}3 \end{pmatrix} \]

可解得

\[\omega_0=\frac{b-a}6, \omega_1=\frac{4(b-a)}6 , \omega_2=\frac{b-a}6 \]

注. 行(1)$\times(-\frac{b+a}2)+$行(2)，可以得到

\[-\omega_0+\omega_2=0 \]

行(1)$\times(-(\frac{b+a}2)^2)+$行(3)，可以得到

\[\omega_0+\omega_2=\frac{b-a}3 \]

或者，取$f(x)$为$1$, $x-\frac{a+b}2$, $(x-\frac{a+b}2)^2$，则可以得到线性方程组

\[\begin{cases} \omega_0+\omega_1+\omega_2&=b-a \\ -\frac{b-a}2 \omega_0 + \frac{b-a}2 \omega_2&=0 \\ (\frac{b-a}2)^2 \omega_0+ (\frac{b-a}2)^2 \omega_2&=\frac{(b-a)^3}{12} \\ \end{cases} \]

容易验证，对于$(x-\frac{a+b}2)^3$，数值积分为

\[-(\frac{b-a}2)^3\omega_0+(\frac{b-a}2)^3\omega_2=0 \]

而积分值也为0。

而对于$(x-\frac{a+b}2)^4$，数值积分公式有误差。

因此，公式是3阶代数精度。

\[\begin{pmatrix} 1 & 1 & 1 \\ -\frac{b-a}2 & 0 & \frac{b-a}2 \\ (\frac{b-a}2)^2 & 0 & (\frac{b-a}2)^2 \end{pmatrix} \begin{pmatrix} \omega_0 \\ \omega_1 \\ \omega_2 \end{pmatrix} =\begin{pmatrix} b-a \\ 0 \\ \frac{(b-a)^3}{12} \end{pmatrix} \]

插值多项式的积分

节点$x_i$上的插值多项式为$L_n(x)$，则它的积分也可以作为函数积分的近似

\[\int_a^b f(x)dx=\int_a^bL_n(x)dx+\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}w_n(x)dx \]

这样，得到了数值积分公式为

\[\begin{aligned} I_n(f)&=\int_a^bL_n(x)dx =\int_a^b\sum_{i=0}^nl_i(x)f(x_i)dx \\ &=\sum_{i=0}^n\int_a^bl_i(x)dxf(x_i) \end{aligned} \]

并且，误差为$\displaystyle\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}w_n(x)dx$

当节点$x_i$，积分区间$[a,b]$给定后，$\displaystyle\int_a^bl_i(x)dx$就是一个与原函数$f(x)$无关的量。

  • 取$\displaystyle a_i=\int_a^bl_i(x)dx$时，得到的数值积分公式就是插值多项式的积分。
  • 由误差表达式知，这个数值积分公式至少具有$n$阶代数精度。

若公式$I_n(f)=\displaystyle\sum_{i=0}^n\omega_if(x_i)$具有至少$n$阶代数精度，则有

\[\begin{cases} \omega_0(x_0)^0+\omega_1(x_1)^0+\cdots+\omega_n(x_n)^0&=\int_a^bx^0dx=(b-a) \\ \omega_0(x_0)^1+\omega_1(x_1)^1+\cdots+\omega_n(x_n)^1&=\int_a^bx^1dx=\frac{b^2-a^2}2 \\ %\omega_0(x_0)^2+\omega_1(x_1)^2+\cdots+\omega_n(x_n)^2&=\int_a^bx^2dx \\ \cdots \\ \omega_0(x_0)^n+\omega_1(x_1)^n+\cdots+\omega_n(x_n)^n&=\int_a^bx^ndx \\ \end{cases} \]

由Vandermonde行列式的特性，此时系数存在且唯一。这组系数就是

\[a_i=\int_a^bl_i(x)dx \]

这样，有如下结论

定理 1.
积分系数取为相应Lagrange插值基函数的积分的数值积分公式具有最高阶的代数精度。

此时，数值积分就是Lagrange插值多项式的积分，误差是插值多项式误差的积分。

Newton-Cotes 积分

在等距节点上得到的数值积分公式称为Newton-Cotes积分。

在等距节点上计算函数$f(x)$在$[a,b]$上积分，

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

其中$n\geq 1$。则节点$x_k$对应的积分系数为

\[\omega_k =\int_a^b \frac{(x-x_0)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)}{(x_k-x_0)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}dx \]

作变量代换$x=a+th$，则有

\[\begin{aligned} \omega_k=&\int_a^b \frac{(th)\cdots ((t-(k-1))h) ((t-(k+1))h)\cdots((t-n)h)} {(kh)\cdots(h) (-h)\cdots(-(n-k)h)}dx \\ =&\frac{(-1)^{n-k}h}{k!(n-k)!} I^{(n)}_k =nh\frac{(-1)^{n-k}}{n\cdot k!(n-k)!} I^{(n)}_k \end{aligned} \]

其中

\[I^{(n)}_k=\int_0^n t(t-1)\cdots(t-(k-1))(t-(k+1))\cdots(t-n)dt \]

注意到$nh=(b-a)$，记

\[\begin{aligned} C_k^{(n)}=\frac{(-1)^{n-k}}{n\cdot k!(n-k)!}I^{(n)}_k \end{aligned} \]

则积分系数可以表达为$\omega_k=(b-a) C_k^{(n)}$

得到$n\geq 1$时的Newton-Cotes积分公式为

\[\int_a^b f(x)dx\approx (b-a)\sum_{k=0}^n C_k^{(n)} f(x_k) \]

这个公式具有至少$n$阶代数精度。

$n=0$时，积分公式为

\[\int_a^b f(x)dx \approx (b-a)f(x_0) \]

因此，可以定义$C^{(0)}_0=1$。

定理 2.
$C^{(n)}_k$是一个只与$n$和$k$相关的数，具有如下特点

（1） $\displaystyle \sum_{k=0}^nC^{(n)}_k=1$

（2） $C^{(n)}_k=C^{(n)}_{n-k}$，即有对称性

(1) 取 $f(x)=1$ ，则有

\[\int_a^b 1 dx=\sum_{k=0}^nC^{(n)}_k\times 1 \]

(2) 令$s=n-t$，则有

\[\begin{aligned} I^{(n)}_k=&\int_0^n t(t-1)\cdots(t-(k-1))(t-(k+1))\cdots(t-n)dt \\ =&\int_n^0 (n-s)\cdots(n-s-(k-1))\\ &(n-s-(k+1))\cdots(n-s-n)d(n-s) \\ =&(-1)^n\int_0^n (s-n)\cdots(s-(n-k+1))(s-(n-k-1))\cdots s ds \\ =&(-1)^n\int_0^n t\cdots(t-(n-k-1))(t-(n-k+1))\cdots(t-n)dt \\ =&(-1)^n I^{(n)}_{n-k} \end{aligned} \]

从而，

\[\begin{aligned} C_k^{(n)}=&\frac{(-1)^{n-k}}{n\cdot k!(n-k)!}I^{(n)}_k \\ =&\frac{(-1)^{n-k}}{n\cdot k!(n-k)!}(-1)^nI^{(n)}_{n-k} \\ =&\frac{(-1)^{2n-k}}{n\cdot k!(n-k)!}I^{(n)}_{n-k} \\ =&\frac{(-1)^k}{n\cdot(n-k)!(n-(n-k))!}I^{(n)}_{n-k} \\ =&C_{n-k}^{(n)} \end{aligned} \]

n=1时，

\[\begin{aligned} C^{(1)}_0&=\frac{(-1)^1}{1\cdot 0!1!}\int_0^1t-1dt=\frac12 \\ C^{(1)}_1&=\frac{(-1)^0}{1\cdot 1!0!}\int_0^1tdt=\frac12 \end{aligned} \]

可以得到

定义 3.
梯形公式(trapezoid formula)

\[\int_a^bf(x)dx\approx \frac{b-a}{2}(f(a)+f(b)) \]

例 4. 分析梯形公式的误差

解. 数值积分的误差为插值多项式误差的积分，则有

\[R=\int_a^b\frac{f''(\xi)}{2!}(x-a)(x-b)dx \]

注意到$(x-a)(x-b)$在区间$[a,b]$内不变号，由积分中值定理

\[\begin{aligned} R&=\frac{f''(\eta)}{2!}\int_a^b(x-a)(x-b)dx \\ &=\frac{f''(\eta)}{2!}\frac{-1}6(b-a)^3 \\ &=-\frac{f''(\eta)}{12}(b-a)^3 , \xi\in(a,b) \end{aligned} \]

注. 梯形公式是由1次多项式的积分得到，具有至少1阶代数精度。这可以由误差表达式中看出。

也可以用Taylor展开来分析误差。

注意到所有奇数阶的积分为0，即

\[\begin{aligned} \int_a^b t^{2k-1} dx=&\int_a^b (x-\frac{a+b}2)^{2k-1}dx=0, \quad \forall k>0 \\ \int_a^b t^{2k} dx=&\int_a^b (x-\frac{a+b}2)^{2k}dx=2\frac{h^{2k+1}}{2k+1} \end{aligned} \]

取$m$为偶数，则积分可以展开为

\[\begin{aligned} \int_a^bf(x)=&(b-a)f(c)+\int_a^b\frac{f''(c)}{2!}t^2dx+\int_a^b\frac{f^{(4)}(c)}{4!}t^4dx \\ &+\cdots+\int_a^b\frac{f^{(m)}(\xi)}{m!}t^mdx \\ =&(b-a)f(c)+\frac{f''(c)}{2!}\frac23h^3+\frac{f^{(4)}(c)}{4!}\frac25h^5 \\ &+\cdots+\int_a^b\frac{f^{(m)}(\xi)}{m!}t^mdx \end{aligned} \]
\[f(a)=f(c)-f'(c)h+\frac{f''(c)}{2!}h^2-\frac{f'''(c)}{3!}h^3+\cdots+\frac{f^{(m)}(\xi_1)}{m!}(-h)^m \]
\[f(b)=f(c)+f'(c)h+\frac{f''(c)}{2!}h^2+\frac{f'''(c)}{3!}h^3+\cdots+\frac{f^{(m)}(\xi_2)}{m!}h^m \]

取$m$为偶数，梯形公式可以展开为

\[\begin{aligned} \frac{b-a}2(f(a)+f(b))=&(b-a)f(c)+2h\frac{f''(c)}{2!}h^2+2h\frac{f^{(4)}(c)}{4!}h^4\\ &+\cdots+h\frac{f^{(m)}(\xi_1)+f^{(m)}(\xi_2)}{m!}h^m \end{aligned} \]

可以得到梯形公式的误差为

\[R=\frac{f''(c)}{2!}2h^3(\frac13-1)+\frac{f''(c)}{4!}2h^3(\frac15-1)+\cdots \]

n=2时，

\[\begin{aligned} C^{(2)}_0&=\frac{(-1)^2}{2\cdot 0!2!}\int_0^2(t-1)(t-2)dt=\frac16 \\ C^{(2)}_1&=\frac{(-1)^1}{2\cdot 1!1!}\int_0^2t(t-2)dt=\frac46 \\ C^{(2)}_2&=\frac{(-1)^0}{2\cdot 2!0!}\int_0^2t(t-1)dt=\frac16 \\ \end{aligned} \]

可以得到

定义 4.
Simpson公式

\[\int_a^bf(x)dx\approx \frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b)) \]

例 5. Simpson公式的误差？

解. Simpson公式是由2次多项式的积分得到，具有至少2阶代数精度。

  • 对于$x^3$，它的Simpson公式计算结果为

    \[\frac{b-a}{6}(a^3+4(\frac{a+b}{2})^3+(b)^3)=\frac{b^4-a^4}{4} \]

  • 而对于$x^4$，Simpson公式计算结果有误差。

因此，Simpson公式具有3阶代数精度。

记$I(f)$为函数$f(x)$的积分值， $S(f)$为Simpson公式计算的值，则误差

\[\begin{aligned} R&=I(f)-S(f) \\ &=I(f)-I(P_3)+I(P_3)-S(f) \\ &=I(f)-I(P_3)+S(P_3)-S(f) \\ \end{aligned} \]

当$P_3$是个3次多项式，且是函数$f(x)$在点$a$, $\frac{a+b}2$, $b$上的插值时，进一步有

\[R=I(f)-I(P_3) \]

显然这样的3次多项式有无穷多个，找一个好用的。

增加中点的导数条件$f'(\frac{a+b}2)$，则有

\[R=I(f-P_3)=\int_a^b\frac{f^{(4)}(\xi)}{4!}(x-a)(x-\frac{a+b}2)^2(x-b)dx \]

由积分中值定理，有

\[\begin{aligned} R&=\frac{f^{(4)}(\eta)}{4!}\int_a^b(x-a)(x-\frac{a+b}2)^2(x-b)dx \\ &=-\frac{(b-a)^5}{2880}f^{(4)}(\eta) \end{aligned} \]

可以从误差表达式中看到，Simpson公式具有3阶代数精度。

同样，可以用Taylor展开分析误差。

在$c=\frac{a+b}2$中，做Taylor展开，则Simpson公式可以展开为

\[\begin{aligned} &\frac{b-a}6(f(a)+f(b)+4f(c)) \\ =&(b-a)f(c)+\frac{f''(c)}{2!}2h^3\frac13+\frac{f^{(4)}(c)}{4!}2h^5\frac13+\cdots \end{aligned} \]

因而，误差为

\[R=\frac{f^{(4)}(c)}{4!}2h^5(\frac13-\frac15) +\frac{f^{(6)}(c)}{6!}2h^7(\frac13-\frac17)+\cdots \]

n=4时，可以得到Boole公式

\[ \frac{b-a}{90}(7f(a)+32f(x_1)+12f(x_2)+32f(x_3)+7f(b)) \]

其中$x_i=a+i\frac{b-a}4, i=1,2,3$

注. Boole公式具有5阶代数精度。它的误差是$-\frac{(b-a)^7}{1935360}f^{(6)}(\xi)$

定理 3.
对于$n\leq7$，所有的系数$C^{(n)}_k$都是正的。这些Newton-Cotes积分格式都是稳定的。

证明. 假定每个节点带有误差$\delta_k$，$k=0,1,\cdots,n$，则有

\[\begin{aligned} \epsilon=&I_n(\tilde f)-I_n(f) \\ =&(b-a)\sum_{k=0}^nC^{(n)}_k(f(x_k)+\delta_k)-(b-a)\sum_{k=0}^nC^{(n)}_k(f(x_k)) \end{aligned} \]

记$\delta=\max|\delta_k|$，则有

\[\begin{aligned} |\epsilon|=&\left|(b-a)\sum_{k=0}^nC^{(n)}_k\delta_k\right| \\ \leq&(b-a)\sum_{k=0}^n\left|C^{(n)}_k\delta_k\right| \\ =&(b-a)\sum_{k=0}^n C^{(n)}_k\left|\delta_k\right| \\ \leq&(b-a)\delta\sum_{k=0}^n C^{(n)}_k = (b-a)\delta \end{aligned} \]

即，格式是稳定的。

定理 4.
当$n$为偶数时，Newton-Cote公式至少有$n+1$阶代数精度

证明. 由数值积分公式的误差表达式

\[R_n=\frac{1}{(n+1)!}\int_a^b f^{(n+1)}(\xi)\omega_n(x)dx \]

其中

\[\begin{aligned} \omega_n(x)=(x-x_0)(x-x_1)\cdots(x-x_n) \\ x_i=a+i h, h=\frac{b-a}n, i=0,1,2,\cdots,n \end{aligned} \]

做变量代换$x=a+th$，则有

\[R_n=\frac{h^{n+2}}{(n+1)!}\int_0^n f^{(n+1)}(\xi)\prod_{j=0}^n(t-j)dt \]

由

\[R_n=\frac{h^{n+2}}{(n+1)!}\int_0^n f^{(n+1)}(\xi)\prod_{j=0}^n(t-j)dt \]
  • 当$f(x)$为不高于$n$次的多项式时，$f^{(n+1)}(\xi)=0$，从而误差为$0$。
  • 当$f(x)=x^{n+1}$时，$f^{(n+1)}(\xi)=(n+1)!$，则

    \[R_n=h^{n+2}\int_0^nt(t-1)\cdots(t-n)dt \]

当$n$为偶数时，令$t=u+\frac{n}2$，则

\[\begin{aligned} R_n=h^{n+2}\int_{-\frac{n}2}^{\frac{n}2} & (u+\frac{n}2)(u+\frac{n}2-1)\cdots(u+1) \\ &u(u-1)\cdots(u-\frac{n}2+1)(u-\frac{n}2)du \\ =h^{n+2}\int_{-\frac{n}2}^{\frac{n}2} & u(u^2-1)(u^2-2^2)\cdots(u^2-\frac{n^2}4)du \end{aligned} \]

注意到被积函数是奇函数，且积分区间关于$0$对称。

因而当$n$为偶数时，对$n+1$次多项式的误差为$0$。

定理 5.
Newton-Cotes积分公式的误差为

(1)$n$为奇数时，$f\in C^{n+1}[a,b]$，有

\[R_n(f)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\int_a^bw(x)dx \]

(2)$n$为偶数时，$f\in C^{n+2}[a,b]$，有

\[R_n(f)=\frac{f^{(n+2)}(\xi)}{(n+2)!}\int_a^bxw(x)dx \]

目录

本节读完

例 6.

6.