highly oscillatory integral

Here's a manual implementation of the extrapolating oscillatory strategy. NIntegrate balks at trying it because the cosine is buried inside other functions.

ClearAll[psum];
psum[i_?NumberQ, opts___] := 
 NIntegrate[(Sign[Cos[x]] Sqrt[Abs[Cos[x]]])/(x + 100),
  {x, (i + 1/2) π, (i + 3/2) π}, opts];

res = NSum[psum[i], {i, 0, ∞},
    Method -> {"AlternatingSigns", Method -> "WynnEpsilon"},
    "VerifyConvergence" -> False] +
  NIntegrate[(Sign[Cos[x]] Sqrt[Abs[Cos[x]]])/(x + 100), {x, 0, (1/2) π}]
(*
  0.00010967(-)
*)

This agrees with Chip Hurst's result to almost 5 digits. But the method is reasonably quick, so we can examine the convergence as the WorkingPrecision (and therefore the PrecisionGoal) is raised:

ClearAll[res];
mem : res[wp_: MachinePrecision] := mem =
  NSum[psum[i, WorkingPrecision -> wp], {i, 0, ∞}, 
     Method -> {"AlternatingSigns", Method -> "WynnEpsilon"}, 
     "VerifyConvergence" -> False, WorkingPrecision -> wp] + 
   NIntegrate[(Sign[Cos[x]] Sqrt[Abs[Cos[x]]])/(x + 100),
     {x, 0, (1/2) π}, WorkingPrecision -> wp];

Table[res[wp], {wp, 16, 40, 8}]
(*
  {0.0001096696058468,
   0.000109676059727856341260, 
   0.00010967605972782927874728238851, 
   0.0001096760597278292787470847687426733221}
*)

% // Differences
(*
  {6.4538811*10^-9, -2.7062513*10^-17, -1.9761977*10^-25}
*)

At each step we extend the number of digits that agree, and we can see that the result seems reliable.


Alternatively, one can multiply and divide by cosine, which enables NIntegrate to use the "ExtrapolatingOscillatory" strategy. It's a bit slower because it introduces (integrable) singularities in the non-oscillatory part of the integrand, but again we get similar results.

ClearAll[res2];
mem : res2[wp_: MachinePrecision] := mem =
  NIntegrate[Cos[q] (1/Sqrt[Abs[Cos[q]]])/(q + 100), {q, 0, Infinity},
    Method -> "ExtrapolatingOscillatory", 
    MaxRecursion -> Ceiling[10 + 1.5 wp],  (* increase in MaxRecursion for singularities *)
    WorkingPrecision -> wp]

res2[]
(*
  0.000109673
*)

Table[res2[wp], {wp, 16, 40, 8}]
(*
  {0.0001096731867831,
   0.000109676059727866251185, 
   0.00010967605972782927874731895632, 
   0.0001096760597278292787470847687424806474}
*)

I think the answer is about 0.000109672. This throws no messages, so I think we can assume it's accurate (?).

cres = NIntegrate[Sqrt[Cos[q]]/(q + 100), {q, 0, ∞}, Method -> "LocalAdaptive"]
0.171013 + 0.170903 I
Re[cres] - Im[cres]
0.000109672

Here's some verification by integrating piecewise.

partials = Prepend[ParallelTable[
  Total[{1, -1} ReIm[NIntegrate[Sqrt[Cos[q]]/(q+100), {q, (2k+1)π/2, (2k+3)π/2}]]], 
  {k, 0, 100000}], 
  NIntegrate[Sqrt[Cos[q]]/(q+100), {q, 0, π/2}]];

Total[partials]
0.000105864
ListPlot[Accumulate[partials], 
  GridLines -> {None, {Total[partials]}}, 
  GridLinesStyle -> Directive[Red, Opacity[0.5], Dashed]]

enter image description here

So I think the first result looks to be fairly accurate, but I'm not sure by how much. My guess would be by at least 4-5 digits. I'd be interested to see what an expert in this area has to say.


another approach: rearrange the integral as:

NIntegrate[ Sqrt[Cos[q]]/(q + 100), {q, 0, Pi/2}] +
 Sum[ NIntegrate[ Sqrt[Sin[qp]]
    (1/(qp + 100 + Pi + (2 n + 1/2) Pi) - 
      1/ (qp + 100 + (2 n + 1/2) Pi)),
   {qp, 0 , Pi }], {n, 0, Infinity}]

(Does not eval, but you can put in a big number for Infinity to get a good approximation), eg at 10000 we get 0.000128712

Now take the infinite sum inside the integral:

p[qp_] = Assuming[{0 < qp < Pi}, 
  Sum[(1/(qp + 100 + Pi + (2 n + 1/2) Pi) - 
     1/ (qp + 100 + (2 n + 1/2) Pi)), {n, 0, Infinity}]]

(PolyGamma[0, (200 + [Pi] + 2 qp)/(4 [Pi])] - PolyGamma[0, (200 + 3 [Pi] + 2 qp)/(4 [Pi])])/(2 [Pi])

NIntegrate[ Sqrt[Cos[q]]/(q + 100), {q, 0, Pi/2}] +
 NIntegrate[ Sqrt[Sin[qp]] p[qp], {qp, 0 , Pi }]

0.000109676

with WorkingPrecision->50

0.000109676059727829278747084768743039807725993572325

further if we let that 100 be a parameter m and use a series about m->Infinity :

Assuming[{m > 0},
  Integrate[Sqrt[Cos[q]] (1/m - q/m^2), {q, 0, Pi/2}] +
   Integrate[
    Sqrt[Sin[qp]] (-(1/(2 m)) + qp/(2 m^2) ), {qp, 0, 
     Pi}]] // Simplify

we get an analytic approximation:

(3 Sqrt[2 [Pi]] Gamma[3/4]^2 - 4 HypergeometricPFQ[{1/2, 1/2, 1}, {3/2, 7/4}, 1])/(6 m^2)

 %/. m->100

0.000109743