马修特征值

如何解决马修特征值

我试图获得特定问题的 Mathieu 特征值。我在获取它们时没有任何问题，并且我已经阅读了 Scipy 中关于这些函数的文档。问题是我知道我得到的分数是不正确的。我获取所需特征值的脚本如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import mathieu_a,mathieu_b,mathieu_cem,mathieu_sem


M = 1.0
g = 1.0
l = 1.0
h = 0.06
U0 = M * g * l
q = 4 * M * l**2 * U0 / h**2

def energy(n,q):
    if n % 2 == 0:
        return (h**2 / (8 * M * l**2)) * mathieu_a(n,q) + U0
    else:
        return (h**2 / (8 * M * l**2)) * mathieu_b(n + 1,q) + U0

n_list = np.arange(0,80,1)
e_n = [energy(i,q) for i in n_list]

plt.plot(n_list,e_n,'.')

这些值的结果图是 this one。有一个区域似乎是“噪音”或数字错误，我知道这些跳跃不能发生。实际上，在 x= 40 到 x > 40 附近，这些点应该表现得像两个连续点的阶梯，类似于在 70

看到Mathieu函数的实现有一些问题，见here。但这是六年前的事了！在回答这个问题时，他们使用了 Python 的 NAG 库，但它并不是完全开源的。

有没有什么方法可以让我仍然使用 Scipy 的这些功能而不会出现这个问题？还是与我用来获取 Mathieu 特征值的精度有关？