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 またはこの1つですscipy.integrate.odeint fails depending on time steps SciPyコア開発者からの返信あり。