I have an integrand that looks like this:

enter image description here

the details of computation are complicated but I only know the integrand numerically (I use NDSolve to solve second order ODE). The integrand is not simply the solution of my ODE either; calling the two solutions of my ODE osc1[s], osc2[s] then schematically the integrand I have looks something like exp(-is)[g(s)osc1[s]osc2*[C-s]+f(s)osc2[s]osc2*[C-s]]. The exp bit is only very slowly oscillating over my integration range, it is really osc1,osc2 that give wild oscillation, as a certain parameter they depend on gets larger.

More explicitely

rstar[r_] := r + 2 M Log[r/(2 M) - 1];
rH = 200001/100000;
r0 = 10;
\[Lambda][l_] = l (l + 1);

eq[\[Omega]_,l_] := \[CapitalPhi]''[r] + (2 (r - M))/(
r (r - 2 M)) \[CapitalPhi]'[
r] + ((\[Omega]^2 r^2)/(r - 2 M)^2 - \[Lambda][l]/(
r (r - 2 M))) \[CapitalPhi][r] == 0;
init=-0.0000894423075560122420468703835499 + 
0.0000447222944185058822813688948339 I;
dinit=-4.464175354293244250869336196691640386266791`30.*^-6 - 
8.950483248390306670770345406047835993931665`30.*^-6 I;

osc1 := \[CapitalPhi] /. 
Block[{$MaxExtraPrecision = 100}, 
NDSolve[{eq[1/10, 1], \[CapitalPhi][rinf] == 
init, \[CapitalPhi]'[rinf] == dinit}, \[CapitalPhi], {r, r0, 
rinf}, WorkingPrecision -> wp, AccuracyGoal -> ac, 
MaxSteps -> \[Infinity]]][[1]];

osc2 is obtained simiarly. Note these are for non problematic params and it will run quick quickly, and not be too badly behaved.

The problem I have is that I only know the integrand to maybe 6-12 digits of precision (dp), depending on the parameters. This is computing the NDSolve with a WorkingPrecision of 50-60, AccuracyGoal->42-52 and it takes around 2 hrs. I want to integrate this with NIntegrate, but when my parameters are large (and the oscillation is very high) I usually only know the integrand around the 6 dp end of scale, and NIntegrate wants a greater WorkingPrecision than this otherwise it complains (since oscillation is also getting very large).

I can force it to do the integral by making the WorkingPrecision higher, but I think this is cheating if I don't believe my integrand any higher than 6 dp?

The only ideas I've had so far are to try different rules. Are there any rules people would recommend for doing such oscillatory integrands? So far I've tried "LevinRule", "ClenshawCurtisRule", "GaussKronrodRule" but none seem to compute it any quicker than just the default. They all agree up to a reasonable number of dp, so no idea if I should just stick to the default, or if there is something better one could do with such an integrand. Speed is not a concern just accuracy.


Let's say I managed to split my integral into a few different integrals. First give the definitions:

NumericQ] := (4 M) ((tau/tauh)^(1/3) + 1) Exp[-(tau/tauh)^(1/3) + 
 1/2 (tau/tauh)^(2/3) - 1/3 (tau/tauh)];
NumericQ] := -(4 M) ((tau/tauh)^(1/3) - 1) Exp[(tau/tauh)^(1/3) + 
 1/2 (tau/tauh)^(2/3) + 1/3 (tau/tauh)];
rtau[tau_?NumericQ] := (2 M) (tau/tauh)^(2/3);

in addition to those made above, then I think I can give my integral as a sum of integrands that look like this

Exp[-I s] (ubar[tau_f - s])^(-i 4/10)Exp[+i 1/10 rstar[tau_f - s]]osc1[rtau[tauf-s]]*

