导数的有限差分逼近¶
简介¶
导数的有限差分逼近在数值分析和计算物理学中非常重要。在本教程中，我们将展示如何使用 SymPy 计算不同精度的逼近。希望这些笔记对那些正在使用某种语言开发代码且需要能够有效地生成各种逼近的有限差分公式的实际研究人员有所帮助。
为了建立符号，我们首先假设存在一个单变量 x 的连续函数 F，F 具有任意阶导数。我们在实数轴上以 h 为间隔均匀地采样 x 值。在大多数情况下，我们希望 h 在某种意义上很小。F(x) 可以通过通常的泰勒级数展开式在某个点 \(x_{0}\) 处展开。设 \(x = x_{0} + h\)。那么泰勒展开式为
为了简化符号，我们现在定义一组系数 \(c_{n}\)，其中
所以现在我们的级数形式为
在下面，我们只使用一组有限的 \(x_{i}\) 值，其中 \(i\) 从 \(1,...,N\) 开始，并且我们函数 F 在这些网格点上的对应值用 \(F_{i}\) 表示。所以问题是如何在使用大小为 N 的有限集合 \((x_{i},F_{i})\) 的子集的约束下生成 F 的导数的近似值。
以下是使用 SymPy 对给定阶导数进行近似和评估其精度的操作。首先，我们使用 SymPy 通过一种相当蛮力的方法来推导出近似值，这种方法通常在入门教程中介绍。稍后我们将使用其他 SymPy 函数，它们可以更有效地完成工作。
使用 SymPy 矩阵的直接方法¶
如果我们令 \(x_{0} = x_{i}\)，在 \(x_{i+1}=x_{i}+ h\) 处评估级数并截断所有高于 \(O(h^1)\) 的项，我们可以求解单个系数 \(c_{1}\) 并得到一阶导数的近似值
其中 \(O(h)\) 指的是 \(h\) 中级数的最低阶项。这表明导数近似值是一阶精度的。换句话说，如果我们决定只能使用两个对 \((x_{i},F_{i})\) 和 \((x_{i+1},F_{i+1})\)，我们得到一个“一阶精度”的导数。
除了 \((x_{i},F_{i})\) 之外，我们接下来使用两个点 \((x_{i+1},F_{i+1})\) 和 \((x_{i+2},F_{i+2})\)。然后我们有两个方程
如果我们再次想要找到一阶导数 (\(c_{1}\))，我们可以通过消除两个方程中包含 \(c_{2}\) 的项来做到这一点。我们展示了如何使用 SymPy 来做到这一点。
>>> from __future__ import print_function
>>> from sympy import *
>>> x, x0, h = symbols('x, x_0, h')
>>> Fi, Fip1, Fip2 = symbols('F_{i}, F_{i+1}, F_{i+2}')
>>> n = 3 # there are the coefficients c_0=Fi, c_1=dF/dx, c_2=d**2F/dx**2
>>> c = symbols('c:3')
>>> def P(x, x0, c, n):
... return sum( ((1/factorial(i))*c[i] * (x-x0)**i for i in range(n)) )
右侧向量
>>> R = Matrix([[Fi], [Fip1], [Fip2]])
现在我们构建一个矩阵，该矩阵包含 nth 次多项式 P 中 c_i 的系数。
在 \(x_i\) 处评估的 \(c_i\) 的系数
>>> m11 = P(x0 , x0, c, n).diff(c[0])
>>> m12 = P(x0 , x0, c, n).diff(c[1])
>>> m13 = P(x0 , x0, c, n).diff(c[2])
在 \(x_i + h\) 处评估的 \(c_i\) 的系数
>>> m21 = P(x0+h, x0, c, n).diff(c[0])
>>> m22 = P(x0+h, x0, c, n).diff(c[1])
>>> m23 = P(x0+h, x0, c, n).diff(c[2])
在 \(x_i + 2*h\) 处评估的 \(c_i\) 的系数
>>> m31 = P(x0+2*h, x0, c, n).diff(c[0])
>>> m32 = P(x0+2*h, x0, c, n).diff(c[1])
>>> m33 = P(x0+2*h, x0, c, n).diff(c[2])
系数矩阵在这种情况下是 3x3 的
>>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
三个方程关于 \(c_i\) 的矩阵形式为 M*X = R
解是通过直接求解 3x3 矩阵 M 的逆得到的
>>> X = M.inv() * R
请注意，所有三个系数都构成了解。所需的导数是系数 \(c_1\)，即 X[1]。
>>> print(together(X[1]))
(4*F_{i+1} - F_{i+2} - 3*F_{i})/(2*h)
计算另一个三点近似于一阶导数很有启发性，除了将近似值集中在 \(x_i\) 处，因此使用 \(x_{i-1}\)、\(x_{i}\) 和 \(x_{i+1}\) 处的点。所以以下是使用“蛮力”方法完成此操作的方法
>>> from __future__ import print_function
>>> from sympy import *
>>> x, x0, h = symbols('x, x_i, h')
>>> Fi, Fim1, Fip1 = symbols('F_{i}, F_{i-1}, F_{i+1}')
>>> n = 3 # there are the coefficients c_0=Fi, c_1=dF/h, c_2=d**2F/h**2
>>> c = symbols('c:3')
>>> # define a polynomial of degree n
>>> def P(x, x0, c, n):
... return sum( ((1/factorial(i))*c[i] * (x-x0)**i for i in range(n)) )
>>> # now we make a matrix consisting of the coefficients
>>> # of the c_i in the nth degree polynomial P
>>> # coefficients of c_i evaluated at x_i
>>> m11 = P(x0 , x0, c, n).diff(c[0])
>>> m12 = P(x0 , x0, c, n).diff(c[1])
>>> m13 = P(x0 , x0, c, n).diff(c[2])
>>> # coefficients of c_i evaluated at x_i - h
>>> m21 = P(x0-h, x0, c, n).diff(c[0])
>>> m22 = P(x0-h, x0, c, n).diff(c[1])
>>> m23 = P(x0-h, x0, c, n).diff(c[2])
>>> # coefficients of c_i evaluated at x_i + h
>>> m31 = P(x0+h, x0, c, n).diff(c[0])
>>> m32 = P(x0+h, x0, c, n).diff(c[1])
>>> m33 = P(x0+h, x0, c, n).diff(c[2])
>>> # matrix of the coefficients is 3x3 in this case
>>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
现在我们有了系数矩阵，接下来我们形成右侧并通过求解 \(M\) 的逆来求解
>>> # matrix of the function values...actually a vector of right hand sides
>>> R = Matrix([[Fi], [Fim1], [Fip1]])
>>> # matrix form of the three equations for the c_i is M*X = R
>>> # solution directly inverting the 3x3 matrix M:
>>> X = M.inv() * R
>>> # note that all three coefficients make up the solution
>>> # the first derivative is coefficient c_1 which is X[1].
>>> print("The second-order accurate approximation for the first derivative is: ")
The second-order accurate approximation for the first derivative is:
>>> print(together(X[1]))
(F_{i+1} - F_{i-1})/(2*h)
这两个例子说明了如何使用 SymPy 直接找到二阶精度的一阶导数。第一个例子在所有三个点 \(x_i\)、\(x_{i+1}\) 和 \(x_{i+2}\) 处使用 \(x\) 和 \(F\) 的值，而第二个例子只使用 \(x\) 在两个点 \(x_{i-1}\) 和 \(x_{i+1}\) 处的值，因此效率更高。
从这两个简单的例子中，一个普遍的规则是，如果要使一阶导数的精度达到 \(O(h^{n})\)，则需要近似多项式中的 n+1 个函数值（这里通过函数 \(P(x,x0,c,n)\) 提供）。
现在让我们评估中心差结果的精度问题，看看我们如何确定它确实是二阶的。为此，我们取 \(dF/dx\) 的结果，代入更高阶多项式的多项式展开，看看我们得到什么。为此，我们创建了一组八个系数 d 并使用它们来执行检查
>>> d = symbols('c:8')
>>> dfdxcheck = (P(x0+h, x0, d, 8) - P(x0-h, x0, d, 8))/(2*h)
>>> print(simplify(dfdxcheck)) # so the appropriate cancellation of terms involving `h` happens
c1 + c3*h**2/6 + c5*h**4/120 + c7*h**6/5040
因此，我们看到导数实际上是 \(c_1\)，级数中的下一项是 \(h^2\) 阶的。
然而，当试图生成高阶导数近似值（例如 6 或 8）时，直接推广上述直接方法可能会变得相当繁琐，尽管该方法确实有效，并且使用当前方法比手动执行计算要省事得多。
正如我们在上面的讨论中所看到的，一阶导数的简单中心近似值只使用 \((x_{i},F_{i})\) 对的两个点值。这在遇到域中的最后一个点（例如在 \(i=N\) 处）之前工作得很好。由于我们的中心导数近似值将使用点 \((x_{N+1},F_{N+1})\) 处的数据，因此我们看到导数公式将不起作用。那么，该怎么做呢？好吧，一个简单的解决方法是为这个最后一点设计一个不同的公式，该公式使用我们确实有值的点。这就是所谓的向后差分公式。为了得到它，我们可以使用相同的方法，只是现在使用三个点 \((x_{N},F_{N})\)、\((x_{N-1},F_{N-1})\) 和 \((x_{N-2},F_{N-2})\) 并将近似值集中在 \((x_{N},F_{N})\) 处。以下是使用 SymPy 完成此操作的方法
>>> from __future__ import print_function
>>> from sympy import *
>>> x, xN, h = symbols('x, x_N, h')
>>> FN, FNm1, FNm2 = symbols('F_{N}, F_{N-1}, F_{N-2}')
>>> n = 8 # there are the coefficients c_0=Fi, c_1=dF/h, c_2=d**2F/h**2
>>> c = symbols('c:8')
>>> # define a polynomial of degree d
>>> def P(x, x0, c, n):
... return sum( ((1/factorial(i))*c[i] * (x-x0)**i for i in range(n)) )
现在我们构建一个矩阵，该矩阵包含 d 次多项式 P 中 \(c_i\) 的系数，即 \(c_i\) 在 \(x_i, x_{i-1},\) 和 \(x_{i+1}\) 处评估的系数
>>> m11 = P(xN , xN, c, n).diff(c[0])
>>> m12 = P(xN, xN, c, n).diff(c[1])
>>> m13 = P(xN , xN, c, n).diff(c[2])
>>> # coefficients of c_i evaluated at x_i - h
>>> m21 = P(xN-h, xN, c, n).diff(c[0])
>>> m22 = P(xN-h, xN, c, n).diff(c[1])
>>> m23 = P(xN-h, xN, c, n).diff(c[2])
>>> # coefficients of c_i evaluated at x_i + h
>>> m31 = P(xN-2*h, xN, c, n).diff(c[0])
>>> m32 = P(xN-2*h, xN, c, n).diff(c[1])
>>> m33 = P(xN-2*h, xN, c, n).diff(c[2])
接下来，我们构建 \(3 \times 3\) 系数矩阵
>>> M = Matrix([[m11, m12, m13], [m21, m22, m23], [m31, m32, m33]])
>>> # matrix of the function values...actually a vector of right hand sides
>>> R = Matrix([[FN], [FNm1], [FNm2]])
然后我们求解 \(M\) 的逆并写出 \(3 \times 3\) 系统的解。
三个方程关于 c_i 的矩阵形式为 \(M*C = R\)。解是通过直接求解 \(M\) 的逆得到的
>>> X = M.inv() * R
一阶导数是系数 \(c_1\)，即 \(X[1]\)。因此，一阶导数的二阶精度近似值为
>>> print("The first derivative centered at the last point on the right is:")
The first derivative centered at the last point on the right is:
>>> print(together(X[1]))
(-4*F_{N-1} + F_{N-2} + 3*F_{N})/(2*h)
当然，我们可以针对 \((x_{1},F_{1})\) 处的一组点中左端点的导数值，设计一个类似的公式，该公式基于 \((x_{2},F_{2})\) 和 \((x_{3},F_{3})\) 处的数值。
此外，我们注意到，可以在上面给出的示例中输出适合 Fortran、C 等的格式。
接下来，我们将展示如何执行这些以及许多其他导数离散化，但使用更有效的方法，该方法最初是由 Bengt Fornberg 提出的，现在已被整合到 SymPy 中。