有限差分法近似导数¶
介绍¶
导数的有限差分近似在数值分析和计算物理中非常重要。在本教程中，我们展示了如何使用 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}\)。 因此，问题是如何在约束条件下生成 F 导数的近似值，即我们使用有限对集合 \((x_{i},F_{i})\) 的子集，大小为 N。
接下来是使用 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]])
现在我们构造一个矩阵，该矩阵由 n 次多项式 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\) 和 \(F\) 在所有三个点 \(x_i\)、\(x_{i+1}\) 和 \(x_{i+2}\) 处的值，而第二个例子仅使用了 \(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)) )
现在我们构造一个矩阵，该矩阵由 \(c_i\) 在 d 次多项式 P 系数中组成，这些系数在 \(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 imes 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 imes 3\) 系统的解。
三个方程的矩阵形式为 \(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 的更高效的方法。
微积分-有限差分