衍生品在python中爆炸

Question

我试图找到数据集(x,y)的高阶导数.x和y是长度为N的1D阵列.

假设我将它们生成为:

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个订单变得非常嘈杂!早期的衍生品都是预期的正弦或余弦函数.然而,第7顺序是一个嘈杂的正弦.因此,之后所有衍生品都会爆炸.

知道发生了什么事吗？

Answer 1

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

要解决这个问题，您必须使用非等距点，例如拉让德多项式的根。不稳定的发生是由于边界处的信息不可用，因此需要在边界处集中更多的点，根据拉根德多项式或其他具有类似属性（例如切比雪夫多项式）的根。