Python numpy gradient用法及代碼示例

本文簡要介紹 python 語言中 numpy.gradient 的用法。

用法:

numpy.gradient(f, *varargs, axis=None, edge_order=1)

返回 N 維數組的梯度。

使用內部點中的二階精確中心差異和邊界處的一階或二階精確one-sides(前向或後向)差異計算梯度。因此，返回的梯度具有與輸入數組相同的形狀。

參數：

f： array_like

包含標量函數樣本的 N 維數組。

varargs： 標量或數組列表，可選

f 值之間的間距。所有尺寸的默認單一間距。可以使用以下方式指定間距：

1. 用於指定所有維度的樣本距離的單個標量。

2. N個標量為每個維度指定一個恒定的樣本距離。即 dx, dy, dz, ...

3. N 個數組，用於指定沿 F 的每個維度的值的坐標。數組的長度必須與對應維度的大小匹配

4. 具有 2. 和 3 含義的 N 個標量/數組的任意組合。

如果給定了軸，則可變參數的數量必須等於軸的數量。默認值：1。

edge_order： {1, 2}，可選

梯度使用N-th 計算邊界處的精確差異。默認值：1。

axis： 無或int 或整數元組，可選

梯度僅沿給定的一個或多個軸計算默認值(axis = None)是計算輸入數組的所有軸的梯度。軸可能是負數，在這種情況下，它從最後一個軸計數到第一個軸。

返回：

gradient： ndarray 或 ndarray 列表

對應於 f 對每個維度的導數的 ndarrays 列表(如果隻有一個維度，則為單個 ndarray)。每個導數具有與 f 相同的形狀。

注意：

假設\(f\in C^{3}\)(即\(f\)至少有3個連續導數)並讓\(h_{*}\)是一個非齊次步長，我們最小化真實梯度和線性組合估計之間的“consistency error”\(\eta_{i}\)相鄰的grid-points：

\[\eta_{i} = f_{i}^{\left(1\right)} - \left[ \alpha f\left(x_{i}\right) + \beta f\left(x_{i} + h_{d}\right) + \gamma f\left(x_{i}-h_{s}\right) \right]\]

通過用它們的泰勒級數展開替換\(f(x_{i} + h_{d})\) 和\(f(x_{i} - h_{s})\)，這轉化為求解以下線性係統：

\[\begin{split}\left\{ \begin{array}{r} \alpha+\beta+\gamma=0 \\ \beta h_{d}-\gamma h_{s}=1 \\ \beta h_{d}^{2}+\gamma h_{s}^{2}=0 \end{array} \right.\end{split}\]

\(f_{i}^{(1)}\) 的結果近似值如下：

\[\hat f_{i}^{(1)} = \frac{ h_{s}^{2}f\left(x_{i} + h_{d}\right) + \left(h_{d}^{2} - h_{s}^{2}\right)f\left(x_{i}\right) - h_{d}^{2}f\left(x_{i}-h_{s}\right)} { h_{s}h_{d}\left(h_{d} + h_{s}\right)} + \mathcal{O}\left(\frac{h_{d}h_{s}^{2} + h_{s}h_{d}^{2}}{h_{d} + h_{s}}\right)\]

值得注意的是，如果\(h_{s}=h_{d}\)(即數據是均勻分布的)我們會找到標準的二階近似值：

\[\hat f_{i}^{(1)}= \frac{f\left(x_{i+1}\right) - f\left(x_{i-1}\right)}{2h} + \mathcal{O}\left(h^{2}\right)\]

通過類似的過程，可以導出用於邊界的前向/後向近似值。

參考：

1：

Quarteroni A.、Sacco R.、Saleri F. (2007) 數值數學(應用數學課本)。紐約：斯普林格。

2：

Durran D. R. (1999) 地球物理流體動力學波動方程的數值方法。紐約：斯普林格。

3：

Fornberg B. (1988) 任意間隔網格上的有限差分公式的生成，計算數學 51，不。 184：699-706。 PDF 。

例子：

>>> f = np.array([1, 2, 4, 7, 11, 16], dtype=float)
>>> np.gradient(f)
array([1. , 1.5, 2.5, 3.5, 4.5, 5. ])
>>> np.gradient(f, 2)
array([0.5 ,  0.75,  1.25,  1.75,  2.25,  2.5 ])

間距也可以用一個數組來指定，該數組表示值 F 沿維度的坐標。例如均勻間距：

>>> x = np.arange(f.size)
>>> np.gradient(f, x)
array([1. ,  1.5,  2.5,  3.5,  4.5,  5. ])

或非統一的：

>>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=float)
>>> np.gradient(f, x)
array([1. ,  3. ,  3.5,  6.7,  6.9,  2.5])

對於二維數組，返回將是按軸排序的兩個數組。在這個例子中，第一個數組代表行中的梯度，第二個代表列方向的梯度：

>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float))
[array([[ 2.,  2., -1.],
       [ 2.,  2., -1.]]), array([[1. , 2.5, 4. ],
       [1. , 1. , 1. ]])]

在此示例中，還指定了間距：axis = 0 的統一和axis = 1 的非統一

>>> dx = 2.
>>> y = [1., 1.5, 3.5]
>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), dx, y)
[array([[ 1. ,  1. , -0.5],
       [ 1. ,  1. , -0.5]]), array([[2. , 2. , 2. ],
       [2. , 1.7, 0.5]])]

可以使用 edge_order 指定如何處理邊界

>>> x = np.array([0, 1, 2, 3, 4])
>>> f = x**2
>>> np.gradient(f, edge_order=1)
array([1.,  2.,  4.,  6.,  7.])
>>> np.gradient(f, edge_order=2)
array([0., 2., 4., 6., 8.])

axis 關鍵字可用於指定計算梯度的軸的子集

>>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=float), axis=0)
array([[ 2.,  2., -1.],
       [ 2.,  2., -1.]])