2. 一元函数数值导数简要介绍

在这一节中,我们将会简单地介绍在我们工作中可能会使用到的数值导数方法。

在这里我们不引入复杂的化学问题,几乎单纯地从非常基础的数学角度来阐释问题。但作为仅仅关心 XYG3 二阶梯度的文档,这一节将会是后面几节的基础与铺垫。在这份课题中,数值导数的意义是可以用来验证解析导数的正确性。

[1]:
%matplotlib notebook

import numpy as np
from functools import partial
from pyscf import gto, scf
import matplotlib.pyplot as plt

np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.set_printoptions(5, linewidth=150, suppress=True)

我们在这一节中,只讨论一元函数的导数问题。关于更高维度的向量数值导数、我们会结合具体的分子来进行说明。

2.1. 用于讨论的 Morse 函数

我们现在拿一个化学中经常使用的 Morse 势能函数来举例子。其中的参数仅仅是用于演示用途,而没有实际的物理意义。

\[f(x) = \left( 1 - e^{1 - x} \right)^2 - 1\]
[2]:
def f(x):
    return (1 - np.exp(1 - x))**2 - 1

该函数的精确 (解析) 一阶导数可以表示为 f_p(x) \(f'(x)\)

\[f'(x) = \frac{\partial f(x)}{\partial x} = 2 e^{1 - x} \left( 1 - e^{1 - x} \right)\]
[3]:
def f_p(x):
    return 2 * np.exp(1 - x) * (1 - np.exp(1 - x))

该函数 \(f(x)\)\(x = 1\) 处取到极小值 (即 \(f'(1) = 0\))。函数显示为蓝色曲线,一阶导数显示为橙色曲线:

[4]:
xlist = np.arange(0, 5, 0.01)
ylist, yplist = f(xlist), f_p(xlist)

fig, ax = plt.subplots(figsize=(4.5, 3))
ax.plot([0, 5], [0, 0], color="black")
ax.plot(xlist, ylist, label="$f(x)$")
ax.plot(xlist, yplist, label="$f'(x)$")
ax.set_xlim(0, 5), ax.set_ylim(-1, 1)
ax.set_title(r"Plot of $f(x) = x^2$")
ax.set_xlabel(r"$x$")
ax.legend(), ax.grid()
fig.tight_layout()

2.2. 三点差分一阶数值导数

我们首先从最简单的情况开始讨论。对于数值导数而言,最简单的情形是三点差分。

对于一阶导数,如果我们希望得到 \(x = 3\) 的导数即 \(f'(3)\),并且已经知道了整个曲线的形状,我们可以用下面图示的方法:

[5]:
xlist = np.arange(2, 4, 0.01)
ylist = f(xlist)
ytlist = f_p(3) * (xlist - 3) + f(3)

fig, ax = plt.subplots(figsize=(4.5, 3))
ax.plot(xlist, ylist, label="$f(x)$")
ax.plot(xlist, ytlist, label="tangent of $(3, f(3))$", linestyle=":")
ax.scatter([2, 3, 4], [f(2), f(3), f(4)], c="C2")
ax.plot([2, 4], [f(2), f(4)], c="C2", label="segment when $h = 1$")
ax.set_xlim(1.8, 4.2), ax.set_ylim(-0.7,-0.0)
ax.set_title(r"Plot of tangent of $(3, f(3))$")
ax.set_xlabel(r"$x$")
ax.legend(), ax.grid()
fig.tight_layout()

上面绿色的三个点与线段,可以用来估计 \(f(x)\)\(x = 3\) 处的导数 (点 \((3, f(3))\) 处曲线切线的斜率)。我们定义函数 f_d3(x, h) \(f^\mathrm{d} (x, h)\)

\[f^\mathrm{d} (x, h) = \frac{f(x + h) - f(x - h)}{2h}\]
[6]:
def f_d3(x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

我们说,该函数的意义是取绿色线段的斜率。对于上图而言,\(h\) 的取值是 \(1\),那么绿色线段左右两个点就是 \((1 - h, f(1 - h)) = (2, f(2))\)\((1 + h, f(1 + h)) = (4, f(4))\)。这两个线段端点的斜率即

\[f^\mathrm{d} (3, 1) = \frac{f(4) - f(2)}{4 - 2} \simeq 0.2517\]
[7]:
f_d3(3, 1)
[7]:
0.2516641072736052

我们可以估计地认为,函数 \(f(x)\)\(x = 3\) 处切线的斜率 (数值导数) 接近于 \(0.2517\)。从肉眼上,我们也确实会觉得绿色的线段 (估计的导数) 与橙色的切线 (解析的导数) 是非常接近的。实际上,解析的导数是 \(f'(3) \simeq 0.2340\)

[8]:
f_p(3)
[8]:
0.23403928869575705

我们会说,若 \(h\) 的值越小,或者说绿色的点之间的距离越小,那么我们所求得的斜率 (数值导数) 会越精确。事实上,

\[\lim_{h \rightarrow 0} f^\mathrm{d} (x, h) = f'(x)\]

并且数学分析或数值分析可以告诉我们,\(|f^\mathrm{d} (x, h) - f'(x)| = o(h)\),因此函数 \(f^\mathrm{d}(x, h)\)\(f'(x)\) 的逼近精度随 \(h\) 越小而线性地更精确。

但事实上,从程序的结果而言并非如此。我们可以对 \(x = 3\) 时,\(|f^\mathrm{d} (3, h) - f'(3)|\) 的值作图,发现当 \(h = 10^{-5}\) 左右时,才能获得比较好的数值导数的精度 (\(10^{-10}\) 左右的精度)。因此,一般而言,用于求数值导数的差分值 \(h\)

  • \(h\) 较大的情况,数值导数的精度受制于 \(|f^\mathrm{d} (x, h) - f'(x)|\)
  • \(h\) 较小的情况,数值导数的精度受制于计算机本身的精度 (机器精度);若计算机可以提供相当高的机器精度,那么 \(h\) 较小时也可以得到良好的数值精度。

但我们一般在运算时只使用双浮点数,因此我们我们也不能选择太过小的 \(h\) 值来计算数值导数。选取一个合适的 \(h\) 值会是之后始终会遇到的问题。

[9]:
hlist = np.logspace(-17, 1, 1000)
ylist = np.abs(f_d3(3, hlist) - f_p(3))

fig, ax = plt.subplots(figsize=(4.5, 3))
ax.plot(hlist, ylist)
ax.set_xscale("log"), ax.set_yscale("log")
ax.set_xlabel(r"$h$")
ax.set_title("Log plot of $|f^\mathrm{d} (3, h) - f'(3)|$")
fig.tight_layout()

我们指出,尽管三点差分声称是“三点”,但一阶梯度的计算只使用了两个点。但它被称为三点差分,是由于其数学推导而致使的。

2.3. 五点差分一阶数值导数

五点差分并不像三点差分一样,有直观的几何解释。我们仍然用 \(f^\mathrm{d} (x, h)\) 记号表示数值导数,但程序的函数写为 f_d5(x, h)

\[f^\mathrm{d} (x, h) = \frac{f(x - 2h) - 8 f(x - h) + 8 f(x + h) - f(x + 2h)}{12 h}\]
[10]:
def f_d5(x, h):
    return (f(x - 2 * h) - 8 * f(x - h) + 8 * f(x + h) - f(x + 2 * h)) / (12 * h)

我们仍然先拿 \(x = 3, h = 1\) 的情况来讨论。我们看到,\(f'(3) \simeq 0.2340\),而计算代价低的三点差分大约是 \(0.2517\),代价较高的五点差分大约是 \(0.2552\)

[11]:
f_p(3), f_d3(3, 1), f_d5(3, 1)
[11]:
(0.23403928869575705, 0.2516641072736052, 0.25524346096060413)

可见代价更高的五点差分并不一定会给出更好的结果。但是五点差分有可能达到更低的逼近精度:在 \(x = 3, h = 10^{-3}\) 时,五点差分的数值精度可以达到 \(10^{-13}\) 量级。

[12]:
f_d3(3, 1e-3) - f_p(3), f_d5(3, 1e-3) - f_p(3)
[12]:
(2.0690889585006644e-08, -1.3128387266192476e-14)

我们可以绘制 \(x = 3\) 时,不同的逼近参数 \(h\) 下,三点差分与五点差分的比较:

[13]:
hlist = np.logspace(-17, 1, 1000)
y3list = np.abs(f_d3(3, hlist) - f_p(3))
y5list = np.abs(f_d5(3, hlist) - f_p(3))

fig, ax = plt.subplots(figsize=(4.5, 3))
ax.plot(hlist, y3list, label="3-point stencil")
ax.plot(hlist, y5list, label="5-point stencil")
ax.set_xscale("log"), ax.set_yscale("log")
ax.set_xlabel(r"$h$")
ax.set_title("Log plot of $|f^\mathrm{d} (3, h) - f'(3)|$")
ax.legend()
fig.tight_layout()

一般而言,

  • 五点差分的计算代价相对三点差分来说更大。对于一元函数的一阶导数而言,五点差分的代价是三点差分的两倍;但随着函数的未知元数量增大,其计算量还会进一步增大。
  • 五点差分在逼近参数 \(h\) 较大时,未必有更好的逼近表现。
  • 五点差分从数学的角度来讲,\(|f^\mathrm{d} (x, h) - f'(x)| = o(h^2)\);因此函数 \(f^\mathrm{d}(x, h)\)\(f'(x)\) 的逼近精度随 \(h\) 越小而平方地更精确 (三点差分是线性地更精确)。因此,若逼近精度不受制于机器精度,那么在 \(h\) 较小的情况下,五点差分可以更快地逼近到解析导数。
  • 但五点差分也更容易受制于机器精度影响,因此 \(h\) 太小时,误差也未必比三点差分更小。在当前的例子中,\(h\) 较小的情况下,五点与三点差分的误差近似相等;但对于以后实际的例子,太小的 \(h\) 很可能带来更为灾难性的结果。
  • 希望用较大的 \(h\) 来逼近时,可以选用五点差分;用较小的 \(h\) 逼近时,选用三点差分。
  • 尽管称为五点差分,但实际上只使用了四个点。这也是由于五点差分的推导所导出的结果。

在以后,我们通常只使用三点查分来解决问题。

2.4. 三点差分二阶数值导数

现在我我们考虑二阶导数的计算。我们首先指出,二阶导数 f_pp(x) \(f''(x)\)

\[f''(x) = 4 e^{2 - 2x} - 2 e^{1 - x}\]
[14]:
def f_pp(x):
    return 4 * np.exp(2 - 2 * x) - 2 * np.exp(1 - x)

二阶数值导数的一种做法,可以是通过一阶数值导数 \(f^\mathrm{d} (x, h)\) 再求一次数值导数获得 (连续两次对函数 \(f(x)\) 的数值导数),其公式可以表达为 f_dd(x, h) \(f^\mathrm{dd} (x, h)\)

\[f^\mathrm{dd} (x, h) = \frac{f(x - h) - 2 f(x) + f(x + h)}{h^2}\]
[15]:
def f_dd(x, h):
    return (f(x - h) - 2 * f(x) + f(x + h)) / h**2

但还有一种做法:如果我们已知一阶解析导数,那么可以通过对解析导数再求一次数值导数获得 (仅一次数值导数),其公式可以表达为 f_pd(x, h) \(f'{}^\mathrm{d} (x, h)\)

\[f'{}^\mathrm{d} (x, h) = \frac{f'(x + h) - f'(x - h)}{2h}\]
[16]:
def f_pd(x, h):
    return (f_p(x + h) - f_p(x - h)) / (2 * h)

我们仍然可以绘制数值导数与解析导数的在 \(x = 3\) 时,不同的逼近参数 \(h\) 的误差图:

[17]:
hlist = np.logspace(-9, 1, 1000)
dlist = np.abs(f_d3(3, hlist) - f_p(3))
ddlist = np.abs(f_dd(3, hlist) - f_pp(3))
pdlist = np.abs(f_pd(3, hlist) - f_pp(3))

fig, ax = plt.subplots(figsize=(4.5, 3))
ax.plot(hlist, dlist, label=r"$f^\mathrm{d} (x, h)$")
ax.plot(hlist, ddlist, label=r"$f^\mathrm{dd} (x, h)$")
ax.plot(hlist, pdlist, label=r"$f'{}^\mathrm{d} (x, h)$")
ax.set_xscale("log"), ax.set_yscale("log")
ax.set_xlabel(r"$h$")
ax.set_title("Error of 3-Point Stencil Derivatives")
ax.legend()
fig.tight_layout()

我们在上图展示了三个数值导数的误差:

  • \(f(x)\) 的两次数值导数 (橙色曲线) 与对 \(f'(x)\) 的一次数值导数 (绿色曲线) 在稍大的 \(h\) 值 (\(h > 10^{-3}\)) 时,误差相近;但两次数值导数 (橙色) 应机器精度导致的误差来得更快;
  • 一般来说,我们很难判断不同函数的数值导数的误差大小 (拿当前的例子而言,我们不能判断对 \(f(x)\) 还是对 \(f'(x)\) 的数值导数误差是哪个更大),但由于机器精度导致的数值导数振荡的 \(h\) 值都在 \(10^{-4.5}\) 左右。

基于上面的原因,我们以后不太会使用两次数值导数,并且通常 \(h\) 的值设在 \(3 \times 10^{-4}\)