如何计算python中数值导数的边界点？

Question

我正在尝试编写一个函数来获取任何通用函数/数字数组的导数。具体来说，我使用的是 Central difference formula。问题是，我无法计算导数的边界点，因为中心差分公式使用了超出范围的索引。我的代码如下

import numpy as np
n = 20000 # number of points in array
xs = np.linspace(start=-2*np.pi, stop=2*np.pi, num=n) # x values
y = np.array([np.sin(i) for i in xs]) # our function, sine

def deriv(f, h):
    """
    Calauclate the numerical derivative of any function
    :param f: numpy.array(float), the array of numbers we differentiate
    :param h: step size
    :rtype d: numpy.array(float)
    """
    d = np.zeros_like(f)
    # this loop misses the first and last points in f
    for i in range(1, f.shape[0]-1):
        # 2-point formula
        d[i] = (f[i+1] - f[i-1])/(2*h)

    return d

h = abs(xs[0] - xs[1]) # step size
y1 = deriv(y, h) # first derivative
y2 = deriv(y1, h) # second derivative
y3 = deriv(y2, h) # third derivative

当我绘图时 y,y1,y2,y3 你可以看到它在终点爆炸了

我尝试做的是将端点设置为 deriv 中最近的邻居，如下所示。虽然这适用于低阶导数（一阶和二阶），但在高阶导数（三阶和更高阶）时开始失效。

...
d = np.zeros_like(f)
for i in range(1, f.shape[0]-1):
    d[i] = (f[i+1] - f[i-1])/(2*h)

d[0] = d[1]
d[-1] = d[-2]
...

中间的导数，远离边界，计算正常。问题在于边界。

这里的边界条件应该怎么处理呢？不同的数值微分方案会比中心差分方案更好吗？

编辑： 我正在寻找一种解决此问题的通用方法，而不仅仅是一种可以应用于正弦函数或任何其他周期函数的方法，因为我已经习惯了在这里说明问题。

Answer 1

与编程问题相比，这更像是一个数值方法问题。

无论如何，如果您的函数具有周期性边界条件（它看起来是一个正弦波，所以在这种情况下您具有周期性）只需创建一个包含 2 个附加元素的新数组：新数组起始元素将是您的最后一个元素原始数组的元素和新数组的结束元素将是原始数组的开始元素。这是一种方法

f_periodic = np.zeros(f.size+2)
f_periodic[1:-1], f_periodic[0], f_periodic[-1] = f, f[-1], f[0]

您现在可以区分 f_periodic，其中 d[1] 和 d[-2] 将是边界上的正确导数值（忽略 d[0] 和 d[-1] ).

根据OP的新要求编辑...

对于更一般的边界条件，比如边界处的特定值，可以采用不同的方法：

使用重影值：

再次扩展新边界的函数和 extrapolate 值。根据数值微分的顺序，将需要更多的幽灵细胞。对于当前的方案，一个简单的线性外推就可以了（每个边界只需要 1 个重影值）：

f_new = np.zeros(f.size+2)
f_new[1:-1] = f
f_new[-1] = f[-2] + (f[-2]-f[-3])/(x[-2]-x[-3])*(x[-1]-x[-2])
f_new[0] = f[1] + (f[1]-f[2])/(x[1]-x[2])*(x[0]-x[1])

请注意，您还必须扩展 x。但是，由于您的间距恒定，因此只需使用 h 而不是空间差异，例如x[-2]-x[-3]。您现在可以微分 f_new 并且您将获得边界上导数的一阶近似值（因为您使用线性外推来找到重影值）。

在边界上使用前向和后向方案

我不会在这里显示代码，但基本上您需要分别使用边界值和左右边界的右（向前）或左（向后）值来区分。这是一阶近似值。

Answer 2

边界点可以使用2阶的前向后向微分方案。本质上，我们知道

(f(x+h)-f(x))/h = f'(x) + h/2*f''(x) + O(h²)                       (I)

和

(f(x+2h) - 2f(x+h) + f(x))/h² = f''(x+h) + O(h²) = f''(x) + O(h)   (II)

用最后一个用二阶导数代替一阶项，即计算(I)-h/2*(II)得到

(-1/2*f(x+2h) + 2*f(x+h) -3/2*f(x))/h = f'(x) + O(h²)

请注意，一阶导数中的 O(h²) 误差通常会导致除差的第二次迭代中的 O(h) 误差和第三次迭代中的 O(1) 误差。有人可能会争辩说误差项适当地取消了，但这只会发生在内部点上，单边导数将 "spoil" 该模式与边界的距离增加。

如何计算python中数值导数的边界点？

How to calculate boundary points of numerical derivative in python?

python

arrays

math

numerical-methods