# Package to evaluate oscillating integrals with increasing frequency

**URL:** <https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565>\
**Category:** General Usage\
**Tags:** quadrature\
**Created:** [January 9, 2024, 1:58pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565 "2024-01-09T13:58:54Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [January 9, 2024, 1:58pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/1 "2024-01-09T13:58:54Z")

</div>

I’m looking to numerically integrate a function of the form

`int[g(x)*sin(r(x)*x) dx, x=a...infinity]`

Where `g(x)` is exponentially decaying and the frequency `r(x) = A+B*sqrt(1+x^2)` is increasing with `x`.  
I saw [GitHub - xzackli/OscillatoryIntegralsODE.jl](https://github.com/xzackli/OscillatoryIntegralsODE.jl) and [GitHub - machakann/DoubleExponentialFormulas.jl: One-dimensional numerical integration using the double exponential formula](https://github.com/machakann/DoubleExponentialFormulas.jl), but they seem to require a constant oscillation frequency.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [January 9, 2024, 2:12pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/2 "2024-01-09T14:12:34Z")

</div>

Maybe a change of variable

(A+B\sqrt{1+x^2})x=Ay

could work out ? (can be solved for x^2 to express dx and g(x) in terms of y)

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [January 9, 2024, 2:27pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/3 "2024-01-09T14:27:08Z")

</div>

I’ll see if that works.

---

<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:** [January 9, 2024, 2:31pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/4 "2024-01-09T14:31:15Z")

</div>

If `g` decays quickly enough, you probably don’t need to do anything special. Just integrate normally until `g<tol`

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [January 9, 2024, 2:36pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/5 "2024-01-09T14:36:20Z")

</div>

Do you already know that QuadGK.jl doesn’t work? Since it uses adaptive routines, I’d think it would be able to do it unless the oscillation was very very extreme. It may be that it is not the most efficient method though since it will not take advantage of the knowledge of the integrand.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [January 9, 2024, 2:55pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/6 "2024-01-09T14:55:26Z")

</div>

I’m new to this, but my understanding was that a quadrature made for specifically oscillating functions was always much more efficient. I’ll experiment and report back.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2024, 3:01pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/7 "2024-01-09T15:01:48Z")

</div>

Have you tried just brute-forcing it with adaptive integration? For example, QuadGK.jl works just fine with a few ~~hundred~~ thousand evaluations for a simple example of the form you proposed:

```julia
julia> using QuadGK

julia> quadgk_count(x -> exp(-x) * sin((1 + 2hypot(1,x))*x), 1, Inf) # default 1e-8 tolerance
(-0.05977814307911682, 8.548339483678494e-10, 2115)

julia> quadgk_count(x -> exp(-x) * sin((1 + 2hypot(1,x))*x), 1, Inf, rtol=1e-13)
(-0.05977818481933722, 5.807105271695107e-15, 11175)

```

Of course, if you analyze your problem and are clever there are almost always ways to do integrals much more efficiently, but it’s not clear from your question whether this is actually necessary for you.

---

<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:** [January 9, 2024, 3:06pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/8 "2024-01-09T15:06:32Z")

</div>

In general, you only really need specific oscillating function integrator when your decay term is slow. If your decay term is `1/sqrt(x)` you need to do special stuff because you need to integrate really far out to capture all the mass, but with exponential decay all the mass is close to 0 anyway so you don’t need to care about the oscillations because after not very long the function is basically equal to 0.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [January 9, 2024, 6:52pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/9 "2024-01-09T18:52:47Z")

</div>

This program will need to call this function a few thousand times, so I would like it to be efficient. But I think I’ll start with the QuadGK version and only get more clever if and when required.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2024, 7:08pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/10 "2024-01-09T19:08:34Z")

</div>

