python - 衍生品在 python 中爆炸

Question

我试图找到数据集 (x,y) 的高阶导数。 x 和 y 是长度为 N 的一维数组。

假设我将它们生成为:

xder0=np.linspace(0,10,1000)
yder0=np.sin(xder0)

我定义了导数函数，它接受 2 个数组 (x,y) 并返回 (x1, y1)，其中 y1 是在每个索引处计算的导数:(y[i+1]-y[i])/(x[i+1]-x[i])。 x1 只是 x[i+1] 和 x[i] 的平均值

这是执行此操作的函数:

def deriv(x,y):
    delx =np.zeros((len(x)-1), dtype=np.longdouble)
    ydiff=np.zeros((len(x)-1), dtype=np.longdouble)
    for i in range(len(x)-1):
            delx[i]  =(x[i+1]+x[i])/2.0
            ydiff[i] =(y[i+1]-y[i])/(x[i+1]-x[i])
    return delx, ydiff

现在计算一阶导数，我称这个函数为:

xder1, yder1 = deriv(xder0, yder0)

与二阶导数类似，我将此函数称为以一阶导数作为输入:

xder2, yder2 = deriv(xder1, yder1)

然后继续:

xder3, yder3 = deriv(xder2, yder2)
xder4, yder4 = deriv(xder3, yder3)
xder5, yder5 = deriv(xder4, yder4)
xder6, yder6 = deriv(xder5, yder5)
xder7, yder7 = deriv(xder6, yder6)
xder8, yder8 = deriv(xder7, yder7)
xder9, yder9 = deriv(xder8, yder8)

当我到达第 7 阶后发生了一些奇怪的事情。第 7 阶变得非常嘈杂!正如预期的那样，较早的导数都是正弦函数或余弦函数。然而，七阶是一个嘈杂的正弦。因此，在那次爆炸之后的所有衍生品都崩溃了。

知道发生了什么吗？

Answer 1

这是使用等距点 进行数值插值的众所周知的稳定性问题。在 http://math.stackexchange.com 阅读答案.

要克服这个问题，您必须使用非等距点，例如Lagendre polynomial 的根.由于边界处的信息不可用，会出现不稳定性，因此需要在边界处更多地集中点，根据拉让德多项式或其他具有类似性质的多项式(如切比雪夫多项式)的根。