任何获得虚误差函数的数值稳定差的方法

如何解决任何获得虚误差函数的数值稳定差的方法

我正在尝试使用 numpy 和 scipy 在 Python 中绘制半高斯的傅立叶逆变换，我不感兴趣。该函数可以解析导出，我的实现如下所示，

import numpy as np
from scipy.special import erfi

def one_sided_gauss_fourier(freq,T,a):
    ift = 0.5j * np.sqrt(np.pi/a) * np.exp(-freq**2/(4.0*a)) * ( 
        erfi(freq/(2*np.sqrt(a))) - erfi(freq+T*2j*a/(2*np.sqrt(a)))
        )
    return ift/T

这是形式函数的相应逆傅立叶变换

def gauss(a,t):
    return np.exp(-a*t**2)

t = np.linspace(0,4,10**3)
a = 1.0
dft_gauss_ifft = np.fft.ifft(gauss(a,t))
# Should be the same as 
dt = t[1] - t[0]
freq = np.fft.fftfreq(len(t),dt) * np.pi * 2
analytic_gauss_ifft = one_sided_gauss_fourier(freq,t.max(),a)

有问题的部分是 erfi(freq/(2*np.sqrt(a))) - erfi(freq+T*2j*a/(2*np.sqrt(a))) 包含角频率轴的数组 freq 可以具有 >> 2*np.sqrt(a) 的值，这会导致对带有参数 erfi 的 >>1 函数求值，从而导致溢出错误。它适用于 freq/(2*np.sqrt(a)) 不太大的值，但对于该范围之外的频率，我得到了错误，尽管结果应该是一个幅度小于 1 的数字。有问题的部分是导致加法的中间 erfi 评估/np.inf 项的减法。

是否有可能以稳定的方式评估这两个假想误差函数的差异，无论是数值上还是通过进一步操作方程来消除差异？