I have an array x of shape (N, d).
If I create the arrays r_1 and r_2 as follows, the outputs should be equal, which is true.
import numpy as np
N = 100
d = 2
x = np.random.rand(N, d)
r_1 = x.reshape(N, 1, d) - np.repeat(x[np.newaxis, :, :], N, axis=0)
r_2 = x.reshape(N, 1, d) - x.reshape(1, N, d)
print((r_1 == r_2).all())
> True
However, if I use the arrays in a deterministic simulation, the two methods yield different end results due to (I assume) floating point errors. I did some digging into the simulation and everything is equal until the eighth iteration, in which I get differences such as:
Method 1: 0.4759682049387506
Method 2: 0.47596820493875036
And these floating point errors can propagate to quite different end results of the simulation:
Method 1:
E_kin: 11651.935893586868
E_pot: -22091.50550903025
Method 2:
E_kin: 12059.81250826829
E_pot: -21037.60706699791
What is the underlying difference between these two arrays that computations with them don't end up with the same results? Everything in the simulation is held constant except for whether or not method one or two is used for computing r.
Edit: This is using numpy version 1.18.1.