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

Question

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

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]
...

中间的导数，远离边界，计算得很好。问题在于边界。

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

编辑：我正在寻找一种通用的方法来解决这个问题，而不仅仅是一种可以应用于正弦函数或任何其他周期函数的方法，就像我在这里用来说明问题一样。

标签： pythonarraysmathnumerical-methods

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的新要求之后编辑...

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

1. 使用幻影值：

再次扩展函数并推断新边界的值。根据数值微分的顺序，将需要更多的幻影单元。对于当前方案，可以进行简单的线性外推（每个边界只需要 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，您将获得边界上导数的一阶近似值（因为您使用线性外推法来查找重影值）。

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

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