数值分析 笔记 7

数值微分与数值积分

数值积分

指的是用数值方法计算定积分，在一些解析解非常复杂，或根本没有解析解的情况下使用

根据定积分的定义：

$$I(f) = \int_a^b f(x)dx = \lim\limits_{n \to \infty}\sum\limits_{i = 0}^{n}(x_{i + 1} - x_{i})f(\xi_i)$$

则一种近似是取一个充分大的 $n$，计算：

$$I_n(f) = \sum\limits_{i = 0}^{n}A_if(x_i)$$

其中 $A_i$ 称为积分系数，$x_k$ 称为积分节点

插值型求积公式

使用多项式函数 $p(x)$ 来插值原函数 $f(x)$，则将 $I(f)$ 转化为了易于计算的多项式积分形式

具体来说，在区间 $[a, b]$ 上取 $n + 1$ 个插值点，之后使用多项式进行 Lagrange 插值，则求积公式为：

$$I_n(f) = \sum\limits_{k = 0}^{n}\Big[f(x_k)\int_a^bl_k(x)dx\Big]$$

当 $n = 0, 1$ 时：

$$\begin{align*} I_0(f) &= (b - a)f(\frac{a + b}{2}) \\ I_1(f) &= \frac{b - a}{2}[f(a) + f(b)] \end{align*}$$

数值积分的性质

积分余项

定义积分余项为：

$$R_n(x) = I(f) - I_n(f)$$

对于插值型求积公式，有：

$$R_n(x) = \int_a^b \frac{f^{(n + 1)}(\xi)}{(n + 1)!}\omega_{n + 1}(x)dx$$

代数精度

定义一个求积公式的代数精度为：如果 $I_n(f)$ 在 $f$ 为次数不大于 $m$ 时均准确，在 $f$ 为次数 $m + 1$ 时不准确，则称该公式的代数精度为 $m$

可以证明的是：

两个推论：

1. 由于代数精度不低于 $0$，因此取 $f(x) = 1$，则有：

$$\sum\limits_{k = 0}^{n}A_k = b - a$$

2. $I_n(f)$ 至少有 $n$ 阶代数精度当且仅当 $I_n(f)$ 是插值型求积公式

数值性质

若对于求积公式 $I_n(f)$，有：

$$\lim\limits_{n \to\infty}I_n(f) = I(f)$$

则称该求积公式是收敛的

考虑对积分的输入进行微扰，$f(x) \to \tilde{f}(x)$，则扰动大小为：

$$\delta = ||f(x) - \tilde{f}(x)||_\infty = \max\limits_{a\leq x \leq b}|f(x) - \tilde{f}(x)|$$

因此结果的误差为：

$$|I(f) - I(\tilde{f})| = (b - a)|f(\xi) - \tilde{f}(\xi)| \leq (b - a)\delta$$

这说明求积分的问题一般是不敏感的

而使用求积公式时，输入进行微扰会导致：

$$\begin{align*} |I_n(f) - I_n(\tilde{f})| &= \Big|\sum\limits_{k = 0}^{n}A_k\big[f(x_k) - \tilde{f}(x_k)\big]\Big| \\ &\leq \sum\limits_{k = 0}^{n}\big|A_k\big|\big|f(x_k) - \tilde{f}(x_k)\big| \\ &\leq \varepsilon\sum\limits_{k = 0}^{n}\big|A_k\big| \\ \end{align*}$$

其中 $\varepsilon = \max\limits_{0\leq k \leq n}|f(x_k) - \tilde{f}(x_k)| \leq \delta$

如果我们确定所有的积分系数都为正，则：

$$\begin{align*} |I_n(f) - I_n(\tilde{f})| \leq \varepsilon(b - a) \\ \end{align*}$$

也即在 $A_k \geq 0$ 的情况下，求积公式是稳定的

Newton-Cotes 公式

Newton-Cotes 是插值型求积公式的一种，将一个区间 $[a, b]$ 均匀划分为 $n$ 个小区间，其积分系数为：

