# \[RFC/ANN\] OscillatoryIntegralsODE.jl: Levin method + OrdinaryDiffEq

**URL:** <https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601>\
**Category:** Package Announcements\
**Tags:** package, announcement, integral, oscillatory\
**Created:** [February 19, 2021, 5:29am UTC](https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601 "2021-02-19T05:29:14Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [February 19, 2021, 5:29am UTC](https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601/1 "2021-02-19T05:29:14Z")

</div>

I’m requesting feedback on [OscillatoryIntegralsODE.jl](https://github.com/xzackli/OscillatoryIntegralsODE.jl), which numerically evaluates integrals

I = \int\_a^b f(x) S(rx) \, dx

where f(x) is smooth and not oscillatory, but S(x) is highly oscillatory for large frequencies r. Examples of oscillatory functions include Bessel functions J\_{\nu}(rx), spherical Bessel j\_{\nu}(rx), \cos(\omega x), and the harmonic transform e^{i \omega x}. Gaussian quadrature is inefficient for these integrals, since the number of nodes required for accurate integration scales with frequency.

Two popular methods for highly-oscillatory quadrature are Filon and Levin-type integration (see [Huybrechs and Olver 2012](https://people.cs.kuleuven.be/~daan.huybrechs/research/HOQ.pdf)), both relying on polynomial interpolation in different ways. I believe @dlfivefifty has had some influential papers in this field, and wrote his thesis about such integrals. Within the Julia package ecosystem, Sheehan Olver and Mikael Slevinsky (the real experts) have worked on [ApproxFun approaches](https://github.com/JuliaApproximation/OscillatoryIntegrals.jl) for the harmonic transform.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/5/a578dd647035436028a0dc7fec545ae7ae1803cb.png)  
_It only gets worse. My average integrand contains frequencies 50-100x higher than this._

I need to efficiently compute thousands of these integrals, where S is typically a spherical Bessel function and f contains an interpolation from an ODE solution. This package attempts to deal with this using a Levin method combined with OrdinaryDiffEq.jl. The basic idea ([Levin 1994](https://www.sciencedirect.com/science/article/pii/0377042794001189)) is that you can transform something gross, like the integral of the function in the figure, into something easier, like a non-oscillatory bivariate ODE

\mathbf{p}^\prime + \mathbf{A}^T \mathbf{p} = \mathbf{f}.

Here \mathbf{f} contains the function you want to integrate, and \mathbf{A} is a time-dependent matrix. Under specific conditions, the integral is easily expressed in terms of the dot product of the ODE solution \mathbf{p} and the oscillatory functions, evaluated at the endpoints of the integration interval.

The literature typically approximates the solution of this ODE using collocation on a polynomial basis, which reduces the problem even further to a linear solve per subinterval. In my package, I just pass this to an ODE solver – I think this is equivalent to collocation, but I’m just making @ChrisRackauckas et al do my matrix factorizations. Every collocation scheme corresponds to a Runge-Kutta method, so something like `Vern9` should _basically_ be the 9-point collocation given in Levin’s original paper. In this case, the subinterval sizes are determined by the ODE timestep.

Here’s an example. The package is currently unregistered, so you’ll need to get it from GitHub.

```julia
using OscillatoryIntegralsODE
using OrdinaryDiffEq
f(x) = exp(-x^2 / 16) # we integrate the Bessel with this f
# set up problem with output type Float64, Bessel order nu=100, frequency r=100
bi = BesselJIntegral{Float64}(100., 100.) # nu, r
# now integrate over (a, b) = (1, 5). kwargs are passed to the ODE solve
levintegrate(bi, f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)

```

This integrates the function in the figure over the interval (a,b) = (1,5), and for reference, Mathematica gets `0.006311630277650583`. For nice f, a nice property is that the local tolerances used in the ODE solve propagate linearly to the global error on the integral.

```julia
@btime levintegrate($bi, $f, 1.0, 5.0, Vern9(); abstol=1e-6, reltol=1e-6)
    120.600 μs (19 allocations: 8.22 KiB)

```

Anyway, I’m hoping for both numerical analysis and code commentary! Is there a better scheme out there for adaptive oscillatory quadrature?

**Remark 1.** Levin mentions in his paper that collocation at the Chebyshev points would be good for reducing the Runge phenomenon. I think this should be something like the high stability RK methods? Do any of them work with StaticArrays?

**Remark 2.** 0 cannot be in the interval [a,b] for Bessel functions when using a Levin method, as A(x) diverges at x=0. If you need to reach 0 like me, you can use Gaussian quadrature from 0 to the first peak of the Bessel function. There are plenty of asymptotics for finding it, I have one in the docs.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 19, 2021, 5:49am UTC](https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601/2 "2021-02-19T05:49:28Z")

</div>

One scheme that seems potentially interesting is to replace

```julia
I &= \int_a^b f(x) S(rx) \, dx

```

with

```julia
I = .5\int_a^b (f(x) S(rx) +f(x+p) S(r(x+p)) ) \, dx + \int_a^{a+p} (f(x) S(rx)) \, dx - \int_b^{b+p} (f(x) S(rx)) \, dx 

```

where `p` is 1/2 the period of `S`. This will cancel out most of the high frequency oscillation, so should make the problem much better behaved.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 19, 2021, 7:46am UTC](https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601/3 "2021-02-19T07:46:12Z")

</div>

Cool stuff!

I would recommend checking out the work by Pawel Keller

[http://matwbn.icm.edu.pl/ksiazki/zm/zm25/zm2535.pdf](http://matwbn.icm.edu.pl/ksiazki/zm/zm25/zm2535.pdf)

> **[A method for indefinite integration of oscillatory and singular functions -...](https://link.springer.com/article/10.1007/s11075-007-9134-y)**
>
> We propose a general method for computing indefinite integrals of the form $$I(y) = {\\int\_0^y {g(t)k} }(t)dt\\;\\;\\;\\;(0 \\leqslant y \\leqslant y\_{{{\\text{max}}}} ),$$ where g is a smooth function, and k is a function that contains a singular factor or...

These are like Levins method but using Chebyshev polynomials to get sparse and stable systems. That yes improve with accuracy in the high frequency regime.

For the exponential oscillatory (Fourier integral) case I did a quick implementation a while ago:

> **[GitHub - JuliaApproximation/OscillatoryIntegrals.jl: Calculate oscillatory...](https://github.com/JuliaApproximation/OscillatoryIntegrals.jl)**
>
> Calculate oscillatory integrals using Julia. Contribute to JuliaApproximation/OscillatoryIntegrals.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [February 20, 2021, 3:12am UTC](https://discourse.julialang.org/t/rfc-ann-oscillatoryintegralsode-jl-levin-method-ordinarydiffeq/55601/4 "2021-02-20T03:12:21Z")

</div>

Wow, that 2007 paper by Keller seems really great. I was totally set on trying to avoid computing any Chebyshev coefficients, but reading this is making me change my mind…