| RMVL 1.5.2 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 创建项目

1. 添加源文件 `main.cpp`
   #include <cstdio>
   #include <rmvl/core/numcal.hpp>
   // 自定义函数 f(x)=x²+4x-3
   inline double quadratic(double x) { return x * x + 4 * x - 3; }
   int main()
   {
   double dydx = rm::derivative(quadratic, 1);
   printf("f'(1) = %f\n", dydx);
   }
   rm::derivative
   double derivative(Func1d func, double x, DiffMode mode=DiffMode::Central, double dx=1e-3)
   计算一元函数的导数
   numcal.hpp
   Numerical Calculation Module 数值计算模块
2. 添加 `CMakeLists.txt`
   cmake_minimum_required(VERSION 3.10)
   project(DerivativeDemo)
   find_package(RMVL COMPONENTS core REQUIRED)
   add_executable(demo main.cpp)
   target_link_libraries(demo PRIVATE rmvl_core)

### 3.2 构建、运行

在项目根目录打开终端，输入

mkdir build
cd build
cmake ..
cmake --build .
./demo

可以看到运行结果

f'(1) = 6