JuliaMath / Cubature.jl

One- and multi-dimensional adaptive integration routines for the Julia language
Other
126 stars 21 forks source link

pquadrature for oscillating integrals #16

Closed mzaffalon closed 7 years ago

mzaffalon commented 8 years ago

I get wrong results when I use pquadrature to compute the Fourier transform of a Gaussian, but only when abstol is set to a value other than the default value.

The commented out lines (changing the integration limits, leaving abstol to its default value and using quadgk) all give the correct result.

using PyPlot
using Cubature

function test_numerical_integration_FT()
    t_c, σ = 0.05, .7
    freq = linspace(-.2, .2, 21)
    y_imag = zeros(freq)
    for i = 1:length(freq)
        f_red = 2π * freq[i]
        y_imag[i] = -pquadrature(t -> exp(-(t - t_c)^2 / 2σ^2) * sin(f_red * t), -5., 5., abstol=1e-10)[1]
        #y_imag[i] = -pquadrature(t -> exp(-(t - t_c)^2 / 2σ^2) * sin(f_red * t), -5., 5.1, abstol=1e-10)[1]
        #y_imag[i] = -pquadrature(t -> exp(-(t - t_c)^2 / 2σ^2) * sin(f_red * t), -5., 5.)[1]
        #y_imag[i] = -quadgk(t -> exp(-(t - t_c)^2 / 2σ^2) * sin(f_red * t), -5., 5., abstol=1e-3, reltol=1e-3)[1]
    end
    # analytical FT
    y = sqrt(2π)σ * exp(-2π * freq .* (π * freq * σ^2 + 1im * t_c))
    figure(); grid(true)
    plot(freq, y_imag, label="numerical, imag")
    plot(freq, imag(y), label="analytical, imag")
    legend(frameon=false)
end