python - Numpy 返回分析函数的意外结果

Question

当我尝试计算 d_j(x) ，定义如下，基于 Numpy 的算法会导致意外的值。我相信它与数值精度有关，但我不确定如何解决这个问题。
功能是:

在哪里

和

当 j>10 时代码失败.例如，当 j=16 ，函数d_j(x)从 x=1 附近返回错误的值，而预期的结果是一条平滑的、几乎是周期性的曲线。0<x<1.5的图:

代码是:

#%%
import numpy as np
import matplotlib.pyplot as plt

#%%
L = 1.5  # Length [m]

def eta(j):
    if j == 1:
        return 1.875 / L
    if j > 1:
        return (j - 0.5) * np.pi / L


def D(etaj):
    etajL = etaj * L
    return (np.cos(etajL) + np.cosh(etajL)) / (np.sin(etajL) - np.sinh(etajL))


def d(x, etaj):
    etajx = etaj * x
    return np.sin(etajx) - np.sinh(etajx) + D(etaj) * (np.cos(etajx) - np.cosh(etajx))


#%%
aux = np.linspace(0, L, 2000)
plt.plot(aux, d(aux, eta(16)))
plt.show()

Answer 1

TL;博士:问题来自数值不稳定性 .
首先，这是一个简化的代码，出现完全相同的问题(具有不同的值):

x = np.arange(0, 50, 0.1)
plt.plot(np.sin(x) - np.sinh(x) - np.cos(x) + np.cosh(x))
plt.show()

这是另一个没有出现问题的示例:

x = np.arange(0, 50, 0.1)
plt.plot((np.sin(x) - np.cos(x)) + (np.cosh(x) - np.sinh(x)))
plt.show()

虽然这两个代码在数学上与实数等效，但由于固定大小的浮点精度，它们并不等效。确实，np.sinh(x)和 np.cosh(x)结果都在 x 时的巨大值(value)大与 np.sin(x) 相比和 np.cos(x) .不幸的是，当两个固定大小的浮点数相加时，精度会下降。当相加数字的数量级非常不同时，精度损失可能会很大(如果不是很严重的话)。例如，在 Python 和主流平台上(IEEE-754 64 位浮点数也是如此)，0.1 + 1e20 == 1e20由于数字表示的精度有限，因此为真。因此(0.1 + 1e20) - 1e20 == 0.0也是如此，而 0.1 + (1e20 - 1e20) == 0.0不是真的(结果值为 0.1)。 浮点加法既不是结合也不是交换 .在这种特定情况下，准确度可以达到阈值，此时结果中不再有重要数字。有关浮点精度的更多信息，请阅读 this post .
关键是你应该减去浮点数时要非常小心 .一个好的解决方案是放置括号，以便添加/减去的值具有相同的数量级。可变大小和更高的精度也有帮助。但是，最好的解决方案是analyses the numerical stability你的算法。比如学习condition number算法中使用的数值运算是一个良好的开端。
在这里，一个比较好的解决方案就是使用第二个代码而不是第一个。