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