Skip to main content
4 of 8
added 931 characters in body
fpghost
  • 2.1k
  • 2
  • 21
  • 25

Strategies to solve an oscillatory integrand only known numerically

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];
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.

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.

fpghost
  • 2.1k
  • 2
  • 21
  • 25