4. Rutishauser–Romberg 数值导数收敛三角的计算#
创建时间：2021-02-26
这一份简单笔记会介绍 (Generalized) Rutishauser–Romberg 数值导数收敛三角的 Python 程序的使用与原理。后文会简称 Rutishauser–Romberg 为 RR，Generalized RR 为 GRR。
由于 GRR 方法与 RR 方法仅仅差别在等比数列的比例上，因此不加赘述。我们只讨论 RR 收敛三角。
这篇文档的复现对象是 Medved et al. [1]。
import numpy as np
import findiff
import pandas as pd
np.set_printoptions(precision=4, suppress=True, linewidth=120)
4.1. 问题阐述#
现在我们想要对函数 \(f(x) = \sin(x - 0.5)\) 求其在零点处的三阶导数 \(f^{(3)}(0)\) 的数值结果。
我们假设上述问题有如下限制：
\(f(x)\) 的值可以直接获得，但 \(f\) 的任意阶导数不能直接获得；
\(f(x)\) 只能通过单浮点精度计算机计算 (即存在较大机器精度误差)。
理论上，最理想的结果会是 \(- \cos(- 0.5) \simeq -0.877583\)。但在只能通过数值导数计算的情况下，这个值并不是轻易可以得到并且信任的。
这样模型问题看起来比较蠢。但现实中，在遇到需要数值地求取高阶导数，却无法保证函数本身精度的情况下，这种分析模式就非常关键。例子是量子化学中的高阶极化率数值计算问题。
def f(x):
return np.array(np.array(np.sin(x-0.5), dtype=np.float32), dtype=float)
4.1.1. 问题的困难：过低或过高的间隔的数值误差#
对于三阶导数问题，采用五点差分时，
原则上，数值导数的间隔 \(h\) (interval) 越小越好。如果 \(h\) 过大 (譬如取 \(h = 0.5\))，那么结果就会与理论值 (-0.877583) 相差很大：
x, h = 0, 0.5
1/h**3 * (-0.5*f(x-2*h) + f(x-h) - f(x+h) + 0.5*f(x+2*h))
-0.8240854740142822
但如果间隔 \(h\) 太小 (譬如取 \(h = 0.005\))，那么又会出于机器精度，反而会有更大的偏差：
x, h = 0, 0.005
1/h**3 * (-0.5*f(x-2*h) + f(x-h) - f(x+h) + 0.5*f(x+2*h))
-0.9536743164062499
因此，间隔 \(h\) 必须要适中 (譬如取 \(h = 0.05\))，才有可能得到正确的结果：
x, h = 0, 0.05
1/h**3 * (-0.5*f(x-2*h) + f(x-h) - f(x+h) + 0.5*f(x+2*h))
-0.877261161804199
不仅如此，如果采用下述七点差分公式时，精度或许 (也或许不) 会有些许提升：
差分点数量越多，一般结果越精确。但这同时也受限于方才提到的大间隔或机器精度的问题。
1/h**3 * 1/48 * (f(x-4*h) - 34*f(x-2*h) + 64*f(x-h) - 64*f(x+h) + 34*f(x+2*h) - f(x+4*h))
-0.8778721094131468
因此，随即而来的问题是，
间隔 \(h\) 在什么情况下可以被信任，可能得到正确的结果？
采用几点差分会有更好的效果？
4.1.2. 问题的解决：RR 收敛三角#
通过上面简单的分析，我们知道合理的间隔在 \(0.005\) 到 \(0.5\) 之间。因此，我们设计以两倍为间隔的等比数列，称为偏移值数列 (英文一般称为 offsets) offsets_half：
offsets_half = np.array([0.004 * 2**n for n in range(10)])
offsets_half
array([0.004, 0.008, 0.016, 0.032, 0.064, 0.128, 0.256, 0.512, 1.024, 2.048])
我们对这些数，求取 \(f(x)\)，生成的数组为 fx_pos \([f(0.004), f(0.008), \cdots, f(2.048)]\)：
fx_pos = f(offsets_half)
fx_pos
array([-0.4759, -0.4724, -0.4653, -0.4511, -0.4223, -0.3635, -0.2416, 0.012 , 0.5003, 0.9997])
我们同时求取 \(f(-x)\)，生成的数组为 fx_neg \([f(-0.004), f(-0.008), \cdots, f(-2.048)]\)：
fx_neg = f(-offsets_half)
fx_neg
array([-0.4829, -0.4864, -0.4934, -0.5073, -0.5346, -0.5875, -0.686 , -0.8479, -0.9989, -0.5593])
最后需要求取 \(f(0)\) f0：
f0 = f(0)
f0
array(-0.4794)
通过下述程序给出并绘制 RR 收敛三角。其中 grrtrig 程序可供下载 grrtrig.py，而大多数函数也会在本文档重新说明。
import grrtrig
grr = grrtrig.calculate_GRR_trig(offsets_half, 3, fx_pos, fx_neg, f0)
df, _ = grrtrig.output_pd_grr_trig(offsets_half, grr)
df
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|---|
| 0.004 | -0.931323 | -0.941024 | -0.943126 | -0.943630 | -0.943755 | -0.943786 | -0.943793 | -0.943795 |
| 0.008 | -0.902219 | -0.909495 | -0.911364 | -0.911835 | -0.911953 | -0.911982 | -0.911989 | |
| 0.016 | -0.880391 | -0.881452 | -0.881722 | -0.881791 | -0.881808 | -0.881813 | ||
| 0.032 | -0.877208 | -0.877397 | -0.877388 | -0.877386 | -0.877386 | |||
| 0.064 | -0.876639 | -0.877527 | -0.877527 | -0.877527 | ||||
| 0.128 | -0.873975 | -0.877533 | -0.877554 | |||||
| 0.256 | -0.863299 | -0.877214 | ||||||
| 0.512 | -0.821555 |
绿色的值中，存有精确的三阶数值导数 \(f^{(3)}(0)\) 的概率比较大。
4.2. 实现原理#
从这里开始仅仅是程序笔记。如果只关心如何使用 RR 收敛三角，上面的文档与程序应当已经足够了。
4.2.1. 给定微小偏移函数值下的任意阶数值导数#
首先是考察下述问题。在给定任意偏移列表 (offsets) \(\{ a_m \}\) 或写作向量 \(\boldsymbol{a}\) 的情况下，其任意阶数值导数 \(f^{(n)}\) 应该如何给出？
数值导数的形式会是 (其中 \(\{ b_m \}\) 为待求量，它与导数阶数 \(n\) 有关，我们称 \(\{ b_m \}\) 为数值差分系数 Finite Difference Coefficients)
举例来说，五点差分的三阶导数情况下，\(m = 5, n = 3\)，并且 \(\{ a_m \} = \{ -2, -1, 0, 1, 2 \}\)。那么由此得到的 \(\{ b_m \} = \{ - 1/2, 1, 0, -1, 1/2 \}\)。但这是在间隔 \(h = 1\) 的情形下给出的。如果是其他间隔，那么
这类数值差分系数可以很容易地从 Wikipedia: Finite difference coefficient 上获得。
但我们实际会遇到的是更复杂的差分情形。譬如方才提到的三阶导数的七点差分。一般的七点差分是从 -3 到 3 的等差数列；但我们方才用到的却是从 -4 到 4 的正负两条等比数列。这要求我们能对任意情形的偏移列表 \(\{ a_m \}\) 进行系数 \(\{ b_m \}\) 的计算。
其计算原理这里不详细展开。详情参考 Crtaylor 的交互网页。实现过程非常简单。
\(m \times m\) 矩阵 \(\mathbf{M}\)，其中矩阵元 \(M_{ij} = a_j^i\)。
\(m\) 维度向量 \(\boldsymbol{r}\)，其中 \(r_n = n!\)，其余值均为零。
随后求解 \(\mathbf{M} \boldsymbol{b} = \boldsymbol{r}\) 即可。
def calculate_findiff_coefs(offsets, deriv):
"""
Calculate Finite Difference Coefficients
Parameters
----------
offsets : array_like
deriv : int
Returns
-------
coefs : ndarray
Examples
--------
>>> calculate_findiff_coefs([-2, -1, 0, 1, 2], 4)
array([ 1., -4., 6., -4., 1.])
>>> calculate_findiff_coefs([-2, -1, 0, 1, 2, 1.99], 3)
array([ -0.1867, -0.6722, 3.7688, -6.0505, -124.5 , 127.6406])
References
----------
https://web.media.mit.edu/~crtaylor/calculator.html
https://github.com/maroba/findiff/blob/e8ca33707e3e25d76bf0f93b2391e466209287b1/findiff/coefs.py
"""
if len(offsets) < deriv + 1:
raise ValueError("Length of offsets should be larger than derivative order plus 1.")
if len(offsets) != len(set(offsets)):
# Note that this program could not handle pathological cases, only make simple check instead.
raise ValueError("Possibly exactly same offset value is given. Please check `offsets'.")
offsets = np.asarray(offsets)
matrix = np.array([offsets**n for n in range(len(offsets))])
rhs = np.zeros(len(offsets))
rhs[deriv] = np.math.factorial(deriv)
return np.linalg.solve(matrix, rhs)
4.2.2. GRR 收敛三角具体实现#
这里根据 Medveď 文章式 (15) 来实现。GRR 收敛三角的首列通过上述程序进行数值差分。如果我们令 GRR 收敛三角为矩阵 \(\mathbf{P}\)，其矩阵元用 \(P_{r, c}\) 表示，那么 \(c >= 1\) 时，
以此方式迭代即可。
事实上，如果收敛三角任意列可以全部使用上述 calculate_findiff_coefs 函数求得，并且与迭代表达式等价。但出于数值稳定性考虑，仍然选择使用迭代式。
def calculate_GRR_trig(offsets_half, deriv, fx_pos, fx_neg, f0=None):
"""
Calculate (Generalized) Rutishauser–Romberg Triangle
Parameters
----------
offsets_half : array_like
Positive half part of offsets. For example, if all offsets are [-2, -1, 0, 1, 2],
then `offsets_half' should be [1, 2].
Must be a geometric sequence. This function does not make double check on this.
Do not contain zero in this array.
deriv : int
fx_pos : array_like
Value list of f(offsets_half). Dimension should be the same to `offsets_half'.
fx_neg : array_like
Value list of f(-offsets_half). Dimension should be the same to `offsets_half'.
f0 : float or None
Value of f(0). May leave as None if derivative order is odd number.
Returns
-------
grr_trig : ndarray
(Generalized) Rutishauser–Romberg triangle.
"""
if deriv % 2 == 0 and f0 is None:
raise ValueError("`f0' should be provided if derivative order is even number.")
else:
f0 = 0 if f0 is None else f0
if len(offsets_half) < 2:
raise ValueError("length of `offsets_half' must >= 2 in order to calculate ratio.")
comp_len = (deriv + 1) // 2
mat_size = len(offsets_half) - comp_len
grr_trig = np.zeros((mat_size, mat_size))
ratio = offsets_half[-1] / offsets_half[-2]
for r in range(mat_size):
i_end = r + comp_len
offsets = np.concatenate([offsets_half[r:i_end], -offsets_half[r:i_end], [0]])
coef_list = calculate_findiff_coefs(offsets, deriv)
val_list = np.concatenate([fx_pos[r:i_end], fx_neg[r:i_end], [f0]])
grr_trig[r, 0] = (coef_list * val_list).sum()
for c in range(1, mat_size):
for r in range(mat_size-c):
grr_trig[r, c] = (ratio**(2*c) * grr_trig[r, c-1] - grr_trig[r+1, c-1]) / (ratio**(2*c) - 1)
for r in range(mat_size):
for c in range(mat_size-r, mat_size):
grr_trig[r, c] = np.nan
return grr_trig
对于文档开头 \(f(x) = \sin(x - 0.5)\) 的问题，其 RR 收敛三角为
calculate_GRR_trig(offsets_half, 3, fx_pos, fx_neg, f0)
array([[-0.9313, -0.941 , -0.9431, -0.9436, -0.9438, -0.9438, -0.9438, -0.9438],
[-0.9022, -0.9095, -0.9114, -0.9118, -0.912 , -0.912 , -0.912 , nan],
[-0.8804, -0.8815, -0.8817, -0.8818, -0.8818, -0.8818, nan, nan],
[-0.8772, -0.8774, -0.8774, -0.8774, -0.8774, nan, nan, nan],
[-0.8766, -0.8775, -0.8775, -0.8775, nan, nan, nan, nan],
[-0.874 , -0.8775, -0.8776, nan, nan, nan, nan, nan],
[-0.8633, -0.8772, nan, nan, nan, nan, nan, nan],
[-0.8216, nan, nan, nan, nan, nan, nan, nan]])
4.3. 实现补充#
4.3.1. GRR 收敛三角的合理收敛值确定#
如果不考虑机器精度所导致的误差，那么收敛三角的最右上方一定是精度 (accuracy) 阶数最高的结果。但很显然，事实是这个值反而是偏差最大的点。因此，我们需要问，上述三角的那个结果是可以信任的？
一般来说，与周围数值偏差最小的点是可以信任的。上述三角的第 1-4 列、第 3-6 行看起来比较可靠。但我们需要一种方法较为系统地评价值是否可以信任。
这里的评价标准是与下方和左方的数值进行比较；如果误差较小，那么就可以信任。因此，首列与对角线没有被纳入评判。事实上，对角线值往往更准确，因此最终评价的时候还是需要人看看周围的数值是否合理。
def check_grr_trig_converge(grr_trig):
"""
Convergence Check Matrix of (Generalized) Rutishauser–Romberg Triangle
Parameters
----------
grr_trig : ndarray
(Generalized) Rutishauser–Romberg triangle.
Return
------
mat_chk : ndarray
"""
n = len(grr_trig)
mat_chk = np.zeros((n, n))
for r in range(0, n-1):
for c in range(1, n-r-1):
mat_chk[r, c] = np.abs(grr_trig[r, c] - grr_trig[r+1, c]) + np.abs(grr_trig[r, c] - grr_trig[r, c-1])
for r in range(n):
mat_chk[r, 0] = np.nan
for r in range(0, n):
for c in range(n-r-1, n):
mat_chk[r, c] = np.nan
return mat_chk
grr = calculate_GRR_trig(offsets_half, 3, fx_pos, fx_neg, f0)
check_grr_trig_converge(grr)
array([[ nan, 0.0412, 0.0339, 0.0323, 0.0319, 0.0318, 0.0318, nan],
[ nan, 0.0353, 0.0315, 0.0305, 0.0303, 0.0302, nan, nan],
[ nan, 0.0051, 0.0046, 0.0045, 0.0044, nan, nan, nan],
[ nan, 0.0003, 0.0001, 0.0001, nan, nan, nan, nan],
[ nan, 0.0009, 0. , nan, nan, nan, nan, nan],
[ nan, 0.0039, nan, nan, nan, nan, nan, nan],
[ nan, nan, nan, nan, nan, nan, nan, nan],
[ nan, nan, nan, nan, nan, nan, nan, nan]])
4.3.2. GRR 收敛三角在 Pandas 的绘制#
上述结果可以通过图像表格呈现。浅绿色打底的单元格是上面给出的误差最小的若干单元格附近的格点；其数值导数可以被相信的概率较大。这些单元格外围一圈也有值得考虑的价值。
def output_pd_grr_trig(offsets_half, grr_trig, tolerance=3):
"""
Pandas Presentation of (Generalized) Rutishauser–Romberg Triangle
Parameters
----------
offsets_half : array_like
grr_trig : ndarray
(Generalized) Rutishauser–Romberg triangle.
tolerance : int
Number of minimum difference cells in convergence check matrix.
Return
------
df : pandas.io.formats.style.Styler
Pandas show of GRR triangle.
df_check : pandas.io.formats.style.Styler
Pandas show of convergence check matrix of GRR triangle.
"""
n = len(grr_trig)
df = pd.DataFrame(grr_trig, columns=range(n), index=offsets_half[:n])
df_check = pd.DataFrame(check_grr_trig_converge(grr_trig), columns=range(n), index=offsets_half[:n])
df.replace(np.nan, "", regex=True, inplace=True)
df_check.replace(np.nan, "", regex=True, inplace=True)
t = check_grr_trig_converge(grr_trig).flatten()
t = t.argsort()[:3]
t = np.array([t // n, t % n]).T
def highlight_cells(x):
df = x.copy()
df.loc[:,:] = ''
for r, c in t:
df.iloc[r, c] = "background-color: lightgreen"
df.iloc[r, c-1] = "background-color: lightgreen"
df.iloc[r+1, c] = "background-color: lightgreen"
return df
df = df.style.apply(highlight_cells, axis=None)
df_check = df_check.style.apply(highlight_cells, axis=None)
return df, df_check
df, df_check = output_pd_grr_trig(offsets_half, grr)
df
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|---|
| 0.004 | -0.931323 | -0.941024 | -0.943126 | -0.943630 | -0.943755 | -0.943786 | -0.943793 | -0.943795 |
| 0.008 | -0.902219 | -0.909495 | -0.911364 | -0.911835 | -0.911953 | -0.911982 | -0.911989 | |
| 0.016 | -0.880391 | -0.881452 | -0.881722 | -0.881791 | -0.881808 | -0.881813 | ||
| 0.032 | -0.877208 | -0.877397 | -0.877388 | -0.877386 | -0.877386 | |||
| 0.064 | -0.876639 | -0.877527 | -0.877527 | -0.877527 | ||||
| 0.128 | -0.873975 | -0.877533 | -0.877554 | |||||
| 0.256 | -0.863299 | -0.877214 | ||||||
| 0.512 | -0.821555 |