Synopsis

#include <boost/math/quadrature/ooura_fourier_integrals.hpp>
namespace boost { namespace math { namespace quadrature {

template<class Real>
class ooura_fourier_sin {
public:
    ooura_fourier_sin(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real));

    template<class F>
    std::pair<Real, Real> integrate(F const & f, Real omega);

};


template<class Real>
class ooura_fourier_cos {
public:
    ooura_fourier_cos(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real))

    template<class F>
    std::pair<Real, Real> integrate(F const & f, Real omega);
};


}}}

Ooura's method for Fourier integrals computes

∫₀^∞ f(t)sin(ω t) dt

and

∫₀^∞ f(t)cos(ω t) dt

by a double exponentially decaying transformation. These integrals arise when computing continuous Fourier transform of odd and even functions, respectively. Oscillatory integrals are known to cause trouble for standard quadrature methods, so these routines are provided to cope with the most common oscillatory use case.

The basic usage is shown below:

using boost::math::quadrature::ooura_fourier_sin;
auto integrator = ooura_fourier_sin<double>();
auto f = [](double x) { return 1/x; };
double omega = 1;
auto [Is, relative_error] = integrator.integrate(f, omega);

The computed value should be π/2. Note that this integrator is more insistent about examining the error estimate, than (say) tanh-sinh, which just returns the value of the integral.

A classical cosine transform is presented below:

using boost::math::quadrature::ooura_fourier_cos;
auto integrator = ooura_fourier_cos<double>();
auto f = [](double x) { return 1/(x*x+1); };
double omega = 1;
auto [Ic, relative_error] = integrator.integrate(f, omega);

The value of Ic should be π/2e.

The integrator precomputes nodes and weights, and hence can be reused for many different frequencies with good efficiency. The integrator is pimpl'd and hence can be shared between threads without a memcpy of the nodes and weights.

Ooura and Mori's paper identifies criteria for rapid convergence based on the position of the poles of the integrand in the complex plane. If these poles are too close to the real axis the convergence is slow. It is not trivial to predict the convergence rate a priori, so if you are interested in figuring out if the convergence is rapid compile with -DBOOST_MATH_INSTRUMENT_OOURA and some amount of printing will give you a good idea of how well this method is performing.

References