我正在尝试寻找数据集(x,y)的高阶导数。X和y是长度为N的一维阵列。
假设我将它们生成为:
xder0=np.linspace(0,10,1000)
yder0=np.sin(xder0)我定义了导数函数,它接受2个数组(x,y)并返回(x1,y1),其中y1是在每个索引上计算的导数:(yi+1-yi)/(xi+1-xi)。x1只是xi+1和xi的意思
下面是执行此操作的函数:
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号订单变得非常嘈杂!正如预期的那样,早期的导数都是正弦或cos函数。然而,7阶是一个有噪声的正弦。因此,在那次爆炸之后,所有的导数。

知道这是怎么回事吗?
发布于 2015-12-04 14:13:49
这是一个众所周知的稳定性问题,使用等间距点进行数值插值。请在http://math.stackexchange.com上阅读答案。
为了解决这个问题,你必须使用非等间距的点,就像Lagendre polynomial的根一样。该不稳定性是由于在边界处的信息不可用而发生的,因此需要在边界处更多地集中点,例如根据拉金德多项式的根或具有类似性质的其它如切比雪夫多项式。
https://stackoverflow.com/questions/34061541
复制相似问题