1. 基于 Taylor 公式的数值微分

利用以下两个 2 阶 Taylor 公式：

$\begin{align}f(x+h)&=f(x)+f'(x)h+\frac{f''(x)}2h^2+\frac{f'''(x)}{3!}h^3\tag{1-1a}\\ f(x-h)&=f(x)-f'(x)h+\frac{f''(x)}2h^2-\frac{f'''(x)}{3!}h^3\tag{1-1b}\end{align}$

可以得到两个数值微分公式及其对应的截断误差（余项） $R$ ：

$\begin{align}f'(x)&\approx\frac{f(x+h)-f(x)}h,\qquad R=-\frac{f''(\xi)}2h+o(h^2)\tag{1-2a}\\ f'(x)&\approx\frac{f(x)-f(x-h)}h,\qquad R=\frac{f''(\eta)}2h+o(h^2)\tag{1-2b}\end{align}$

其中 $\fml{1-2a}$ 称为， $\fml{1-2b}$ 称为。由于截断误差的存在，差商法的精度不高，因此我们可以考虑使用更高阶的差商法。我们将 $\fml{1-1a}$ 和 $\fml{1-1b}$ 相减，得到：

$f(x+h)-f(x-h)=2f'(x)h+\frac2{3!}f'''(\xi)h^3+o(h^3)\tag{1-3a}$

整理得到

$f'(x)=\frac{f(x+h)-f(x-h)}{2h},\qquad R=-\frac{f'''(\xi)}{3!}h^2+o(h^2)\tag{1-3b}$

如果导数 $f'(x)$ 用 $\fml{1-3b}$ 计算，那么截断误差为 $o(h^2)$ ，比 $\fml{1-2a}$ 和 $\fml{1-2b}$ 的截断误差 $O(h)$ 更小。这种方法称为中心差商。

2. Richardson 外推原理

2.1 公式推导

将公式 $\fml{1-3b}$ 扩展，可以得到更完全的写法：

$f'(x)=\green{\frac{f(x+h)-f(x-h)}{2h}}-\left[\frac1{3!}f^{(3)}(x)h^2+\frac1{5!}f^{(5)}(x)h^4+\frac1{7!}f^{(7)}(x)h^6+\cdots\right]\tag{2-1}$

一般的，可以记作

$L=\green{\varphi(h)}+\left[a_2h^2+a_4h^4+a_6h^6+\cdots\right]\tag{2-2}$

正如上一节所说，如果 $L$ 用 $\varphi(h)$ 来近似表示，则截断误差是按 $h^2$ 展开的幂级数，误差阶为 $o\left(h^2\right)$ 。

事实上，对于公式 $\fml{2-2}$ ，求解 $L$ 的过程还可以继续向前推进（外推），不妨以 $\frac h2$ 替换 $\fml{2-2}$ 中的 $h$ ，可得

$L=\varphi\left(\frac h2\right)+\left[a_2\frac{h^2}4+a_4\frac{h^4}{16}+a_6\frac{h^6}{64}+\cdots\right]\tag{2-3}$

由 $4\times\fml{2-3}-\fml{2-2}$ 可以得到

$L=\left[\frac43\varphi\left(\frac h2\right)-\frac13\varphi(h)\right]-\left[a_4\frac{h^4}4+5a_6\frac{h^6}{16}+\cdots\right]\tag{2-4}$

上式表达了外推过程的第 1 步，说明 $\varphi(h)$ 和 $\varphi\left(\frac h2\right)$ 的一个线性组合提供了导数即 $L$ 的新的计算公式，其至了 $o\left(h^4\right)$

继续进行第 2 步，令

$\psi(h)=\frac43\varphi\left(\frac h2\right)-\frac13\varphi(h)=\frac1{4^1-1}\left[4^1\varphi\left(\frac h2\right)-\varphi(h)\right]\tag{2-5}$

那么 $\fml{2-4}$ 可以改写成

$L=\psi(h)+\left[b_4h^4+b_6h^6+\cdots\right]\tag{2-6}$

与第 1 步的外推过程类似，以 $\frac h2$ 替换 $\fml{2-6}$ 中的 $h$ ，可得

$L=\psi\left(\frac h2\right)+\left[b_4\frac{h^4}{16}+b_6\frac{h^6}{64}+\cdots\right]\tag{2-7}$

由 $4^2\times\fml{2-7}-\fml{2-6}$ 可以得到

$L=\left[\frac{16}{15}\psi\left(\frac h2\right)-\frac1{15}\psi(h)\right]-\left[b_6\frac{h^6}{20}+\cdots\right]\tag{2-8}$

此式表达了外推过程的第 2 步，说明 $\psi(h)$ 和 $\psi\left(\frac h2\right)$ 的一个线性组合提供了导数即 $L$ 的新的计算公式，其至了 $o\left(h^6\right)$

如果我们令

$\theta(h)=\frac{16}{15}\psi\left(\frac h2\right)-\frac1{15}\psi(h)\tag{2-9}$

不难猜出，新的导数计算公式为

$L\approx\frac{64}{63}\theta\left(\frac h2\right)-\frac1{63}\theta(h)\tag{2-10}$

读者可以自行验证，按照这个思路可以进一步外推，可执行任意多步得到精度不断提高的新公式，这就是 Richardson 外推原理。

2.2 外推算法总结

令 $\varphi(h)=\frac{f(x+h)-f(x-h)}{2h}$ ，外推 $M$ 步，则外推算法的步骤如下

① 选取一个适当的 $h$ 值，计算 $M+1$ 个数，记为 $T(*,*)$

$T(n,0)=\varphi\left(\frac h{2^n}\right),\qquad n=0,1,2,\cdots,M\tag{2-11}$

② 按公式计算

$T(n,k)=\frac1{4^k-1}\left[4^kT(n,k-1)-T(n-1,k-1)\right],\qquad\left\{\begin{align} k&=1,2,\cdots,M\\n&=k,k+1,\cdots,M \end{align}\right.\tag{2-12}$

按照以上步骤计算得到的 $T(n,k)$ ，满足等式

$L=T(n,k)+o(h^{2k+2})\qquad(h\to0)\tag{2-13}$

即 $T(n,k)$ 具备 $2k+2$ 阶的精度。

3. 如何使用

RMVL 中提供了一元函数以及多元函数的微分工具，求解一元函数导数时可使用 rm::derivative ，求解多元函数梯度是需要使用 rm::grad 。

以下展示了使用自动求导、数值微分的例子。

3.1 一元函数导数示例代码

#include <cstdio>
#include <rmvl/algorithm/numcal.hpp>
 
// 自定义函数 f(x) = x²+4x-3
inline double quadratic(double x) { return x * x + 4 * x - 3; }
 
int main()
{
    // 使用外推原理计算 f'(1)
    double dydx = rm::derivative(quadratic, 1, rm::DiffMode::Ridders);
    printf("f'(1) = %f\n", dydx);
    //  f'(1) = 6
}

3.2 多元函数梯度示例代码

#include <cstdio>
#include <rmvl/algorithm/numcal.hpp>
 
int main()
{
    // 自定义函数 f(x,y) = x²+y²+4x-3y
    rm::FuncNd func([](const std::valarray<double> &x) { return x[0] * x[0] + x[1] * x[1] + 4 * x[0] - 3 * x[1]; });
 
    // 计算 f'(1, 2)
    auto dydx = rm::grad(func, {1, 2});
    printf("f'(1, 2) = (%f, %f)\n", dydx[0], dydx[1]);
    // f'(1, 2) = (6, 1)
}

目录