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()

虽然这两个代码在数学上与实数等价，但由于固定大小的浮点精度，它们并不等价。事实上，与{}和{}相比，{}和{}在{}较大时，都会产生巨大的值。不幸的是，当两个固定大小的浮点数相加时，会出现精度损失。当添加的数字的数量级非常不同时，精度损失可能是巨大的（如果不是关键的话）。例如，在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数值运算是一个良好的开端

在这里，一个相对较好的解决方案是只使用第二个代码而不是第一个代码

Answer 2

使用：

aux = np.linspace(0, L, 2000).astype(np.longdouble)

稍微好一点