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;DR：问题来自于数值不稳定性。

首先，这是一个简化的代码，出现完全相同的问题（具有不同的值）：

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。不幸的是，当两个固定大小的浮点数相加时，精度会下降。当相加数字的数量级非常不同时，精度损失可能会很大（如果不是很严重的话）。例如，在 Python 和主流平台上（对于 IEEE-754 64 位浮点数也是如此），由于数字表示的精度有限，这是正确的。因此也是正确的，而不是正确的（结果值为 0.1）。np.sin(x)np.cos(x)0.1 + 1e20 == 1e20(0.1 + 1e20) - 1e20 == 0.00.1 + (1e20 - 1e20) == 0.0浮点加法既不是结合的也不是可交换的. 在这种特定情况下，准确度可以达到阈值，此时结果中不再有重要数字。有关浮点精度的更多信息，请阅读这篇文章。

关键是你在减去浮点数时应该非常小心。一个好的解决方案是放置括号，以便添加/减去的值具有相同的数量级。可变大小和更高的精度也有帮助。但是，最好的解决方案是分析算法的数值稳定性。例如，研究算法中使用的数值运算的条件数是一个好的开始。

在这里，一个比较好的解决方案就是使用第二个代码而不是第一个。