# highly oscillatory integral

Hi am trying to compute the following integral in Mathematica:

    NIntegrate[(Sign[Cos[q]] Sqrt[Abs[Cos[q]]])/(q+100),{q,0,Infinity}]


But when I run that command I get the following error:

    NIntegrate::ncvb: NIntegrate failed to converge to prescribed
accuracy after 9 recursive bisections in q near {q} = {6.3887*10^56}.
NIntegrate obtained -317.031 and 171.60663912222586 for the integral
and error estimates.


Is there a way to solve this issue? Thanks

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}
*)

• For people who are interested in manually implementing the Longman approach like what Michael did here, note that you always have the option of either integrating between the roots of the oscillatory function, or integrating between extrema. Both should result in an easily extrapolated alternating series. – J. M. will be back soon Jul 26 '17 at 5:08

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


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