# Method for numerical integration of difficult oscillatory integral

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?

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

• What is $x$ in your integral? – Pedro Jan 23 '13 at 12:14
• Any positive, real number. I just updated my post. – torbonde Jan 23 '13 at 12:22
• If you could show some code and errors it is probably not too hard to solve most of them. Of course please try reading the error carefully first and see whether you can make it dissapear on your own. – Dennis Jaheruddin Jan 23 '13 at 15:24
• I will make a comment later today with some code and errors. Or tomorrow. – torbonde Jan 24 '13 at 10:52
• Okay, so I forgot. But now I updated my post with an example (I have split the integral in two by calculating $sinc'$ explicitly). – torbonde Jan 28 '13 at 13:01

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.

• Thank you very much. I will have a look at the different methods. For my purposes, the error does not need to be as small as is standard in e.q. integral (1e-10 I think), but 1.7e+07 is still really, really large. Perhaps another transform will do good, as you mention. – torbonde Jan 23 '13 at 13:24
• @cimrg.joe: Note that the error estimate is an estimate of the absolute error based, amongst others, on the maximum absolute values of the integrand. In some extreme cases, the returned value may actually be quite ok. If you're looking for ten digits of accuracy, then I strongly recommend using the Levin-type methods I mentioned at the end of my post. – Pedro Jan 23 '13 at 13:30
• I maybe don't need ten digits of accuracy, but I think I need at least five. Can your method produce that? – torbonde Jan 28 '13 at 13:34
• The method can't guarantee that kind of precision for your integral since the values at the right end of the interval are several orders of magnitude larger than the integral itself. – Pedro Jan 28 '13 at 14:59

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 +

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.

• Unfortunately I don't have access to Mathematica - only MATLAB. I will just update my question with some added questions, regarding the paper @Pedro linked to. – torbonde Jan 28 '13 at 13:16
• OK, as you say you will have to make do with Matlab. I will add another answer about that. – Andrew Moylan Jan 28 '13 at 16:25

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.

• I have thought about implementing a Levin method myself, but I'm not sure if I'm up for the challenge yet. I think I need to understand the method a little better. Maybe I could talk to my advisor about that. Anyway, the reason that I asked about the Filon methods, is that they appear easier to implement. And since I don't need extremely high accuracy, but this is part of my master thesis, difficulty weighs in. – torbonde Jan 29 '13 at 7:58
• I have had a look at the chebfun library (which is impressive) and the Levin integration-example. But I can't get it running. I have actually posted a question regarding it here. – torbonde Jan 29 '13 at 8:02

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)