python中气候科学计算的数值精度

Question

我目前正在尝试为我的硕士论文重建一篇论文(https://www.researchgate.net/publication/309723672_Evidence_for_wave_resonance_as_a_key_mechanism_for_generating_high-amplitude_quasi-stationary_waves_in_boreal_summer)的发现。我计算了一段时间内Rossby波(大气中的一种特定气流)经向波数平方的经向(北度)分布。该值仅取决于平均纬向(西度和东度)风及其一、二经向导数，以及二维Rossby波的纬向波数。例如，这种波可以被认为是鼓上的波，只有在球形环境中，大气中才会这样。我使用的是python 3.6.5，我怀疑问题出在数值精度上，但是我不确定。

我已经阅读了其他关于数字精度的帖子，并遇到了这个例子：python sine and cosine precision。然而，我还没有尝试过这一点，因为我试图避免编写自己的三角函数。另外，由于我必须处理相当多的数据，所以我尽量不让我的代码变慢。通过实验，我发现在三角函数方面，数学库并不比numpy库更精确。

下面是我关心的代码片段：

Lat = np.linspace(0,90,37)
MeridWN = np.zeros((29,36), dtype='float64')        
        ######################################################
        #define Meridional wavenumber, l^2
        for i in range(5,28,10):
            for j in range(36):
                MeridWN[i,j] = (((2*EarthRot*np.cos(Lat[j]*np.pi/180.0)**3.0)/(EarthRad*ZonMeanZonWiNH[j]))-
                   ((np.cos(Lat[j]*np.pi/180.)**2.)/(EarthRad**2.0*ZonMeanZonWiNH[j]))*
                   ZonMeanZonWiMeridGradGrad[j]+
                   ((np.sin(Lat[j]*np.pi/180.)*np.cos(Lat[j]*np.pi/180.))/(EarthRad**2.0*ZonMeanZonWiNH[j]))*
                   ZonMeanZonWiMeridGrad[j]+(1./(EarthRad**2.0))-(ZonWN[i]/EarthRad)**2.0)
                MeridWNMerge[i,x,j] = MeridWN[i,j]

索引i表示一系列的纬度波数，x是日期(这段代码来自一个更大的循环，该循环在几天内运行)，j是纬度位置。为了计算导数，我使用了如下数值梯度函数：

ZonMeanZonWiMeridGrad = np.gradient(ZonMeanZonWiNH,np.linspace(0,90,37))
ZonMeanZonWiMeridGradGrad = np.gradient(ZonMeanZonWiMeridGrad,np.linspace(0,90,37))

This是计算子午线波数(l)的平方的公式，其中Omega是地球自转，φ是纬度位置，a是地球半径，U是纬向平均值，k是纬向波数，在我的例子中是一个从5.5到8.5的数组。

This是我从6月到8月的纬向平均纬向风场(底部)与论文中的纬向平均风场(顶部)的比较，表明我们有相同的数据，这不是问题。颜色尺度略有不同，但风廓线最突出的特征非常相似，微小的差异不应产生如此不同的profile of the meridional wavenumber (for k = 7)，其中来自纸张的图形再次位于顶部，而我的图形位于底部。在这里，色标再次不同，但仍然应该捕捉到大量的结构相似之处。正如您所看到的，我处理的数字非常小，这导致我怀疑有一个数字不精确的代码。

如果你愿意，我可以上传我的整个代码，但是我认为对于精度的讨论这是足够的。我试着解决这个问题大约2周了，尝试了我的代码中所有不同的更改，其中一些是好的更改，但是没有一个给出了想要的输出。

先谢谢你，

托马斯

Answer 1

如果是数值精度的问题，那么算法不是问题，而是你使用的变量类型。如果您切换到更长的浮点表示形式，如dtype='c16l‘( 128-复数浮点数。https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html