odeint浮点运算

Question

我对使用scipy.integrate.odeint函数理解浮点算术很感兴趣。

我正在处理的案例如下

# data
omega = 136 # rad/s
d = 75 # Nm/s
k = 390000 # N/m
m = 4 # kg
n = 1000 # 
t_0 = 1 # s
t_1 = 5.5 # s
Y = 0.05 # m

# time
t = np.linspace(t_0, t_1, n)

# initial condition
x_0 = np.array([0, 0])

# first function
def fun(x, t, k, d, m, Y, omega):
    y = Y*np.sin(omega*t)
    return np.array([x[1], (y - k*x[0] - d*x[1]) / m])

# second function
def fun2(x, t, k, d, m, Y, omega):
    y = Y*np.sin(omega*t)
    return np.array([x[1], (-k*x[0] - d*x[1] + y)/m])

# results
res = odeint(fun, x_0, t, args=(m, k, d, Y, omega))
res2 = odeint(fun2, x_0, t, args=(m, k, d, Y, omega))

注意，这两个函数在数学上是相同的。唯一的区别是数值运算的顺序。

我想更好地理解结果res - res2中的区别，即：

array([[  0.00000000e+00,   0.00000000e+00],
       [ -1.95628215e-22,   1.91508855e-19],
       [  6.33676391e-19,  -2.16307730e-17],
       ..., 
       [ -8.50849113e-10,   3.04613004e-09],
       [ -8.49843242e-10,  -9.43460353e-10],
       [ -1.00314946e-09,   4.45237878e-09]])

但它应该是一个零的数组。

Answer 1

您所看到的是浮点运算(加法和乘法)的不同顺序的结果。查看经典论文http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html或维基百科https://en.wikipedia.org/wiki/Floating_point#Floating-point_arithmetic_operations

即使差异可能很小，因为odeint解决了给定的精度，所以您留下的差异可以达到该treshold。

我对跨SE站点的副本了解不多，但这是一个密切相关的问题：https://scicomp.stackexchange.com/questions/10506/number-of-equations-and-precision-of-scipys-integrate-odeint或这一个scipy.integrate.odeint fails depending on time steps，由SciPy核心开发人员回复。