 |
RMVL
1.5.0-dev
Robotic Manipulation and Vision Library
|
载入中...
搜索中...
未找到
|
RMVL
1.5.0-dev
Robotic Manipulation and Vision Library
|
包含中心差商以及 Richardson 外推原理的介绍
上一篇教程：常微分方程（组）数值解与 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)}} \]
利用以下两个 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)\)更小。这种方法称为中心差商。
将公式 \(\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 外推原理。
令 \(\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\)阶的精度。
RMVL 中提供了一元函数以及多元函数的微分工具，求解一元函数导数时可使用 rm::derivative ，求解多元函数梯度是需要使用 rm::grad 。
以下展示了使用自动求导、数值微分的例子。
main.cpp CMakeLists.txt 在项目根目录打开终端，输入
可以看到运行结果
1.12.0