python - 使用NumPy的Mittag-Leffler函数的不稳定性

在试图重现the plot on Wolfram MathWorld和帮助this SO question时，我遇到了一些我不理解的数值不稳定性：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

def MLf(z, a):
    """Mittag-Leffler function
    """
    k = np.arange(100).reshape(-1, 1)
    E = z**k / gamma(a*k + 1)
    return np.sum(E, axis=0)

x = np.arange(-50, 10, 0.1)

plt.figure(figsize=(10,5))
for i in range(5):
    plt.plot(x, MLf(x, i), label="alpha = "+str(i))
plt.legend()
plt.ylim(-5, 5); plt.xlim(-55, 15); plt.grid()

您可以在橙色线中看到最不稳定的地方，从大约a = 1开始，但是x = -35也有问题（蓝色线）。改变要求和的项数（即a = 0）会改变发生不稳定性的j。
发生什么事？我怎样才能避免这个？

最佳答案

如果a=0，则仅当| z |案例a=1
这个函数是z**k的，从技术上讲，它是由所有z的幂级数exp(z)表示的。但是对于大的负z，这个幂级数经历了灾难性的loss of significance：单个项是巨大的，例如，z**k / k!大于1e16，但是它们的和很小（(-40)**40/factorial(40)几乎为零）。由于1e16接近双精度的极限，因此输出由截断/舍入操作的噪声控制。
一般来说，从效率和精度的角度来看，通过添加exp(-40)来计算多项式并不是最好的方法。Horner的方案已经在NumPy中实现，使用它可以简化代码：

k = np.arange(100)
return np.polynomial.polynomial.polyval(z, 1/gamma(a*k + 1))

但是，这不会为exp（z）保存序列，它的数值问题超出了NumPy的范围。
您可以使用c(k) * z**k进行计算，以获得精度（mpmath支持任意高精度的浮点运算）和降低速度（无编译代码，无矢量化）。
或者当a=1时，可以从MLf返回mpmath。
例0级数收敛，但又会带来灾难性的精度损失；现在没有可依赖的明确公式。前面提到的exp(z)是一个选项：设置非常高的精度（mpmath），并希望它足以对序列求和。另一种方法是寻找另一种计算函数的方法。
环顾四周，我发现报纸
鲁道夫·戈伦弗洛、朱丽亚·卢奇科和尤里·卢奇科。我从中取了公式（23），用于负z和0

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma
from scipy.integrate import quad

def MLf(z, a):
    """Mittag-Leffler function
    """
    z = np.atleast_1d(z)
    if a == 0:
        return 1/(1 - z)
    elif a == 1:
        return np.exp(z)
    elif a > 1 or all(z > 0):
        k = np.arange(100)
        return np.polynomial.polynomial.polyval(z, 1/gamma(a*k + 1))

    # a helper for tricky case, from Gorenflo, Loutchko & Luchko
    def _MLf(z, a):
        if z < 0:
            f = lambda x: (np.exp(-x*(-z)**(1/a)) * x**(a-1)*np.sin(np.pi*a)
                          / (x**(2*a) + 2*x**a*np.cos(np.pi*a) + 1))
            return 1/np.pi * quad(f, 0, np.inf)[0]
        elif z == 0:
            return 1
        else:
            return MLf(z, a)
    return np.vectorize(_MLf)(z, a)


x = np.arange(-50, 10, 0.1)

plt.figure(figsize=(10,5))
for i in range(1, 5):
    plt.plot(x, MLf(x, i/3), label="alpha = "+str(i/3))
plt.legend()
plt.ylim(-5, 5); plt.xlim(-55, 15); plt.grid()

这里没有数字问题。
"Computation of the Mittag-Leffler function and its derivative"

关于python - 使用NumPy的Mittag-Leffler函数的不稳定性，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/48645381/