LevinRule and SphericalBessels

Question

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.

score 9 · Accepted Answer · edited Jun 04 '13 at 18:57

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:

In[199]:= 
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:

In[200]:= 
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)];

In[201]:= 
NIntegrate`LevinIntegrandReduce[
  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.

In[217]:= 
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.

that was exactly the kind of thing I was after, thanks a lot.. — fpghost, Oct 18 '12 at 21:25

score 4 · Answer 2 · edited Oct 18 '12 at 12:41

4

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:

In[356]:= 
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

edited Oct 18 '12 at 12:41

J. M.'s missing motivation

124,525
11
401
574

answered Oct 18 '12 at 12:37

Sasha

7,373
36
47

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. (http://reference.wolfram.com/legacy/v5_2/Built-inFunctions/NumericalComputation/Integration/NIntegrate.html), 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
Also I'm particularly interested in very high $\omega$ maybe 100, and in reality have a limited WorkinPrecision in my integrand of maybe 10. In which case even turning Ag off->infty, and giving a low precisongoal still gen slwcon – fpghost Oct 18 '12 at 13:22

LevinRule and SphericalBessels

2 Answers2

Linked