使用泰勒展开的运动预测#

在 \(n\) 维空间运动的物体, 其运动方程总可以写成 \(n\) 个关于时间 \(t\) 的参数方程.

进一步地, 将关于 \(t\) 的参数方程泰勒展开为关于 \(t\) 的多项式函数, 可以显式得到其表达式, 并且可以计算出超过当前时刻, 也就是未来时刻的位置. 当超出当前时刻的时间不长, 运动轨迹是平滑的且分辨率足够高时, 可以认为此级数收敛于实际的位置.

泰勒展开所需的各阶导数可以通过数值计算获得.

参数获取#

假如我们已经获得位置关于时间的一系列数据 \(\{r_n\}\), 其中 \(r_n\) 是当前时刻的数据. 它也可以被记作 \(r_n^{(0)}\). 对应的时刻记为 \(t_n\), 或 \(t_n^{(0)}\).

根据函数导数的公式, 自然地我们想到当前时刻位置的一阶导数为 \(r_n^{(1)} = \frac{r_n - r_{n-1}}{\Delta t_n}\), 其中 \(\Delta t_n := t_n - t_{n-1}\), 以及过去时刻的导数同理.

于是, 我们获得了 \(n-1\) 个一阶导数组成的序列 \(\{r_n^{(1)}\}\). 进一步可以计算二阶乃至 \(n-1\) 阶导数:

另外 \(\Delta t_n^{(0)} := \Delta t_n\).

这样便可递归地求出一阶直到 \(n-1\) 阶导数, 如下方的倒金字塔形图:

1^0   2^0   3^0   4^0   5^0   6^0
   2^1   3^1   4^1   5^1   6^1
      3^2   4^2   5^2   6^2
         4^3   5^3   6^3
            5^4   6^4
               6^5

如此一来, 末位置的 \(0\) 至 \(n-1\) 阶导数已经得到. 据此构建泰勒展开便可预测未来时刻的位置.

需要注意的是, 可以看到 \(r_6^{(5)}\) 实际计算的是整个序列中点时刻的导数, 于是序列持续时间越短, 计算越精确.

增量更新#

我们已经得到末位置 \(r_n\) 的 \(0\) 至 \(n-1\) 阶导数与 \(\Delta t_n^{(0)} \sim \Delta t_n^{(n-1)}\), 若又获取到下一时刻的位置 \(r_{n+1}\) 与时间差 \(\Delta t_{n+1}\), 如何利用现有数据计算 \(r_{n+1}\) 的各阶导数?

注意到 \(r_{n+1}^{(1)}\) 由 \(r_{n+1}^{(0)}\) \(r_{n}^{(0)}\) 计算而来, 而这些数据我们已经知道. 进一步, \(r_{n+1}^{(2)}\) 由 \(r_{n+1}^{(1)}\) \(r_{n}^{(1)}\) 计算得出. 于是我们有 \(r_{n+1}^{(m)} = \frac{r_{n+1}^{(m-1)} - r_{n}^{(m-1)}}{\Delta t_{n+1}^{(m-1)}}\). 如下:

               6^0   7^0
            6^1   7^1
         6^2   7^2
      6^3   7^3
   6^4   7^4
6^5   7^5
   7^6

误差分析#

编程计算时, 因为浮点数所能表示的精度有限, 因此求导多次后误差便会逐步累积变大, 直到不可忽略的程度.

例如, IEEE 754 标准下双精度浮点数 double 表示某一实数, 其误差最大为 \(\epsilon = 2^{-52}\). 这被称为 机器精度.

因此, 若原始数据因为浮点舍入造成的误差为 \(\epsilon\), 计算 \(n\) 阶导数的误差即为 \(\frac{\epsilon}{\Delta t^n}\). 给定能接受的最大误差 \(\delta\), 即有

两边取对数有 \(\ln{\delta} \ge \ln{\epsilon} - n \ln{\Delta t}\), 进而 \(n \ge \frac{\ln{\epsilon / \delta}}{\ln{\Delta t}}\). 若 \(\Delta t > 1\), 则不等号取反.

资源#

根据本文编写的代码可以在此处找到 [Web].