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 =


  • $\begingroup$ What is $x$ in your integral? $\endgroup$
    – Pedro
    Jan 23, 2013 at 12:14
  • $\begingroup$ Any positive, real number. I just updated my post. $\endgroup$
    – torbonde
    Jan 23, 2013 at 12:22
  • $\begingroup$ 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. $\endgroup$ Jan 23, 2013 at 15:24
  • $\begingroup$ I will make a comment later today with some code and errors. Or tomorrow. $\endgroup$
    – torbonde
    Jan 24, 2013 at 10:52
  • $\begingroup$ Okay, so I forgot. But now I updated my post with an example (I have split the integral in two by calculating $sinc'$ explicitly). $\endgroup$
    – torbonde
    Jan 28, 2013 at 13:01

4 Answers 4


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 =
err =
npoints =

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.

  • $\begingroup$ 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. $\endgroup$
    – torbonde
    Jan 23, 2013 at 13:24
  • $\begingroup$ @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. $\endgroup$
    – Pedro
    Jan 23, 2013 at 13:30
  • $\begingroup$ I maybe don't need ten digits of accuracy, but I think I need at least five. Can your method produce that? $\endgroup$
    – torbonde
    Jan 28, 2013 at 13:34
  • $\begingroup$ 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. $\endgroup$
    – Pedro
    Jan 28, 2013 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}]]

Plot from x = 0.5 to x=10

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` + 
     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.

  • $\begingroup$ 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. $\endgroup$
    – torbonde
    Jan 28, 2013 at 13:16
  • $\begingroup$ OK, as you say you will have to make do with Matlab. I will add another answer about that. $\endgroup$ Jan 28, 2013 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.

  • $\begingroup$ 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. $\endgroup$
    – torbonde
    Jan 29, 2013 at 7:58
  • $\begingroup$ 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. $\endgroup$
    – torbonde
    Jan 29, 2013 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:
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)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.