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 или wikipedia https://en.wikipedia.org/wiki/Floating_point#Floating-point_arithmetic_operations

Несмотря на то, что разница может быть незначительной, поскольку odeint решает данную точность, вы остаетесь с разницей, которая может быть до этого.

Я не знаю много о дубликатах через сайты 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.