$$A_k = \int_a^b \Bigg[ \prod\limits_{\substack{j = 0 \\j \neq k}}^{n}\frac{(x - x_j)}{(x_k - x_j)}\Bigg] dx$$

令 $x = a + th$，其中 $h = \dfrac{b - a}{n}$，则积分系数为：

$$A_k = \int_0^n \Bigg[ \prod\limits_{\substack{j = 0 \\j \neq k}}^{n}\frac{(t - j)}{(k - j)}\Bigg]\frac{b - a}{n} dt = (b - a)C_k^{(n)}$$

其中 $C_k^{(n)}$ 是与积分区间无关的系数，被称为 Cotes 系数，低阶的 Cotes 系数表为：

注：n = 1 的时候即为梯形公式，且一般不使用 $n = 3$ 的 N-C 公式

几个特例：

1. 梯形公式，$n = 1$ 时的 N-C 公式

   $$T(f) = (b - a)\big[\frac{1}{2}f(a) + \frac{1}{2}f(b)\big]$$

   其余项为：

   $$R_T = \int_a^b \frac{f''(\xi)}{2!}(x - a)(x - b)dx = -\frac{f''(\eta)}{12}(b - a)^{3}$$

2. Simpson 公式，$n = 2$ 时的 N-C 公式

   $$S(f) = (b - a)\big[\frac{1}{6}f(a) + \frac{2}{3}f(\frac{a + b}{2}) + \frac{1}{6}f(b)\big]$$

   其余项为：

   $$R_T = \int_a^b \frac{f^{(4)}(\xi)}{4!}(x - a)(x - \frac{a + b}{2})(x - b)dx = -\frac{f^{(4)}(\eta)}{2880}(b - a)^{5}$$

3. Cotes 公式，$n = 4$ 时的 N-C 公式

   $$C(f) = (b - a)\big[\frac{7}{90}f(x_0) + \frac{16}{45}f(x_1) + \frac{2}{15}f(x_2) + \frac{16}{45}f(x_3) + \frac{7}{90}f(x_4)\big]$$

但是 $n = 8$ 或 $n \geq 10$ 的时候，Cotes系数表中出现了负系数，这代表此时的求积公式不稳定，并且由于 Runge 现象，$n$ 增大并不会代表收敛性一定会增强，因此高阶 N-C 公式是不稳定、不一定准确的

可以证明，当 $n$ 为偶数的时候，$n$ 阶 N-C 公式至少有 $n + 1$ 阶的代数精度

因此，综合考虑计算量、稳定性与准确性，通常情况下只会使用 $1, 2, 4, 6$ 阶四种 N-C 公式

复合求积公式

类似全局插值和分段插值的思路，既然 N-C 公式在高阶的时候很多方面的表现都很差，那就将全局区间划分成若干个小区间，之后在每一个小区间上用低阶的 N-C 公式即可，这种方法被称为复合求积公式

复合梯形公式

将 $[a, b]$ 划分为 $n$ 个小区间，则：

$$\int_{x_k}^{x_{k + 1}}f(x)dx \approx T_{n}^{(k)} = \frac{b - a}{2n}[f(x_k) + f(x_{k + 1})]$$

在全区间上：

$$T_n = \sum\limits_{k = 0}^{n - 1}\Bigg(\int_{x_k}^{x_{k + 1}}f(x)dx\Bigg) = \frac{b - a}{2n}\Big[f(a) + 2\sum\limits_{k = 1}^{n - 1}f(x_k) + f(b)\Big]$$

这被称为复合梯形公式

每个区间上的截断误差，即求积余项为：

$$R_T^{(k)} = -\frac{f''(\eta)}{12}\big(\frac{b - a}{n}\big)^{3}$$

则复合梯形公式的截断误差为：

$$R_T = \sum\limits_{k = 0}^{n - 1}R_{T}^{(k)} = -\frac{f''(\eta)}{12n^{2}}\big(b - a)^{3} = O(h^{2})$$

