【问题标题】:Derivatives blow up in python衍生物在 python 中爆炸
【发布时间】:2015-12-03 08:54:55
【问题描述】:

我正在尝试查找数据集 (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 阶变得非常嘈杂!正如预期的那样,早期的导数都是正弦或余弦函数。然而,7 阶是一个嘈杂的正弦波。因此,在那之后的所有衍生品都爆炸了。

知道发生了什么吗?

【问题讨论】:

  • 近似推导的差分公式都是病态的,因为我们减去两个比较大的(x[i], x[i+1], y[i], y[i + 1] ) 数字多次以获得相对较小的差异。 order = 7 似乎是你的下车点。
  • 有什么想法吗?
  • 寻找一个处理数字部分的库。 numdifftools 浮现在脑海中。如您所见,直接的、不稳定的方法只能让您走这么远。
  • 一种可能性是为您的数据拟合一些函数,然后您可以对其进行分析导数。
  • 我建议阅读Numerical Recipes 的第 5.7 节。它提供了多种计算数值导数的方法,同时最大限度地减少了稳定性问题。您每月可以访问几页。我认为应该足以满足您的目的。

标签: python numpy numerical-methods derivative


【解决方案1】:

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

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

【讨论】:

    猜你喜欢
    • 2023-03-12
    • 2012-07-13
    • 1970-01-01
    • 1970-01-01
    • 2014-02-04
    • 2021-07-17
    • 1970-01-01
    • 2010-09-30
    • 2012-04-07
    相关资源
    最近更新 更多