Really weird issue for a simple calculation with Python code

Question

I have the following two functions to do some simple physics calculation.

def force(r):
    if r == 0:
        return 0
    else:
        return 1/r**13 - 1/r**7

def calculate_forces_and_positions(x, l):
    n = len(x)  # Current number of points
    forces = np.zeros_like(x)
    # Calculate forces
    for i in range(n):
        left_neighbour =  (x[i] - x[i-1]) / l[i] if i > 0 else x[i] / l[i]
        right_neighbour = (x[i+1] - x[i]) / l[i] if i < n - 1 else 0
        forces[i] = -force(right_neighbour) + force(left_neighbour)
        print(f"n= {i} and L ={left_neighbour}  R = {right_neighbour} Force{i} = {forces[i]} force-R = {force(right_neighbour)} force-L= {force(left_neighbour)}")
    # Calculate positions
    print(forces)
    x_new = x + dt * forces
    return forces, x_new

Then I pass x[1,2,3] and l[1,1,1]. Surprisingly, I get this output:

n= 0 and L =0.9900990099009901  R = 0.9900990099009901 Force0 = 0 force-R = 0.06595792832627945 force-L= 0.06595792832627945
n= 1 and L =0.9900990099009901  R = 0.9900990099009901 Force1 = 0 force-R = 0.06595792832627945 force-L= 0.06595792832627945
n= 2 and L =0.9900990099009901  R = 0 Force2 = 0 force-R = 0 force-L= 0.06595792832627945

So for n = 2 Force2 = 0, but it is clearly force-R + force-L = 0.06. I don't know where things go wrong and I am really confused. How can I fix it?

Possibly related: *[What is the best way to compare floats for almost-equality in Python?](https://stackoverflow.com/questions/5595425/what-is-the-best-way-to-compare-floats-for-almost-equality-in-python)* — Peter Mortensen, Jun 01 '23 at 20:39
If `x` value is `[1,2,3]`, then `n==3` and range(n) will yield 0,1,2 (i.e. these are values of `i`. — buran, Jun 01 '23 at 20:42
As with @KellyBundy I can't reproduce this behaviour. Also, there is no value for `dt` so the code fails on running at the second-last line of `calculate_forces_and_positions(x, l)` — Simon, Jun 01 '23 at 20:46

score 3 · Accepted Answer · answered Jun 01 '23 at 20:42

np.zeros_like is creating an int64 array. In your example, when you convert 0.06 to an integer, you get zero. You need to specify the dtype if you want floats.

    forces = np.zeros_like(x, dtype=float)

Output:

n= 0 and L =0.9900990099009901  R = 0.9900990099009901 Force0 = 0.0 force-R = 0.06595792832627945 force-L= 0.06595792832627945
n= 1 and L =0.9900990099009901  R = 0.9900990099009901 Force1 = 0.0 force-R = 0.06595792832627945 force-L= 0.06595792832627945
n= 2 and L =0.9900990099009901  R = 0 Force2 = 0.06595792832627945 force-R = 0 force-L= 0.06595792832627945

And, as a side note, you are not passing l[1,1,1]. You are passing [1.01,1.01,1.01].

Really weird issue for a simple calculation with Python code

1 Answers1