随着 $n$ 的增加，即划分的精细化，误差会显著降低

对于等距节点求积公式，若其截断误差为 $O(h^p)$，则称其准确度阶数为 $p$，其中 $h$ 为两个节点之间的距离

则根据上述推导，复合梯形公式具有 2 阶准确度

同样，复合梯形公式还有几条性质：

• $T_n$ 具有 1 阶代数精度
• 对于任意可积函数 $f$，$T_n$ 是收敛的

由于我们很难预知 $n$ 的值，因此一般采取动态折半的方式，每次在区间折半的过程中，之前的所有点仍然是区间端点，计算也可以复用

记 $x_{k + 1/2} = \dfrac{x_k + x_{k + 1}}{2}$，则：

$$\begin{align*} T_{2n}^{(k)} + T_{2n}^{(k + 1/2)} &= \frac{b - a}{4n}[f(x_k) + f(x_{k + 1/2})] + \frac{b - a}{4n}[f(x_{k + 1/2}) + f(x_{k + 1})] \\ &= \frac{1}{2}T_{n}^{(k)} + \frac{b - a}{2n}f(x_{k + 1/2}) \end{align*}$$

因此：

$$\begin{align*} T_{2n} &= \sum\limits_{k = 0}^{n - 1}\big[T_{2n}^{(k)} + T_{2n}^{(k + 1/2)}\big] \\ &= \frac{1}{2}T_{n} + \frac{b - a}{2n}\sum\limits_{k = 0}^{n - 1}f(x_{k + 1/2}) \end{align*}$$

这说明，我们只需要计算新节点的函数值，即可递推得出区间折半之后的复合梯形公式

复合 Simpson

同理，我们对每一个小区间使用 Simpson 公式，则得到复合 Simpson 公式，最后化简得到的结果为：

$$S_n = \frac{b - a}{6n}\Big[f(a) + 4\sum\limits_{k = 0}^{n - 1}f(\frac{x_k + x_{k + 1}}{2}) + 2\sum\limits_{k = 1}^{n - 1}f(x_k) + f(b)\Big]$$

余项为：

$$R_S = \sum\limits_{k = 0}^{n - 1}R_{S}^{(k)} = -\frac{f^{(4)}(\eta)}{2880n^4}(b - a)^{5} = O(h^{4})$$

复合 Simpson 公式的性质为：

• $S_n$ 是稳定、收敛的
• $S_n$ 具有 3 阶代数精度
• $S_n$ 具有 4 阶准确度

同样，复合 Simpson 也可以递推得到，但是计算复杂度更高

Gauss 求积公式

以上的推导都是基于等距节点的，即 $A_kf(x_k)$ 中的 $x_k$ 是一个已知量，而如果让其成变距节点，则可以得到更高的代数精度

我们定义当 $I_n(f)$ 的代数精度不小于 $2n + 1$ 时，其称为 Gauss 求积公式，积分节点为 Gauss 点

考虑一般的带权积分：

$$I(f) = \int_a^b f(x)\rho(x)dx$$

其对应的 Gauss 求积公式被称为 $\rho(x)$ 对应的 Gauss 积分公式

对于一般的 Gauss 积分公式，$\{x_{k}\}$ 为 Gauss 点的充要条件是：$\omega_{n + 1}(x)$ 与 $\forall P(x) \in \mathbb{P}_{n}$ 带权正交，即：

$$\int_a^b P(x)\omega_{n + 1}(x)\rho(x)dx = 0$$

证明如下：

这引出了如下的求 Gauss 积分公式的方式：

1. 根据权函数，计算 $n + 1$ 次正交多项式的零点，即 Gauss 点
2. 求出带权积分系数：

   $$A_k = \int_a^b l_k(x)\rho(x)dx$$

Gauss 积分公式有如下性质：

