python - 如何在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 的新要求进行编辑...

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

1. 使用幽灵值:

再次扩展功能和 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 进行微分，并且您将获得边界上导数的一阶近似值(因为您使用了线性外推法来查找鬼值)。

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

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