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波的纬向波数。这种波可以被认为像鼓上的波，例如，只有在一个球形的环境中，大气中。 我使用的是python3.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代表纬度位置。 为了计算导数，我使用如下的numpy梯度函数：

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

This是计算经向波数（l）平方的公式，其中ω是地球自转，φ是纬度位置，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