here tau_f constant. The first part is an amplitude, the osc1 satisfies the linear ODE given above. I think this has Levin potential if I can work out how to input the LevinRules given the above second order ODE? (Here and in the above I fix my parameters the ODE depends on to (1/10,1) to simplify giving the ICs but I don't that detracts from the main problem). Would need to work out what the Kernal is from the ODE above.

  • $\begingroup$ What about rewriting your differential equations to directly give the integral from NDSolve? $\endgroup$
    – celtschk
    Commented Oct 16, 2012 at 17:34
  • $\begingroup$ Have you tried a Fourier on that ? $\endgroup$ Commented Oct 16, 2012 at 17:34
  • $\begingroup$ Have you seen this video ? It talks about similar problem and hybrid numeric-symbolic methods to address it. $\endgroup$ Commented Oct 16, 2012 at 17:43
  • 1
    $\begingroup$ "I can force it to do the integral by making the WorkingPrecision higher, but I think this is cheating if I don't believe my integrand any higher than 6 dp?" <-- Definitely. You will need the integrand to higher precision. $\endgroup$ Commented Oct 17, 2012 at 2:23
  • 1
    $\begingroup$ @fpghost: Hint: If you want people to be notified of your answers to their comments, prepend @ to the user name of the user you reply to (as I did with yours in this comment, although in this case it's not strictly necessary because question/answer authors always get notified about comments on their post). $\endgroup$
    – celtschk
    Commented Oct 19, 2012 at 7:55

3 Answers 3


This is how you manually invoke "LevinRule" when you know part of the integrand is a rapidly oscillatory function satisfying a linear ODE:

First, a rapidly oscillatory function:

In[25]:= osc = 
 y /. NDSolve[{y''[x] - (x^2 - 3 x) y'[x] + 10000 y[x] == 0, 
     y[0] == 3, y'[0] == 1}, y, {x, 0, 5}] // First

Out[25]= InterpolatingFunction[{{0.,5.}},<>]

In[26]:= Plot[osc[x], {x, 0, 5}]

enter image description here

The integrand is f[x]*osc[x]:

In[27]:= f[x_] := x^2

Regular NIntegrate (it can't detect that your purely numerical InterpolatingFunction has a special oscillatory form):

In[28]:= NIntegrate[f[x] osc[x], {x, 0, 5}] // Timing

Out[28]= {0.116796, -2.80375}

Manually invoke "LevinRule":

In[29]:= NIntegrate[f[x] osc[x], {x, 0, 5}, 
  Method -> {"LevinRule", "AdditiveTerm" -> 0, 
    "Amplitude" -> {f[x], 0}, "Kernel" -> {osc[x], osc'[x]}, 
    "DifferentialMatrix" -> {{0, 1}, {-10000, x^2 - 3 x}}}] // Timing

Out[29]= {0.032645, -2.80375}

Note that the first step, actually constructing an interpolation of the rapidly oscillatory function, is often the most costly in a scheme like this.

For more information about how "LevinRule" works in NIntegrate see this part of the documentation.

  • $\begingroup$ Thanks. My integrand looks like Exp[-iEs] w[s], where it's actually not the Exp[..] that causes rapid oscillation it's the w[s] bit. This w[s] doesn't directly solve an ODE but it's constructed by taking the modulus of two particular solutions of an ODE (with some other complicated factors). It's these solutions of the ODE that are highly oscillatory. Given this, would it be possible for me to form the "DifferntialMatrix" ? $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 8:08
  • $\begingroup$ Using your example it's a bit like my integrand was Exp[-is] (g(s) osc1[C]osc1*[C-s]+h(s) osc2[C]osc2*[C-s]) , where C is some constant, g(s),h(s) are just functions of integration var, and osc1, osc2 are two distinct solutions to the ODE for different ICs. $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 8:12
  • $\begingroup$ I'm finding the very strange behaviour that for a given integrand, if I stick with WorkingPrecision->5. Then the NIntegrate goes through with PrecisionGoal->3, but gives ::slwcon with PrecisionGoal->1. This is the opposite of what I would have expected, why could this happen? $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 18:13
  • $\begingroup$ I don't think my integrand osc bit satifies a differential Kernal, so don't know what to put there, but If I try Method->{"LevinRule"} alone I see no improvement in timing or ::slwcon errors going away. If I try Method->{"LevinRule","Kernal"->w[s]} where w[s] is the oscillatory bit, I get ::nonlev errors, saying Integrand not a Levin function. If I try Method->{"LevinRule","ExcludedForms"->Exp[-iEs]} then again no increase in speed and ::slwcon errors still present. Not sure what to do as I can't really increase the number of good decimals of my integrand either, as this is limited by NDSolve $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 18:21
  • $\begingroup$ If the oscillatory part of your integrand satisfies a nonlinear rather than linear ODE, then LevinRule can't be applied. $\endgroup$ Commented Oct 18, 2012 at 5:37

One might consider using the simple-minded strategy of splitting the known oscillatory part over its roots (or extrema), evaluating the integral over the intervals determined by the roots, and summing all those integrals to arrive at the actual integral you need.

Now, finding the roots of an oscillatory function that is only known through its differential equation is easily done, thanks to the "EventLocator" functionality built into NDSolve[]. Here's how to apply it to Andrew's example:

{osc, rts} = Reap[y /. First @ NDSolve[
                       {y''[x] - (x^2 - 3 x) y'[x] + 10^4 y[x] == 0, y[0] == 3, y'[0] == 1},
                       y, {x, 0, 5}, 
                       Method -> {"EventLocator", "Event" -> y[x], "EventAction" :> Sow[x],
                                  Method -> "StiffnessSwitching"}]];
rts = First[rts];

One might want to verify that all the roots within the integration interval were captured. Here's one way:

Plot[osc[x], {x, 0, 5}, Axes -> None, Frame -> True,
     Epilog -> {Directive[Red, AbsolutePointSize[2]], Point[{#, osc[#]} & /@ rts]}]

plot of oscillatory function and roots

Having done this, here's how to evaluate the integral $\int_0^5 x^2\mathtt{osc}[x]\;\mathrm dx$:

Total[NIntegrate[x^2 osc[x], {x, ##}] & @@@ Partition[Union[Flatten[{0, rts, 5}]], 2, 1],
      Method -> "CompensatedSummation"]

In the interest of making my post a lot less useless than it seems to be, here's how to adapt the approach given above when the function multiplying the osc[x] in fpghost's answer is also oscillatory (note the increased WorkingPrecision setting):

{osc, rts} = Reap[y /. First@NDSolve[{
                  y''[x] + (2 (x - 1))/(x (x - 2)) y'[x] + ((100 x^2)/(x - 2)^2 -
                  2/(x (x - 2))) y[x] == 0, y[5] == 3, y'[5] == 1},
                  y, {x, 5, 10},
                  Method -> {"EventLocator",
                             "Event" -> {Re[Exp[-I x] y[x]], Im[Exp[-I x] y[x]]},
                             "EventAction" :> Sow[x],
                             Method -> "StiffnessSwitching"}, WorkingPrecision -> 25]];

rts = Union[First[rts], SameTest -> (Chop[#1 - #2] == 0 &)];

Total[NIntegrate[Exp[-I x] osc[x], {x, ##}, WorkingPrecision -> 20] & @@@
      Partition[Union[Flatten[{5, rts, 10}]], 2, 1], Method -> "CompensatedSummation"]
   -0.048159751342842237133 + 0.045326948711103488692 I
  • $\begingroup$ Looks like a nice method, so by splitting between the roots then individually NIntegrating these regions and summing up at end we expect better accuracy for some reason? $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 8:54
  • $\begingroup$ The idea is that since you're integrating over a region where the integrand is always nonnegative or always nonpositive, the component integrals are computed pretty accurately, and any subtractive cancellation will only manifest when you sum up these integrals. $\endgroup$ Commented Oct 17, 2012 at 9:11
  • $\begingroup$ I see, I really like the idea, I will keep it in mind for future. Unfortunately for me I am not directly integrating the solution to my ODE, (the integrand is more like Exp[-is] (g(s) osc1[C]osc1*[C-s]+h(s) osc2[C]osc2*[C-s]) , where C is some constant, g(s),h(s) are just functions of integration var, and osc1, osc2 are two distinct solutions to the ODE for different ICs) so the roots of osc don't always correspond to roots of the ODE.sigh. $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 9:18
  • $\begingroup$ Then, you might be interested in this for finding the zeroes... (the Exp[-I s] is oscillatory in itself, so the zeroes of your integrand will certainly not correspond to the zeroes of the solution(s) of your DE). Could you maybe edit your question to mention your real problem? $\endgroup$ Commented Oct 17, 2012 at 9:28
  • $\begingroup$ I updated with the schematic. $\endgroup$
    – fpghost
    Commented Oct 17, 2012 at 18:11

Here is a very simple way of expressing my problem and gives the answer I found for this setup:

The ODE:

  eq1 := y''[x] + (2 (x - 1))/(x (x - 2))
  y'[x] + ((100 x^2)/(x - 2)^2 - 2/(x (x - 2))) y[x] == 0;

Solve it:

osc = y /. NDSolve[{eq1, y[5] == 3, y'[5] == 1}, y, {x, 5, 10}] // 

Plot[osc[x], {x, 5, 10}]


Give a Levin amplitude

f[x_] := Exp[-I x];

Compare against default methods

NIntegrate[f[x] osc[x], {x, 5, 10}] // Timing

{4.96831,-0.0481038+0.0453335 I}

NIntegrate[f[x] osc[x], {x, 5, 10}, 
Method -> {"LevinRule", "AdditiveTerm" -> 0, 
"Amplitude" -> {f[x], 0}, "Kernel" -> {osc[x], osc'[x]}, 
"DifferentialMatrix" -> {{0, 
   1}, {-((100 x^2)/(x - 2)^2 - 2/(x (x - 2))), (-2 (x - 1))/(
   x (x - 2))}}}] //Timing

{133.84, -0.0481038 + 0.0453335 I}

So I believe that is the Levin implementation of near enough my real problem, it appears slower for this set of parameters, but one can see that putting the '100' in the ODE to '1000' Levin comes into its own.

The problem I now have is that my 'amplitude' factor can also be highly oscillatory. Using the definitions in the 'Update' part of my OP I have

 f[x_,w_]:=Exp[-I (-5-x)] (ubar[x])^(-i 4 w)Exp[+i w rstar[x]]

as an amplitude where x is typically negative.


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.