Method for numerical integration of difficult oscillatory integral

Question

I need to numerically evaluate the integral below:

$$\int_0^\infty \mathrm{sinc}'(xr) r \sqrt{E(r)} dr$$

where $E(r) = r^4 (\lambda\sqrt{\kappa^2+r^2})^{-\nu-5/2} K_{-\nu-5/2}(\lambda\sqrt{\kappa^2+r^2})$, $x \in \mathbb{R}_+$ and $\lambda, \kappa, \nu >0$. Here $K$ is the modified Bessel function of the second kind. In my particular case I have $\lambda = 0.00313$, $\kappa = 0.00825$ and $\nu = 0.33$.

I am using MATLAB, and I have tried the built-in functions integral and quadgk, which gives me a lot of errors (see below). I have naturally tried numerous other things as well, such as integrating by parts, and summing integrals from $kx\pi$ to $(k+1)x\pi$.

So, do you have any suggestions as to which method I should try next?

UPDATE (added questions)
I read the paper @Pedro linked to, and I don't think it was too hard to understand. However, I have a few questions:

• Would it be okay to use $x^k$ as the basis-elements $\psi_k$, in the univariate Levin method described?
• Could I instead just use a Filon method, since the frequency of the oscillations is fixed?

Example code
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89

ans =

3.3197e+06

Answer 1

I've written my own integrator, quadcc, which copes substantially better than the Matlab integrators with singularities, and provides a more reliable error estimate.

To use it for your problem, I did the following:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

The function f is now your integrand. Note that I've just assigned any old value to x.

In order to integrate on an infinite domain, I apply a substitution of variables:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

i.e. integrating g from 0 to 1 should be the same as integrating f from 0 to $\infty$. Different transforms may produce different quality results: Mathematically all transforms should give the same result, but different transforms may produce smoother, or more easily integrable gs.

I then call my own integrator, quadcc, which can deal with the NaNs on both ends:

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

Note that the error estimate is huge, i.e. quadcc doesn't have much confidence in the results. Looking at the function, though, this is not surprising as it oscillates at values three orders of magnitude above the actual integral. Again, using a different interval transform may produce better results.

You may also want to look at more specific methods such as this. It's a bit more involved, but definitely the right method for this type of problem.

Answer 2

As Pedro points out, Levin-type methods are the best established methods for these kinds of problems.

Do you have access to Mathematica? For this problem, Mathematica will detect and use them by default:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

Here is a plot over a range of values of x:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

You can also manually specify the particular Levin-type method to apply, which in this case can yield a slight performance improvement:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

See the documentation for details of Levin-type methods in Mathematica.

Answer 3

If you don't have access to Mathematica, you could write a Levin-type (or other specialized oscillatory) method in Matlab as Pedro suggests.

Do you use the chebfun library for Matlab? I just learned it contains an implementation of a basic Levin-type method here. The implementation is written by Olver (one of the experts in the oscillatory quadrature field). It doesn't deal with singularities, adaptive subdivision etc, but it may be just what you need to get started.

Answer 4

The transformation recommended by Pedro is a great idea. Have you tried to play around with the parameters in Matlab's "quadgk" function? For example, using Pedro's transformation, you can do the following:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Using this gives me a solution of:
-2184689.50220729
and only takes 0.8 seconds (using the values mentioned above: x=10)
Walter Gander and Walter Gautschi have a paper on adaptive quadrature with Matlab code you can use as well (link here)