Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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

share|improve this question
    
What is $x$ in your integral? –  Pedro Jan 23 '13 at 12:14
    
Any positive, real number. I just updated my post. –  cimrg.joe 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. –  cimrg.joe 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). –  cimrg.joe Jan 28 '13 at 13:01
show 2 more comments

4 Answers

up vote 12 down vote accepted

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.

share|improve this answer
    
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. –  cimrg.joe 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? –  cimrg.joe 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
add comment

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

share|improve this answer
    
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. –  cimrg.joe 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
add comment

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.

share|improve this answer
    
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. –  cimrg.joe 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. –  cimrg.joe Jan 29 '13 at 8:02
add comment

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)

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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