Why does np.array**3 lead to a different solution than (np.array/1)**3 in my code?-CodePudding

I wrote this code to calculate the integral in time of: Q = (Q_max*(1 - (time/t_0 - 1)**2)) which I derived analytically.

import numpy as np
Q_max       = 400                   # W m-2
t_0         = 6*3600                # seconds
dt          = 60                    # seconds
time        = np.arange(0,2*t_0,dt) 
Q_integral_A = Q_max*((time)**2/(t_0) - (time)**3/(3*(t_0)**2))

However, I found that Q_integral_A gives the wrong solution. After trying a lot of stuff, I found out that doing the following leads to the right solution (dividing the second "time" by 1):

Q_integral_B = Q_max*((time)**2/(t_0) - (time/1)**3/(3*(t_0)**2))

What is happening here? Why is there a difference between Q_integral_A and Q_integral_B?

screenshot of output

Versions used: Python 3.8.5 Numpy 1.20.3 Spyder 4.2.5

CodePudding user response：

I looked into problem myself, and I get the same result. So at first time is a int32, but when you do time / 1 it becomes float64. It shouldn't bring problems by itself, but time contains some big numbers, and raising them to 3rd power results in overflowing (here's what i get), but it doesn't effect float, because it works in different way.

To solve it just pass dtype="int64" time = np.arange(0, 2*t_0, dt, dtype = "int64"), but it won't solve the problem for even larger numbers.

CodePudding user response：

I suspect what is going on is that you are on Windows so:

time = np.arange(0,2*t_0,dt)

Defaults to dtype=np.int32, then the / operator results in promotion to np.float64.

I can reproduce your error if I use

time = np.arange(0,2*t_0,dt, dtype=np.int32)

And I bet it would be fixed for you if you use:

time = np.arange(0,2*t_0,dt, dtype=np.int64)

In general, you should be explicit with dtypes. If only to get more reproducible behavior.