# Strategies to solve an oscillatory integrand only known numerically

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.46417535429324425086933619669164038626679130.*^-6 -
8.95048324839030667077034540604783599393166530.*^-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. - What about rewriting your differential equations to directly give the integral from NDSolve? – celtschk Oct 16 '12 at 17:34 Have you tried a Fourier on that ? – belisarius Oct 16 '12 at 17:34 Have you seen this video ? It talks about similar problem and hybrid numeric-symbolic methods to address it. – Vitaliy Kaurov Oct 16 '12 at 17:43 "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. – Andrew Moylan Oct 17 '12 at 2:23 @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). – celtschk Oct 19 '12 at 7:55 ## 3 Answers 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. - 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" ? – fpghost Oct 17 '12 at 8:08 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. – fpghost Oct 17 '12 at 8:12 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? – fpghost Oct 17 '12 at 18:13 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 – fpghost Oct 17 '12 at 18:21 If the oscillatory part of your integrand satisfies a nonlinear rather than linear ODE, then LevinRule can't be applied. – Andrew Moylan Oct 18 '12 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]}]  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

-
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? –  fpghost Oct 17 '12 at 8:54
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. –  Ｊ. Ｍ. Oct 17 '12 at 9:11
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. –  fpghost Oct 17 '12 at 9:18
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? –  Ｊ. Ｍ. Oct 17 '12 at 9:28
I updated with the schematic. –  fpghost Oct 17 '12 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}] //
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 params 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 defintions 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.

-