七、数值微分与数值积分

1 数值微分

1.1 差商型求导公式

待补充

两点公式（n=1）

过节点 $x_0,x_1=x_0+h$ 的Lagrange插值多项式为：

$$L_1(x)=\frac{x-x_1}{-h}f(x_0)+\frac{x-x_0}{h}f(x_1)$$

故有：

$$\{\begin{array}{l}{f}^{\prime }(x_0)\approx L_1(x_0)=\frac{f(x_1)-f(x_0)}{h}\\ f^{\prime }(x_1)\approx L_1^{\prime }(x_1)=\frac{f(x_1)-f(x_0)}{h}\end{array}(7-9)$$

截断误差为：

待补充

三点公式

过节点 $x_i=x_0+ih(i=0,1,2)$ , $f(x)$ 的Lagrange插值多项式为：

$$L_2(x)=\frac{(x-x_1)(x-x_2)}{2h^2}f(x_0)+\frac{(x-x_0)(x-x_2)}{-h^2}f(x_1)f(x_2)+\frac{(x-x_0)(x-x_1)}{2h^2}f(x_2)$$

求导得

$$$$

1.2 插值型求导公式

待补充

下面给出等距节点情况下，求一阶导数的几个常用公式。

两点公式（n=1）

过节点 $x_0,x_1=x_0+h$ 的Lagrange插值多项式为：

$$L_1(x)=\frac{x-x_1}{-h}f(x_0)+\frac{x-x_0}{h}f(x_1)$$

故有：

$$\{\begin{array}{l}{f}^{\prime }(x_0)\approx L_1^{\prime }(x_0)=\frac{f(x_1)-f(x_0)}{h}\\ f^{\prime }(x_1)\approx L_1^{\prime }(x_1)=\frac{f(x_1)-f(x_0)}{h}\end{array}(7-9)$$

截断误差（？）为：

$$\{\begin{array}{l}{R}_1^{\prime }(x_0)=-\frac{h}{2}{f}^{\prime \prime }({\xi }_0)\\ R_1^{\prime }(x_1)=\frac{h}{2}{f}^{\prime \prime }({\xi }_1)\end{array}({\xi }_0,{\xi }_1\in (x_0,x_1)(7-10)$$

三点公式

过节点$x_i=x_0+ih(i=0,1,2),f(x)$的Lagrange插值多项式（？）为：

待补充

分别代入$x_0,x_1,x_2$，得三点公式：

$$\{\begin{array}{l}{f}^{\prime }(x_0)\approx \frac{1}{2h}[-3f(x_0)+4f(x_1)-f(x_2)]\\ f^{\prime }(x_1)\approx \frac{1}{2h}[-f(x_0)+f(x_2)]\\ f^{\prime }(x_2)\approx \frac{1}{2h}[f(x_0)-4f(x_1)+3f(x_2)]\end{array}(7-11)$$

待补充

2 牛顿-科茨（Newton-Cotes）求积公式

2.1 数值积分的基本思想

定积分的定义表明，函数$f(x)$在$[a,b]$上的积分值是和式的极限：

$$I(f)={\int }_a^{b}f(x)\text{d}x=\underset{\lambda \to 0}{\lim }\sum _{i=1}^{n}f({\xi }_i)\Delta x_i$$

其中$\lambda =\underset{1\le i\le n}{\max }\Delta x_i$，积分和式$\sum _{i=1}^{n}f({\xi }_i)\Delta x_i$是一些点上函数值的线性组合。构造数值积分方法的基本思想就是用被积函数在积分区间$[a,b]$上的某些节点$x_k(k=0,1,\cdots ,n)$处的函数值的线性组合作为定积分的近似值，即：

$$I(f)={\int }_a^{b}f(x)\text{d}x\approx \sum _{k=0}^n{A}_{k}f(x_k)(7-16)$$

其中{Ak}\{A_k\}{Ak​}称为求积系数，它只与节点{xk}\{x_k\}{xk​}有关，与被积函数$f(x)$无关。

由节点{xk}\{x_k\}{xk​}和求积系数{Ak}\{A_k\}{Ak​}的不同取法，可以得到不同的求积公式。最直接自然的一种想法是用$f(x)$在$[a,b]$上的插值多项式${\phi }_n(x)$代替$f(x)$，以${\phi }_n(x)$在区间$[a,b]$上的积分值作为所求积分$I(f)$的近似值，即：

$$I(f)\approx {\int }_a^b{\phi }_n(x)\text{d}x$$

这样得到的求积公式称为插值型求积公式，通常采用Lagrange插值。

设 $[a,b]$ 上有 $n+1$ 个互异节点 $x_0,x_1,\dots ,x_n$，$f(x)$ 的 $n$ 次Lagrange插值多项式为：

$$L_N(x)=\sum _{k=0}^n{l}_k(x)f(x_k)$$

其中$l_k(x)=\prod _{j=0,j\ne k}^n\frac{x-x_i}{x_k-x_j}$，插值型求积公式为：

$$I(f)\approx {\int }_a^b{L}_n(x)\text{d}x=\sum _{k=0}^n({\int }_a^b{l}_k(x)\text{d}x)f(x_k)(7-17)$$

显然，求积系数$A_k={\int }_a^b{l}_k(x)\text{d}x(k=0,1,\cdots ,n)$仅与节点{xk}\{x_k\}{xk​}有关，与被积函数$f(x)$的具体形式无关。求积公式（7-17）

待补充

2.2 Newton-Cotes公式

如果取等距节点{xk}\{x_k\}{xk​}，则插值型求积公式更简便易算：

设$x_k=a+kh(k=0,1,\cdots ,n),h=\frac{b-a}{n}$，由式（7-17）：

$$A_k={\int }_a^b{l}_k(x)\text{d}x={\int }_a^b(\prod _{j=0,j\ne k}^n\frac{x-x_j}{x_k-x_j})\text{d}x$$

待补充

令

$$C_K^{(n)}=\frac{(-1{)}^{n-k}}{k!(n-k)!n}{\int }_0^n\prod _{j=0,j\ne k}^n(t-j)dt(7-19)$$

则$A_k=(b-a)C_k^{(n)}$，求积公式（7-17）可简化为：

$$I(f)={\int }_a^{b}f(x)d_x\approx (b-a)\sum _{k=0}^n{C}_K^{(n)}f(x_k)(7-20)$$

式（7-20）称为n阶Newton-Cotes公式，简记为N-C公式，{Ck(n)}\{C_k^{(n)}\}{Ck(n)​}称为Cotes系数。

N-C公式的截断误差为：

$$\begin{gathered} R_n(f)={\int }_a^b\frac{f^{(n+1)}(\xi )}{(n+1)!}\prod _{j=0}^n(x-x_j)dx \\ =\frac{h^{n+2}}{(n+1)!}{\int }_0^n{f}^{(n+1)}(\xi )\prod _{j=0}^n(t-j)dt(7-21) \end{gathered}$$

因为$C_k^n$仅与等分区间数n有关，与积分空间及被积函数均无关，故可事先造好Cotes系数表（？表7-2），这样使用N-C公式计算定积分近似值更加方便。例如，n=4时，N-C公式为：

待补充

式（7-22）是常用的求积公式，称为Cotes公式。

特别地，当n=1时，$C_0^{(1)}=C_1^{(1)}=\frac{1}{2}$，因此有：

$$I(f)={\int }_a^{b}f(x)\text{d}x\approx \frac{b-a}{2}[f(a)+f(b)](7-23)$$

这是以过点$(a,f(a)),(b,f(b))$的直线代替曲线$y=f(x)$，以梯形面积近似曲边梯形面积。故式（7-23）又称为梯形公式。

当n=2时，N-C公式为：

$$I(f)={\int }_a^{b}f(x)\text{d}x=\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)](7-24)$$

这是以过曲线上三点$(a,f(a)),(\frac{a+b}{2},f(\frac{a+b}{2})),(b,f(b))$的抛物线代替曲线$y=f(x)$求曲边梯形面积的近似值，因而此公式称为抛物线求积公式，又称辛普森（Simpson）公式。

2.3 误差估计

首先引入一种衡量数值积分公式近似程度的概念。

求积公式的代数精确度

定义7.1

若当$f(x)$为任意次数不高于m的多项式时，求积公式

$${\int }_a^{b}f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)$$

均精确成立，而对某个m+1次多项式，公式不精确成立，则称该求积公式具有m次代数精确度。

可以验证，梯形公式具有一次代数精确度。因为若$f(x)$是一次多项式，则：

$$R_1(x)=f(x)-L_1(x)=\frac{f^{\prime \prime }(\xi )}{2}(x-a)(x-b)=0$$

从而有：

$${\int }_a^{b}f(x)\text{d}x={\int }_a^b{L}_1(x)dx=\frac{b-a}{2}[f(a)+f(b)]$$

若取$f(x)=x^2$，则：

$${\int }_a^{b}f(x)dx={\int }_a^b{x}^{2}dx=\frac{1}{3}(b^3-a^3)\ne \frac{b-a}{2}(a^2+b^2)$$

一般地，由n次插值多项式导出的求积公式至少有n次代数精确度。对N-C公式，更有以下结论：

定理7.1

2n阶Newton-Cotes公式至少具有2n+1次代数精确度。

容易证明，求积公式（7-16）具有m次代数精确度的充要条件是它对$f(x)=1,x,\cdots ,x^m$都能精确成立，但对$f(x)=x^{m+1}$，式（7-16）不准确成立。这一结论给出了判别一个求积公式的代数精确度的方法。

待补充 215

利用余项公式及积分中值定理，可以导出低阶N-C公式的误差估计。

梯形公式和辛普森公式的误差估计

由式（7-21），梯形公式的截断误差为：

$$R_1(f)={\int }_a^b\frac{f^{\prime \prime }({\xi }_x)}{2}(x-a)(x-b)\text{d}x,{\xi }_x\in (a,b)$$

因为(x-a)(x-b)在区间$[a,b]$内不变号，如果$f(x)$二阶连续可导，则由反常积分中值定理（？），存在$\xi \in (a,b)$，使得：

$$\begin{gathered} R_1(f)=\frac{f^{\prime \prime }(\xi )}{2}{\int }_a^b(x-a)(x-b)\text{d}x \\ =\frac{f^{\prime \prime }(\xi )}{2}[\frac{x^3}{3}-(a+b)\frac{x^2}{2}+abx{]}_a^b \\ =-\frac{f^{\prime \prime }(\xi )}{12}(b-a{)}^3 \end{gathered}$$

定理7.2

（？）这个C是啥意思

若$f(x)\in C^{(2)}[a,b]$，则梯形公式的误差为：

$$R_1(f)={\int }_a^{b}f(x)\text{d}x-\frac{b-a}{2}[f(a)+f(b)]=-\frac{(b-a{)}^3}{12}{f}^{\prime \prime }(\xi ),\xi \in (a,b)(7-25)$$

借助于$f(x)$的Hermite插值与定理7.1（？），可导出Simpson公式的误差估计：

定理7.3

若$f(x)\in C^{(4)}[a,b]$，则Simpson公式的误差为：

$$R_2(f)={\int }_a^{b}f(x)\text{d}x-\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]=-\frac{(b-a{)}^5}{2880}{f}^{(4)}(\xi )=-\frac{h^5}{90}{f}^{(4)}(\xi ),\xi \in (a,b)(7-26)$$

其中$h=\frac{(b-a)}{2}$为步长。

$C^{(2)}[a,b]$应该是对$[a,b]$二等分的含义。

3 复化求积公式

前面导出的误差估计式表明，用N-C公式计算积分近似值时，步长越小，截断误差越小（？）

待补充

3.1 复化梯形公式

用分段线性插值函数近似被积函数，等于把积分区间分成若干小区间，在每个小区间上以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。这样求得的近似值显然比用梯形公式计算精度高。定积分存在定理表明，只要被积函数连续，当小区间长度趋于零时，小梯形面积之和趋于曲边梯形面积的准确值，即定积分的准确值。

将积分区间 $[a,b]$ 分成 $n$ 等份，记$h=\frac{b-a}{n},x_k=a+kh(k=0,1,\cdots ,n)$，在每个小区间$[x_k,x_{k+1}](k=0,1,\cdots ,n-1)$上用梯形公式并求和，得：

$${\int }_a^{b}f(x)\text{d}x=\sum _{k=0}^{n-1}{\int }_k^{k+1}f(x)\text{d}x\approx \sum _{k=0}^{n-1}\frac{h}{2}[f(x_k)+f(x_{k+1})]$$

整理得：

$${\int }_a^{b}f(x)\text{d}x\approx \frac{h}{2}[f(a)+f(b)+2\sum _{k=1}^{n-1}f(x_k)]=T_n(7-28)$$

式（7-28）称为复化梯形公式。

如果$f(x)\in C^{(2)}[a,b]$（？），在小区间$[x_k,x_{k+1}]$上，梯形公式的截断误差为：

$${\int }_{x_k}^{x_{k+1}}dx-\frac{h}{2}[f(x_k)+f(x_{k+1})]=-\frac{h^3}{12}{f}^{\prime \prime }({\xi }_k),{\xi }_k\in (x_k,x_{k+1})$$
？？？？ 12 是哪来的

因此

$$R_T(f)={\int }_a^{b}f(x)dx-T_n=-\frac{h^3}{12}\sum _{k=0}^{n-1}{f}^{\prime \prime }({\xi }_k)$$

因为$f^{\prime \prime }(x)$在[a,b]连续，由介值定理，存在$\xi \in (a,b)$，使得：

$$f^{\prime \prime }(\xi )=\frac{1}{n}\sum _{k=0}^{n-1}{f}^{\prime \prime }({\xi }_k)$$

从而有：

$$R_T(f)={\int }_a^{b}f(x)dx-T_n=-\frac{h^3}{12}n{f}^{\prime \prime }(\xi )=-\frac{b-a}{12}{h}^2{f}^{\prime \prime }(\xi ),\xi \in (a,b)(7-29)$$

这就是复化梯形公式的截断误差。

3.2 复化Simpson公式

如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就导出复化Simpson公式。

将区间[a,b]分成n=2m等份，分点为$x_k=a+kh(k=0,1,\cdots ,n),h=\frac{b-a}{n}$，在每个小区间$[x_{2k-2},x_{2k}]$上，用Simpson公式求积分，则有：

待补充

整理后得到：

$${\int }_a^{b}f(x)\text{d}x\approx \frac{h}{3}[f(a)+f(b)+2\sum _{k=1}^{m-1}f(x_{2k})+4\sum _{k=1}^{m}f(x_{2k-1})](7-30)$$

式（7-30）称为复化Simpson公式。

4 龙贝格（Romberg）求积公式

4.1 理查逊（Richardson）外推法

设用某种数值方法求量$I$的近似值，一般地，近似值是步长$h$的函数，记为$I_1(h)$，相应的误差为：

$$I-I_1(h)=a_1{h}^{p_1}+a_2{h}^{p_2}+\cdots +a_k{h}^{p_k}+\cdots (7-36)$$

其中$0<p_1<p_2<\cdots <p_k<\cdots ,a_i(i=1,2,\cdots )$与h无关。若用$\alpha h$代替式（7-36）中的h，则得：

$$I-I_1(\alpha h)=a_1{\alpha }^{p_1}{h}^{p_1}+a_2{\alpha }^{p_2}{h}^{p_2}+\cdots +a_k{\alpha }^{p_k}{h}^{p_k}+\cdots (7-37)$$

式（7-37）减去式（7-36）乘以${\alpha }^{p_1}$，得：

$$\begin{gathered} I-I_1(\alpha h)-{\alpha }^{p_1}[I-I_1(h)] \\ =a_2({\alpha }^{p_2}-{\alpha }^{p_1})h^{p_2}+a_3({\alpha }^{p_3}-{\alpha }^{p_1})h^{p_3}+\cdots +a_k({\alpha }^{p_k}-{\alpha }^{p_1})h^{p_k}+\cdots \end{gathered}$$

取$\alpha$满足$|\alpha |\ne 1$，以$1-{\alpha }^{p_1}$除上式两端，得：

$$I-\frac{I_1(\alpha h)-{\alpha }^{p_1}{I}_1(h)}{1-{\alpha }^{p_1}}=b_2{h}^{p_2}+b_3{h}^{p_3}+\cdots +b_k{h}^{p+k}+\cdots (7-38)$$

其中$b_i=a_i({\alpha }^{p_i}-{\alpha }^{p_1})/(1-{\alpha }^{p_1})(i=2,3,\dots )$仍与h无关。令：

$$I_2(h)=\frac{I_1(\alpha h)-{\alpha }^{p_1}{I}_1(h)}{1-{\alpha }^{p_1}}$$

由式（7-38），以$I_2(h)$作为I的近似值，其误差至少为$O(h^{p_2})$（？），因此$I_2(h)$收敛于$I$的速度比$I_1(h)$快。不断重复以上做法，可以得到一个函数序列：

$$I_m(h)=\frac{I_{m-1}(\alpha h)={\alpha }^{P_{m-1}}{I}_{m-1}(h)}{1-{\alpha }^{P_{}m-1}}(m=2,3,\cdots )$$

以$I_m(h)$近似I，误差为$I-I_m(h)=O(h^{P_m})$。随着m的增大，收敛速度越来越快，这就是Richardson外推法。将此方法用于梯形值序列{T2k}\{T_{2^k}\}{T2k​}，可以构造出一种收敛很快，计算简便的数值积分方法。

4.2 Romberg求积公式

设 $I={\int }_a^{b}f(x)\text{d}x$ ，记 $2^k$ 等分区间 $[a,b]$ 用复化梯形公式所求得的近似值为 $T_{2^k}=T_0^{(k)}=I_1(\frac{b-a}{2^k})(k=0,1,2,\dots )$ 。由Euler-Maclaurin公式

$$R(f,T_0^{(k)})=a_1(\frac{b-a}{2^k}{)}^2+a_2(\frac{b-a}{2^k}{)}^4+\cdots +a_i(\frac{b-a}{2^k}{)}^{2i}+\cdots$$

取$\alpha =\frac{1}{2}$，由Richardson外推公式：

$$I_2(\frac{b-a}{2^k})=\frac{I_1(\frac{b-a}{2^{k+1}})-(\frac{1}{2}{)}^2{I}_1(\frac{b-a}{2^k})}{1-(\frac{1}{2}{)}^2}=\frac{4T_0^{(k+1)}-T_0^{(k)}}{3}$$

记$I_2=(\frac{b-a}{2^k})=T_1^{(k)}$，则：

$$T_1^{(k)}=\frac{4T_0^{(k+1)}-T_0^{(k)}}{3}(k=0,1,2\cdots )(7-39)$$

且有：

$$R(f,T_1^{(k)})=O((\frac{b-a}{2^k}{)}^4)$$

事实上，容易验证，{T1(k)}\{T_1^{(k)}\}{T1(k)​}就是Simpson值序列{S2k}\{S_{2^k}\}{S2k​}

待补充

记$T_m^{(k)}=I_{m+1}(\frac{b-a}{2^k})$，则上式可记为：

$$T_m^{(k)}=\frac{4^m{T}_{m-1}^{(k+1)}-T_{m-1}^{(k)}}{4^m-1}(m=1,2,\cdots ;k=0,1,2,\cdots )(7-40)$$

5 高斯（Gauss）型求积公式

5.1 一般理论

Gauss型求积公式

Gauss型求积公式的基本思想是：在节点数n固定时，适当地选取节点{xk}\{x_k\}{xk​}与求积系数{Ak}\{A_k\}{Ak​}，使求积公式具有最高的代数精确度。

由定义（7.1），求积公式

$${\int }_a^b\rho (x)f(x)\text{d}x\approx \sum _{k=1}^n{A}_{k}f(x_k)(7-41)$$

其中$\rho (x)\ge 0$为权函数，具有m次代数精确度，充要条件是对$f(x)=x^l(l=0,1,\cdots ,m)$式（7-41）精确成立，即：

$${\int }_a^b\rho (x)x^l\text{d}x\approx \sum _{k=1}^n{A}_k{x}_k^l(7-42)$$

令：

$${\mu }_l={\int }_a^b\rho (x)x^l\text{d}x(l=0,1,\cdots ,m)(7-43)$$

将式（7-43）代入式（7-42），得节点{xk}\{x_k\}{xk​}与求积系数{Ak}\{A_k\}{Ak​}应满足方程组：

$$\{\begin{array}{l}{A}_1+A_2+\cdots +A_n={\mu }_0\\ A_1{x}_a+A_2{x}_2+\cdots +A_n{x}_n={\mu }_1\\ A_1{x}_1^m+A_2{x}_2^m+\cdots +A_n{x}_n^m={\mu }_m\end{array}(7-44)$$

式（7-44）含有2n个未知数，m+1个方程，可以证明，当$m+1\le 2n$时，方程组（7-44）有解。这表明n个节点的求积公式的代数精确度可达到2n-1。但若取2n次多项式：

$$P_{2n}(x)=(x-x_1{)}^2(x-x_2{)}^2\cdots (x-x_n{)}^2$$

显然，对任意n个互异节点$x_1,x_2,\cdots ,x_n$，都有$P_{2n}(x_k)=0(k=1,2,\cdots ,n)$，从而有$\sum _{k=1}^n{P}_{2n}(x_k)=0$。然而：

$${\int }_a^b\rho (x)P_{2n}(x)\text{d}x>0$$

因此，n个节点的求积公式的最高代数精确度为2n-1（？）

定义 7.2

如果一组节点$x_1,x_2,\cdots ,x_n\in [a,b]$能使求积公式：

$${\int }_a^b\rho (x)f(x)\text{d}x\approx \sum _{k=1}^n{A}_{k}f(x_k)(7-41)$$

具有2n-1次代数精确度，则称这组节点{xk}\{x_k\}{xk​}为Gauss点，公式（7-41）称为带权函数$\rho (x)$的Gauss型求积公式。

Gauss点与正交多项式的关系