| RMVL 2.4.0 Robotic Manipulation and Vision Library

载入中...
搜索中...
未找到
自动求导、数值微分

包含中心差商以及 Richardson 外推原理的介绍

作者
  赵曦

日期
  2024/05/06

版本
  1.0

上一篇教程：常微分方程（组）数值解与 Runge-Kutta 算法
下一篇教程：一维最优化方法

\[\def\red#1{\color{red}{#1}} \def\teal#1{\color{teal}{#1}} \def\green#1{\color{green}{#1}} \def\transparent#1{\color{transparent}{#1}} \def\orange#1{\color{orange}{#1}} \def\fml#1{\text{(#1)}} \]

## 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 一元函数导数示例代码

- C++

  #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
  }
  rm::derivative
  double derivative(Func1d func, double x, DiffMode mode=DiffMode::Central, double dx=1e-3)
  计算一元函数的导数
  rm::DiffMode::Ridders
  @ Ridders
  Richardson 外推 ridders82
  定义 numcal.hpp:306

- Python

  import rm
  # 自定义函数 f(x) = x²+4x-3
  def quadratic(x):
  return x ** 2 + 4 * x - 3
  # 使用外推原理计算 f'(1)
  dydx = rm.derivative(quadratic, 1, mode=rm.DiffMode.Ridders)
  print("f'(1) = {:.2f}".format(dydx))
  # f'(1) = 6.00

### 3.2 多元函数梯度示例代码

- C++

  #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)
  }
  rm::FuncNd
  std::function< double(const std::valarray< double > &)> FuncNd
  多元函数
  定义 numcal.hpp:299
  rm::grad
  std::valarray< double > grad(FuncNd func, const std::valarray< double > &x, DiffMode mode=DiffMode::Central, double dx=1e-3)
  计算多元函数的梯度

- Python

  import rm
  # 自定义函数 f(x,y) = x²+y²+4x-3y
  func = lambda x: x[0] ** 2 + x[1] ** 2 + 4 * x[0] - 3 * x[1]
  dydx = rm.grad(func, [1, 2])
  print("f'(1, 2) = ({:.2f}, {:.2f})".format(dydx[0], dydx[1]))
  # f'(1, 2) = (6.00, 1.00)