Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm currently looking at a simplified problem that approximates another problem I'm looking into. In this simplified problem I at least have an analytic integrand and can easily provide all info on here:

Given the definitions

v = 0.6;
g = 1/Sqrt[1 - v^2];

ig[tau1_?NumericQ, tau2_?NumericQ, ωω_?NumericQ, ll_?IntegerQ] :=
(((2 ll + 1) ωω)/(4 Pi^2)) *
 SphericalBesselJ[ll, ωω v g tau1] *
 SphericalBesselJ[ll, ωω v g tau2] *
 Exp[-I ωω g (tau1 - tau2)];

I would like to numerically integrate this:

f[ω_, l_, pg_, wp_] := 
 2 Re[NIntegrate[
    Exp[-I 1 s] ig[100, 100 - s, ω, l], {s, 0, 40}, 
    PrecisionGoal -> pg, WorkingPrecision -> wp, MaxRecursion -> 20]]

But for example f[2, 1, 15, 25] etc I get a host of errors like ::slwcon,::einc.

I was wondering if I could use "LevinRule" here and if so what would the options be?

I believe the problem gets worse at large $\omega$ maybe 100 or more. Where even putting the WorkingPrecision upto 40 and PrecisionGoal down to 10, things still scream, if I generate a table of these from 1 to 100 in omega.

share|improve this question
up vote 8 down vote accepted

There are two oscillatory parts to your integrand: Exp[-I 1 s] and ig[100, 100 - s, ω, l].

However, because you have defined the function ig for strictly numerical values of its arguments, NIntegrate is not able to "see" its symbolic structure and realize that it is oscillatory with a certain form. Therefore, "LevinRule" can only automatically be applied to the Exp[-I 1 s] part.

You can see what "LevinRule" does:

NIntegrate`LevinIntegrandReduce[Exp[-I 1 s] ig[100, 100 - s, 2, 1], s]["Rules"]

Out[199]= {"Variables" -> {s}, "AdditiveTerm" -> 0, 
 "Amplitude" -> {ig[100, 100 - s, 2, 1]}, "Kernel" -> {E^(-I s)}, 
 "DifferentialMatrices" -> {{{-I}}}}

For more information on the LevinIntegrandReduce function, please see this part of the documentation.

However if you define your ig function to allow symbolic evaluation, then NIntegrate can automatically detect the other oscillatory part of your integrand:

igsymbolic[tau1_, tau2_, ωω_, 
   ll_] := (((2 ll + 1) ωω)/(4 Pi^2))*
   SphericalBesselJ[ll, ωω v g tau1]*
   SphericalBesselJ[ll, ωω v g tau2]*
   Exp[-I ωω g (tau1 - tau2)];

  Exp[-I 1 s] igsymbolic[100, 100 - s, 2, 1], s]["Rules"]

Out[201]= {"Variables" -> {s}, "AdditiveTerm" -> 0, 
 "Amplitude" -> {(3 SphericalBesselJ[1, 150])/(2 π^2), 0}, 
 "Kernel" -> {E^(-((7 I s)/2)) SphericalBesselJ[1, (3 (100 - s))/2], 
   E^(-((7 I s)/
     2)) (-(SphericalBesselJ[1, (3 (100 - s))/2]/(3 (100 - s))) + 
      1/2 (SphericalBesselJ[0, (3 (100 - s))/2] - 
         SphericalBesselJ[2, (3 (100 - s))/2]))}, 
 "DifferentialMatrices" -> {{{-((7 I)/2), -(3/2)}, {(
     2 (-2 + 9/4 (100 - s)^2))/(
     3 (100 - s)^2), -((7 I)/2) + 2/(100 - s)}}}}

This results in a more efficient integration.

fsymbolic[ω_, l_, {wp_, pg_, ag_}, mthd_] := 
 2 Re[NIntegrate[
    Exp[-I 1 s] igsymbolic[100, 100 - s, ω, l], {s, 0, 40}, 
    WorkingPrecision -> wp, AccuracyGoal -> ag, PrecisionGoal -> pg, 
    Method -> mthd]]

In[218]:= fsymbolic[2, 1, {25, 15, 18}, Automatic] // Timing

Out[218]= {0.343, -1.925723011031919293641188*10^-6}

Scaling up to large values of l is also fast:

In[219]:= fsymbolic[2, 100, {25, 15, 18}, Automatic] // Timing

Out[219]= {0.608, -0.00002319232924789404095060205}

Note that Method -> Automatic is fine here since NIntegrate automatically detects the oscillatory part of the integrand and selects LevinRule.

share|improve this answer
that was exactly the kind of thing I was after, thanks a lot.. – fpghost Oct 18 '12 at 21:25

The problem is that the numerical value of the integral is small. Thus in order to appease NIntegrate's error estimator you have to lower the precision/accuracy goals.

Then both LevinRule and Automatic methods work well:

f[ω_, l_, {wp_, pg_, ag_}, mthd_] := 
 2 Re[NIntegrate[
    Exp[-I 1 s] ig[100, 100 - s, ω, l], {s, 0, 40}, 
    WorkingPrecision -> wp, AccuracyGoal -> ag, 
    PrecisionGoal -> pg, Method -> mthd]]

In[364]:= f[2, 1, {25, 15, 18}, "LevinRule"]

Out[364]= -1.925723011031859693480016*10^-6

In[366]:= f[2, 1, {25, 15, 18}, Automatic]

Out[366]= -1.925723011031860585537216*10^-6
share|improve this answer
Out of interest I would still quite like to know how to specify the LevinRule options for this integrand, i.e. what one might put for "Kernal" and that kind of thing. – fpghost Oct 18 '12 at 12:48
You just use Method->"LevinRule" as an option of NIntegrate. – Sasha Oct 18 '12 at 13:00
Sasha: I know you can just do that. I'm wondering about refining the options of LevinRule to better suit the integrand if possible.For example such options as ExcludeForms, AdditiveTerms, Kernal.. – fpghost Oct 18 '12 at 13:03
It's interesting to note (for me at least) that the default behaviour (i.e. omitting) AccuracyGoal is to set to Infinity, whereas "Automatic" means processor decides. Whereas PrecisionGoal has the opposite behaviour the default (i.e omitting it from options) is to set it to Auto decide, whereas setting it to infinity would turn it off. – fpghost Oct 18 '12 at 13:14
Thus if you try f[2,1,{25,15,Infinity},Automatic] you set the ag to Infinity (turn it off and use only pg) this gives slwcon. But if you put in f[2,1,{25,15,Automatic},Automatic] you don't get slwcon any more. (…), I guess it chooses a low enough value of Ag as it likes with this option. I had originally expected ag=Auto to reproduce it's absence... – fpghost Oct 18 '12 at 13:17

Your Answer


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.