Strategies to solve an oscillatory integrand only known numerically

Question

I have an integrand that looks like this:

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];
M=1;
rinf=10000;
rH = 200001/100000;
r0 = 10;
wp=40;
ac=wp-8;
\[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.

UPDATE

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

vbar[tau_?
NumericQ] := (4 M) ((tau/tauh)^(1/3) + 1) Exp[-(tau/tauh)^(1/3) + 
 1/2 (tau/tauh)^(2/3) - 1/3 (tau/tauh)];
ubar[tau_?
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.

Answer 1

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}]

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.

Answer 2

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]}]

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"]
   -2.802321164674166

---

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

Answer 3

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}] // 
First;

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.