python - 来自 scipy.special 的 fadeeva 函数的二阶导数

Question

我想计算 Fadeeva 函数 special.wofz 的二阶导数。 Fadeeva 函数与错误函数密切相关。因此，如果有人更熟悉 erf，我们将不胜感激。 下面是求 wofz 二阶导数的代码:

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


def Z(x):
    return wofz(x)

## first derivative of wofz (analytically)
def Zp(x):
    return -2/1j/np.pi**0.5 - 2*x*Z(x)

##second derivative (analytically)
def Zpp(x):
    return (Z(x)+x*Zp(x))*x

x = np.float64(np.linspace(1e4,14e4,1000))

plt.plot(x, Zpp(x).imag,"-")

Zpp_num=np.diff(Zp(x))/np.diff(x)  ##calc numerically the second derivative
plt.plot(x[:-1],Zpp_num.imag)

代码生成下图:

显然，解析计算存在严重错误。我一直使用的公式是正确的。这一定是一些数字问题。

问:有人能告诉我是什么原因导致了这种行为吗？是由于 wofz 函数的精度吗？有谁知道计算 wofz 的算法？产生可靠结果的争论有多大？我找不到任何关于它的信息。另外，我知道我可以使用 wofz 的渐近逼近来找到二阶导数，但如果可能的话，我想使用 scipy。

Answer 1

正如您所怀疑的，问题出在计算导数时的数值起源上。正确的二阶导数，如 @clwainwright已经在评论中指出，是

Zpp = -2*(Z(x) + x*Zp(x))

这两项的虚部表现如下:

显示您有两个几乎相等的小量，然后计算它们的差值。

进入更多细节，

Zpp = -2*(1-2*x**2)*Z(x) - 4j/sqrt(pi)*x

其中的虚部是

Im(Zpp) = - 4*x/sqrt(pi) - 2*(1-2*x**2)*Im(Z)

和Im(Z)正比于道森函数D (scipy.special.dawsn),

Im(Z) = 2/sqrt(pi) * D

问题是你有

Im(Zpp(x)) = -4/sqrt(pi)*( x - 2*x**2*dawsn(x) ) - 4/sqrt(pi)*dawsn(x)

为什么这是一个问题是因为the asymptotic expansion of the Dawson function开始于

D(x) ~ 1/(2x) + ...

它的前导项被 Im(Zpp(x)) 的第一项吃掉了，小的修正赋予函数它的值(实际上，前导项是 1/(2x) 在 Im(Zpp(x)) 的最后一项。

所以问题出在Zpp的解析表达式上。您可以尝试 reshape 解析表达式以摆脱这个数值问题(特别是精度损失)，但这并不容易。您也可以尝试使用 sympy。我已经尝试了一段时间，但没有成功。它可能仍然是可能的。