为何NumPy数组time**3与(time/1)**3计算积分结果存在差异？

我推导了积分的解析解，需计算 ( Q = Q_{max} \times (1 - (\frac{time}{t_0} - 1)^2) ) 的时间积分，编写了如下Python代码：

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))

但Q_integral_A得到错误结果，尝试后发现将time改为time/1的Q_integral_B能得到正确结果：

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

请问这是怎么回事？为何Q_integral_A与Q_integral_B存在差异？

使用版本：Python 3.8.5、NumPy 1.20.3、Spyder 4.2.5

解答 #

核心原因是整数数组的溢出问题：

• 你的环境中，np.arange生成的time数组是32位整数类型（int32），int32的最大值仅为2147483647（约2×10⁹）。当计算time**3时，time的最大值为43140，其立方值约为8×10¹³，远超int32的范围，触发整数溢出，导致time**3的数值完全错误，最终使得Q_integral_A计算结果偏差。
• time/1会将int32数组强制转换为64位浮点数类型（float64），float64的精度足以容纳8×10¹³这样的数值，立方运算不会溢出，因此(time/1)**3能得到正确结果，Q_integral_B的计算自然正确。

你可以通过打印time.dtype和(time**3).dtype验证类型。除了用time/1转换类型外，还可以在生成数组时指定 dtype 避免溢出：

# 指定为64位整数
time = np.arange(0,2*t_0,dt, dtype=np.int64)
# 或者直接生成浮点数组
time = np.arange(0,2*t_0,dt, dtype=np.float64)

内容的提问来源于stack exchange，提问作者Sarah_NW