I'm working on computing the inverse Fourier transform of a function derived from a matrix construction involving a Poisson-weighted system. Here's the relevant setup.The main object is a matrix-valued function matrixFunction[ω, β, τ], whose determinant is a polynomial in ω of degree ranging from 4 to 11.
I define expr[ω] as the diagonal of the matrix product MplusInv[ω] . matrixK . Transpose[MminusInv[ω]], where M±Inv[ω] are the inverses of matrixFunction evaluated at ±ω.
I'm trying to compute the inverse Fourier transform of expr[ω] using NInverseFourierTransform in Mathematica v.14.
The problem is that expr[ω] has poles at ω = 0, due to the determinant of matrixFunction vanishing there. This leads to singularities in the integrand, making the inverse Fourier transform numerically divergent. The full expression can have a denominator that's a high-degree polynomial (from 8 up to degree 22), and the singularity at the origin causes the numerical integration to fail.
pki[v_] := PDF[PoissonDistribution[1], v]; kAvg = 1; k2OverK = 2; be = 0.145; ta = 10; i = 1; matrixFunction[omega_, beta_, tau_] := { {-omega^2 + I*omega + beta*pki[1]/tau, beta*pki[2]/tau, beta*pki[3]/tau, beta*pki[4]/tau, beta*pki[5]/tau, -beta*1}, {beta*pki[1]/tau, -omega^2 + I*omega + beta*pki[2]/tau, beta*pki[3]/tau, beta*pki[4]/tau, beta*pki[5]/tau, -beta*2}, {beta*pki[1]/tau, beta*pki[2]/tau, -omega^2 + I*omega + beta*pki[3]/tau, beta*pki[4]/tau, beta*pki[5]/tau, -beta*3}, {beta*pki[1]/tau, beta*pki[2]/tau, beta*pki[3]/tau, -omega^2 + I*omega + beta*pki[4]/tau, beta*pki[5]/tau, -beta*4}, {beta*pki[1]/tau, beta*pki[2]/tau, beta*pki[3]/tau, beta*pki[4]/tau, -omega^2 + I*omega + beta*pki[5]/tau, -beta*5}, {beta*pki[1]/tau, beta*pki[2]/tau, beta*pki[3]/tau, beta*pki[4]/tau, beta*pki[5]/tau, I*omega + (1 - beta*k2OverK)}}; matrixK := { {1, 0, 0, 0, 0, 1*pki[1]/kAvg}, {0, 1, 0, 0, 0, 2*pki[2]/kAvg}, {0, 0, 1, 0, 0, 3*pki[3]/kAvg}, {0, 0, 0, 1, 0, 4*pki[4]/kAvg}, {0, 0, 0, 0, 1, 5*pki[5]/kAvg}, {1*pki[1]/kAvg, 2*pki[2]/kAvg, 3*pki[3]/kAvg, 4*pki[4]/kAvg, 5*pki[5]/kAvg, 1}}; Mminus[omega_] := matrixFunction[-omega, be, ta]; Mplus[omega_] := matrixFunction[omega, be, ta]; MminusInv[omega_] := Inverse[Mminus[omega]]; MplusInv[omega_] := Inverse[Mplus[omega]]; Mresult[omega_] := MplusInv[omega] . matrixK . Transpose[MminusInv[omega]]; expr[omega_] := Mresult[omega][[i, i]]; << FourierSeries` ft[t_?NumericQ] = NInverseFourierTransform[(\expr[\[Omega]], \[Omega], t]; Plot[Re[ft[t]], {t, -100, 100}] ``` 
InverseFourierTransform[1/omega^2, omega, t] (*-Sqrt[(\[Pi]/2)] t Sign[t]*)$\endgroup$