> [@Oscar\_Smith](#):
>
> In general, you only really need specific oscillating function integrator when your decay term is slow.

Note also that if your decay term is slow (relative to the oscillation period), then direct methods (not specialized for the oscillation term) correspond to very ill-conditioned sums that are subject to large [cancellation errors](https://en.wikipedia.org/wiki/Catastrophic_cancellation) — basically, you have a lot of positive and negative contributions that are supposed to exactly cancel, but don’t because of roundoff and other errors.

> [@Oscar\_Smith](#):
>
> but with exponential decay all the mass is close to 0 anyway so you don’t need to care about the oscillations because after not very long the function is basically equal to 0.

Even with exponential decay, if the rate is slow enough then you could run into trouble with non-specialized integration schemes.

A classic example of this is asymptotics of Fourier transforms — even if you have a nice smooth, rapidly decaying function f(x), if you want to compute high-frequency asymptotics of its Fourier transform F(k) (which goes to zero for large k, but maybe you want to know just how small it is), then you often need analytical tricks such as [saddle-point methods](https://arxiv.org/abs/1508.04376), and/or [specialized quadrature schemes](https://doi.org/10.1093/imanum/24.3.365).

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [January 10, 2024, 12:26pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/11 "2024-01-10T12:26:22Z")

</div>

Thanks for the information.

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [February 18, 2024, 11:21am UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/12 "2024-02-18T11:21:10Z")

</div>

I’ve come back to this as it turns out the decay rate can become very slow in my application. There is a region of stationary phase, and this dominates the integral, but it’s unclear to me how to take advantage of this numerically. Any papers or packages that deal with such functions is greatly appreciated.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [February 20, 2024, 8:30pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/13 "2024-02-20T20:30:19Z")

</div>

Other than the papers mentioned above which seem intended for finite intervals (edit: the numerical support of an exponentially decaying function in floating point arithmetic is finite), another general option for specialized quadrature schemes is called Generalized Chebyshev quadrature (which is a precursor for Generalized Gaussian quadrature). The basic idea is to choose a basis of functions which efficiently captures the behavior of those you would like to integrate and that is in some sense low-rank. Then using the integrals of the basis functions (obtained either analytically or by brute force), sample the basis functions densely and find the most linearly-independent points to use for a specialized quadrature scheme. There is overhead to calculating the quadrature, which could use as many points as you need and be reused for many integrals. For badly behaved integrals, e.g. power-law singularities, I’ve seen I can get 3-5 digits of accuracy in the integral from ~40 quadrature points, often because of the aforementioned cancellation errors.

Moreover, it looks like there are some unregistered packages for exactly these calculations: [GitHub - daanhb/GeneralizedGauss.jl: Package for the computation of generalized Gaussian quadrature rules](https://github.com/daanhb/GeneralizedGauss.jl)

---

<div class="post-metadata">

**Author:** ![weymouth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/weymouth/32/15839_2.png) [@weymouth](https://discourse.julialang.org/u/weymouth)\
**Post date:** [February 21, 2024, 2:14pm UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/14 "2024-02-21T14:14:07Z")

</div>

Thanks for the information and the nod towards that package authored by @daanhb. In his paper related to that repo, he gives an example dealing with boundary element method quadratures, which is what I need. Looks like I’ll need to read up on _numerical_ steepest decent methods to finish this up…

---

<div class="post-metadata">

**Author:** ![daanhb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/daanhb/32/209040_2.png) [@daanhb](https://discourse.julialang.org/u/daanhb)\
**Post date:** [March 4, 2024, 9:34am UTC](https://discourse.julialang.org/t/package-to-evaluate-oscillating-integrals-with-increasing-frequency/108565/15 "2024-03-04T09:34:12Z")

</div>

Checking in late here, but thought I’d mention there is a book on the topic of computing oscillatory integrals: [https://epubs.siam.org/doi/book/10.1137/1.9781611975123](https://epubs.siam.org/doi/book/10.1137/1.9781611975123)

Lots of work in the area happened for highly oscillatory integrals closely related to asymptotics as mentioned above. The goal here is to evaluate an oscillatory integral at a cost that is independent of the frequency parameter involved. Not a lot of software exists, one recent exception being the PathFinder toolbox ([GitHub - AndrewGibbs/PathFinder: A Matlab/Octave package for oscillatory integration](https://github.com/AndrewGibbs/PathFinder), Matlab I’m afraid) which automates the numerical method of steepest descent for polynomial phases (i.e. for integrands of the form f(x) _exp(i_omega\*g(x)) with polynomial g).