迭代和矢量化正弦函数的计算值存在差异

Question

我正在使用包含纬度值的 NetCDF 文件，并且我正在尝试计算纬度正弦值的倒数。我写的代码如下：

import xarray as xr
import numpy as np

data = xr.open_dataset('ERA5_2000_01_01.nc')
lat = data.latitude.values
omega = 7.292e-5

f1 = np.array([2 * omega * np.sin(x * np.pi / 180.0) for x in lat])
f2 = 2 * omega * np.sin(lat * np.pi / 180.0)

diff = 1/f1 - 1/f2
print(diff)

差异是

[-0.00018776 -0.00050296 -0.00059619 -0.00035761  0.00011444 -0.00065717
 -0.00024923 -0.00060329  0.00018545 -0.00014504 -0.00064833 -0.00043263
  0.00061779  0.00013936 -0.00046115 -0.00041766 -0.0009715   0.0005565
 -0.00034907 -0.00123883 -0.00013027  0.00075469  0.00162418 -0.00066007
 -0.00037552 -0.00189229 -0.00101494  0.00089398  0.00063323  0.00030524
  0.00051748  0.00012772  0.00051868  0.00344283 -0.002541    0.00145573
         nan -0.00145573  0.002541   -0.00344283 -0.00051868 -0.00012772
 -0.00051748 -0.00030524 -0.00063323 -0.00089398  0.00101494  0.00189229
  0.00037552  0.00066007 -0.00162418 -0.00075469  0.00013027  0.00123883
  0.00034907 -0.0005565   0.0009715   0.00041766  0.00046115 -0.00013936
 -0.00061779  0.00043263  0.00064833  0.00014504 -0.00018545  0.00060329
  0.00024923  0.00065717 -0.00011444  0.00035761  0.00059619  0.00050296
  0.00018776]

在此代码中，我通过迭代纬度值并计算每个纬度值的正弦来计算

f1

。相反，我直接在矢量化操作中使用整个纬度值数组来计算

f2

。然而，我注意到当我计算差异时，f1 和 f2 的结果会产生不同的结果

diff = 1/f1 - 1/f2

。

我试图理解为什么两个结果不同。用于重现此示例的数据可在此link中找到。

Answer 1

这是由于两件事：

• lat

  是

  np.float32

  类型，因此它只有大约 6 位精度。

• 1/f1 - 1/f2

  会导致灾难性取消，因此在此操作期间大部分精度都会丢失。

您可以通过以下方式提高精度：

• 使用

  lat

   将 

  lat = lat.astype(np.float64)
转换为 64 位数组
• 使用另一种数值方法来计算差异，或者使用扩展精度，例如双精度（这不是灵丹妙药，而是精度和速度之间的良好平衡）。