第 6 章：数值积分和数值微分（一）

主要内容

6.1 数值微分
6.2 数值积分
6.3 复化数值积分
6.4 重积分
6.5 高斯型积分

6.1 数值微分

当函数$f(x)$是以离散点列给出时，或当函数的表达式过于复杂时，常用数值微分近似计算$f(x)$的导数$f'(x)$

我们知道，导数是由差商的极限来定义的，有如下三种

\[\begin{align*} f'(x) &= \lim_{h \to 0} \frac{f(x+h) - f(x)}{h} \\ &= \lim_{h \to 0} \frac{f(x) - f(x-h)}{h} \\ &= \lim_{h \to 0} \frac{f(x+h) - f(x-h)}{2h} \end{align*}\]

当极限无法求得时，一种自然的想法是用差商作为极限的近似值，也就是使用函数的差商作为函数导数的近似值
相应的有三种差商形式的数值微分公式，都是用平均的变化率来近似导数

向前差商

使用向前差商近似导数有

\[f'(x_0) \approx \frac{f(x_0 + h) - f(x_0)}{h} \]

$h$是一个增量，称为步长,$x_0+h$的位置在$x_0$的前面，因此称为向前差商。同理有向后、中心差商的定义
由泰勒公式展开

\[f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2!} f''(\xi), \quad \xi \in [x_0, x_0 + h] \]

得到

\[f'(x_0) = \frac{f(x_0 + h) - f(x_0)}{h} - \frac{h}{2} f''(\xi) \]

则向前差商的截断误差阶

\[R(x) = f'(x_0) - \frac{f(x_0 + h) - f(x_0)}{h} = -\frac{h}{2} f''(\xi) = O(h) \]

向后差商

用向后差商近似导数有

\[f'(x_0) \approx \frac{f(x_0) - f(x_0 - h)}{h} \]

与计算向前差商的方法类似，由泰勒展开得到向后差商的截断误差阶

\[R(x) = f'(x_0) - \frac{f(x_0) - f(x_0 - h)}{h} = O(h) \]

中心差商

\[f'(x_0) \approx \frac{f(x_0 + h) - f(x_0 - h)}{2h} = D(h, x_0) \]

由泰勒展开

\[f(x_0 + h) = f(x_0) + h f'(x_0) + \frac{h^2}{2!} f''(x_0) + \frac{h^3}{3!} f'''(\xi_1) \]

\[f(x_0 - h) = f(x_0) - h f'(x_0) + \frac{h^2}{2!} f''(x_0) - \frac{h^3}{3!} f'''(\xi_2) \]

两式相减得

\[\begin{align*} R(x) &= f'(x_0) - \frac{f(x_0 + h) - f(x_0 - h)}{2h} \\ &= -\frac{h^2}{12} \left( f'''(\xi_1) + f'''(\xi_2) \right) \\ &= -\frac{h^2}{6} f'''(\xi) = O(h^2) \\ &\quad x_0 - h \leq \xi \leq x_0 + h \end{align*}\]

中心差商公式可以视为前两种方式的算数平均值，但误差阶由$O(h)$提高到$O(h^2)$，因此也更加常用

差商的几何意义

定义极限$f'(x) = \lim_{h \to 0} \frac{f(x_0 + h) - f(x_0)}{h}$,表示$f(x)$在$x=x_0$处切线的斜率，差商$\frac{f(x_0 + h) - f(x_0)}{h}$表示过$x_0,f(x_0)$和$(x_0+h,f(x_0+h))$两点直线的斜率，是一条通过$x_0$的割线，用割线的斜率代替切线的斜率

设定最佳步长

由上述对截断误差的分析，$h$越小，逼近效果约好，那么实际计算结果是不是随着$h$的减少，精度越高呢？

计算结果见下表（精确值$f'(2)=0.353553$）
由上表可以看到$h=0.1$时，逼近的效果最好，进一步缩小步长，逼近效果反而更差
一般对导数的近似公式进行分析可得到截断误差的表达式，以中心差商为例子，截断误差不超过

\[\frac{h^2}{8} M_3 = \frac{h^2}{6} |f'''(x)| \]

而对于舍入误差，可设$f(x_0-h)$和$f(x_0+h)$分别有舍入误差$\varepsilon_1$和$\varepsilon_2$，令$\varepsilon = \max\{|\varepsilon_1|, |\varepsilon_2|\}$,若忽略$h$的舍入误差，则：

\[\delta(D(h, x_0)) = |D^*(h, x_0) - D(h, x_0)| \leq \frac{|\varepsilon_1| + |\varepsilon_2|}{2h} \leq \frac{\varepsilon}{h} \]

其中$D(h,x_0)$为中心差商公式。即舍入误差可以用$\frac{\varepsilon}{h}$估计
总误差为

\[\frac{h^2}{6} M_3 + \frac{\varepsilon}{h} \]

当$h = \sqrt[3]{\frac{3\varepsilon}{M_3}}$时，误差达到最小，即

\[\overline{M_3} \]

可以看到用误差的表达式确定步长，难度较大，可行性差
实际应用中，通常用事后估计的方法来确定步长$h$
- 基本思想是，如果步长减半，所得的计算结果与原结果的误差在给定的误差范围内，则认为是合适的步长
- 例如，记 $ D(h, x), D\left(\frac{h}{2}, x\right) $ 是步长为 $ h, \frac{h}{2} $ 时由差商得到的 $ f'(x) $ 的近似值，给定误差界 $ \varepsilon $，当 $ \left|D(h, x) - D\left(\frac{h}{2}, x\right)\right| < \varepsilon $ 时，步长 $ h $就是合适的步长。

差商与数值微分

我们先前讨论了用插值函数来近似函数的方法，自然的对于给定的$f(x)$的函数表，建立插值函数$P(x)$后，可以用插值函数$P(x)$的导数近似函数$f(x)$的导数

与差商法相比，插值型求导方法更灵活，也可以用更复杂的插值技术如三次样条插值作为数值微分

6.2 数值积分

对于积分

\[I = \int_{a}^{b} f(x) \, \text{d}x \]

只要找到被积函数$f(x)$的原函数$F(x)$,根据牛顿莱布尼兹公式有

\[I = \int_{a}^{b} f(x) \, \text{d}x = F(b) - F(a) \]

但这种方法求积分经常遇到困难
- 被积的原函数不能用初等函数表示
- 被积函数是以点列的形式给出
- 被积函数的原函数虽然可以用初等函数表示，但过于复杂

建立数值积分公式的思路很多，最常见的有两种

根据积分中值定理，在区间内$[a,b]$存在一点$\xi$，使得

\[\int_{a}^{b} f(x) \, dx = (b - a) f(\xi) \]
- 即底为$(b-a)$而高为$f(\xi)$的矩形的面积恰好等于所求曲边梯形的面积。但$\xi$的位置一般是未知的，因而$f(\xi)$也是未知的，我们将$f(\xi)$称为$f(x)$在区间$[a,b]$上的平均高度。那么，只要对平均高度$f(\xi)$提供一种算法，相应地就获得一种数值求积的方法
比如，用中点公式代替$\xi$，即$\xi = \frac{a+b}{2}$,则可得到中点公式或称为中矩形公式（或简称矩形公式）