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

Question

我正试图在Python中使用numpy和scipy绘制半个高斯的傅里叶逆变换，我对全高斯不感兴趣。这个函数可以用解析的方法导出，我的实现是这样的

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)的值，这导致使用参数>>1对erfi函数求值，从而导致溢出错误。它适用于freq/(2*np.sqrt(a))不是太大的值，但是对于该范围之外的频率，我会得到错误，尽管结果应该是一个小幅度的数字<；1.有问题的部分是导致np.inf术语加减的中间erfi评估

是否有可能以稳定的方式计算这两个虚函数的差，或者通过数值计算，或者通过进一步操纵方程来消除差