• Gauss 积分公式是稳定的、准确的
• 积分余项为：

  $$R_n[f] = \frac{f^{(2n + 2)}(\eta)}{(2n + 2)!}\int_a^b\omega_{n + 1}^{2}(x)\rho(x)dx$$

稳定性的证明方法：

$$A_{k} = \sum\limits_{j = 0}^{n}A_jl_k^2(x_j) = I_n(l_k^2) = \int_a^bl_k^2(x)\rho(x)dx > 0$$

最后一个等号是因为 $l_k(x) \in \mathbb{P}_{2n + 1}$

Gauss Legendre 公式

使用 Legendre 递推式得到正交多项式，利用这组多项式得到 Gauss 点，进而求出 Gauss 求积公式的方法

Gauss-Legendre积分表

上表是使用 Legendre 多项式得到的积分节点与积分系数的表格，由于 Legendre 多项式是定义在 $[-1, 1]$ 上的，因此对于一般区间 $[a, b]$ 上的积分，直接变量代换即可

数值微分

指在未知函数表达式的情况下，利用数值方法计算函数在某一点处的导数值

基本的有限差分公式为：

$$\begin{align*} D_f(x) &= \frac{f(x + h) - f(x)}{h} \\ D_b(x) &= \frac{f(x) - f(x - h)}{h} \\ D_c(x) &= \frac{f(x + h) - f(x - h)}{2h} \end{align*}$$

利用 Taylor 展开可以得到：

$$D_c(x) = f'(x) + \frac{f'''(\xi)}{6}h^2 = f'(x) + O(h^{2})$$

因此 $D_c(x)$ 的准确度为 2 阶

令 $M = \max|f'''(x)|$，则可以得到总误差限为：

$$\varepsilon_{\text{tot}} = \frac{Mh^{2}}{6} + \frac{\varepsilon}{h}$$

其中 $\varepsilon$ 是单次计算 $f$ 的舍入误差，由于分母上有个 2，因此总舍入误差最后是 $2\varepsilon / 2h = \varepsilon / h$

因此，中心差分公式的最小误差为：

$$\varepsilon_{\text{tot}} = \frac{Mh^{2}}{6} + \frac{\varepsilon}{2h} + \frac{\varepsilon}{2h} \geq 3\sqrt[3]{\frac{9M\varepsilon^{2}}{8}}$$

当且仅当 $h = \sqrt[3]{3\varepsilon / M}$ 时取到

同理，我们可以利用差商计算二阶导数，乃至更高阶的导数：

$$G_c(h) \approx 2f[x - h, x, x + h] = \frac{f(x + h) - 2f(x) + f(x - h)}{h^2}$$

同样 Taylor 可以证明这个公式的准确度也是 2 阶

这两种方法实际上都是插值求导的特例，即先根据函数值进行插值，之后对插值得到的函数求导，来近似原函数的导数，即：

$$f^{(i)}(x_j) = \sum\limits_{k = 0}^{n}f(x_k)l_k^{(i)}(x_j)$$

理查德外推法

一种通过将步长减半，组合后得到更准确近似表达式的方法，可以利用在数值积分和微分上，下面以数值微分为例介绍

令：

$$f'(x) - D_c(h) = \alpha_1 h^2 + \alpha_2 h^4 + \cdots + \alpha_k h^{2k} + \cdots$$

其中 $\alpha_k$ 是与 $h$ 无关的常数，则：

$$f'(x) - D_c(\frac{h}{2}) = \frac{\alpha_1}{4} h^2 + \frac{\alpha_2}{16} h^2 + \cdots + \frac{\alpha_k}{2^k} h^{2k} + \cdots$$

通过这两个表达式可以消除 $h^{2}$ 这一项，得到：

$$f'(x) - \frac{4D_c(\frac{h}{2}) - D_c(h)}{3} = O(h^{4})$$

因此，$\dfrac{4D_c(\frac{h}{2}) - D_c(h)}{3}$ 则是一个准确度阶数更高的估计表达式

并且这种方法是可以不断进行下去的，可以递推的不断增加精度