# NIntegrate not giving the correct value

I have an integrand which cannot be integrated analytically:

pr[x_] = -Pi/4 x^2 Cos[x^2] FresnelC[Sqrt[2] x/Sqrt[Pi]]^2;


but I also have an approximation to it:

g[x_] = Sqrt[\[Pi]/2]/(16 x) - Cos[x^2]/32 + (3 Cos[x^2])/(128 x^4) -
1/16 \[Pi] x^2 Cos[x^2] + (Sqrt[\[Pi]/2] Cos[2 x^2])/(16 x) +
1/32 Cos[3 x^2] - (7 Cos[3 x^2])/(128 x^4) - (3 Sin[x^2])/(
128 x^6) + Sin[x^2]/(32 x^2) + (3 Sqrt[\[Pi]/2] Sin[2 x^2])/(
32 x^3) - 1/8 Sqrt[\[Pi]/2] x Sin[2 x^2] - (3 Sin[3 x^2])/(
128 x^6) + Sin[3 x^2]/(32 x^2);


I try to estimate the error between the integrand and its approximation by doing:

NIntegrate[pr[x] - g[x], {x, 10, ∞},
Method -> "MonteCarlo",
PrecisionGoal -> 10]


but it gives me a large value (-7e8) whereas I know that it should be very small.

How do I correct the numerical integration?

• I dont get what you want to do. Is $g(x)$ and approximation to $pr(x)$? If yes, based on what approach? Why do you want to evaluate the error with an integral? Is this a question about the software Mathematica or is it a math question? Please clarify. – Mauricio Fernández Jun 24 '17 at 12:43
• g(x) is an approximation to pr(x), yes, based on expanding the Fresnel function for large x. Visually, the two functions overlap for x > 1, so the difference between them must converge. I want to verify this by numerically integrating their difference between 10 and infinity. – Morgan Jun 24 '17 at 13:10
• This is not a Mathematica problem since the integral of pr[x_] over the interval {x, 10, [Infinity]} diverges. – user64494 Jun 24 '17 at 17:13
• @user64494 Because g(x) is a Taylor expansion of pr(x) at infinity. Also, someone else has done the same calculation before, only with less terms in the expansion, and managed to integrate it with Mathematica to get a converging result. Unless they published an error in a scientific journal, it should still be valid. – Morgan Jun 24 '17 at 17:40
• Sure, g[x_] is obtained with the following code: -Pi/4 x^2 Cos[ x^2] Series[ FresnelC[Sqrt[2] x/Sqrt[Pi]], {x, [Infinity], 5}]^2 // Normal // TrigReduce // Expand – Morgan Jun 24 '17 at 19:18

The integrand, as an exact function, seems a reasonably well-behaved oscillatory function that can be integrated (numerically) successfully via extrapolation; however, as a numerical function, it is horrible. (Expect the code below to run slowly, from seconds to minutes.) To evaluate at x -> 100. accurately you need 4343 digits of accuracy, or you'll be way off:

pr[x] - g[x] /. {{x -> 100.4342}, {x -> 100.4343}};
% // N
Precision /@ %%
(*
{0.00136186, -5.3138*10^-11}
{26.7328, 4323.03}
*)


All I think of is that there is some internal catastrophic cancellation that goes undetected. Note also that because of the need

Another possible issue, dwarfed by the precision one, is that the linear system constructed by the Levin rule is of dimension 48. Mathematica balks at solving this system at high precision. Perhaps with a significantly higher "TimeConstraint", it might work.

The built-in "ExtrapolatingOscillatory" strategy seems a good candidate, but it turns out that it cannot parse the complicated integrand. However we can imitate, as was done in How to find a more precise value of integral?

Borrowing again from Anton Antonov's answer, we can divide the interval of integration up at x == Sqrt[2 Pi * i] and represent it as a sum integrals between the nodes. I put a Check[] on the integral so I could investigate ones that did not converge. We'll also be using the terms repeatedly, so it's good to memoize their values.

ClearAll[term];
mem : term[i_?NumberQ] := mem = Module[{iRes},
Check[
iRes = NIntegrate[(pr[x] - g[x]), {x, Sqrt[2 Pi*i], Sqrt[2 Pi*(i + 1)]},
PrecisionGoal -> 20, AccuracyGoal -> 20,
WorkingPrecision -> 40 + 3 i,
Method -> {"GaussKronrodRule", "Points" -> 15}],
Print["i = ", i]
];
iRes];


To make it easy to code, I picked x == Sqrt[2 Pi*i], i = 16, 17, ... for the nodes, but it leaves a little bit of the OP's integral from x == 10 to x == Sqrt[2 Pi*16] to be added:

bit = NIntegrate[pr[x] - g[x], {x, 10, Sqrt[2 Pi*16]},
PrecisionGoal -> 10, WorkingPrecision -> 200]
sf[n_] := Total@Table[SetPrecision[term[i], 200], {i, 16, n}]


We can then use Anton's SearchSumValue[] function (see below) to approximate the integral.

res = SearchSumValue[
sf,                        (* partial sum function *)
12,                        (* desired accuracy: abs. error < ~10^-12 *)
30,                        (* number of terms to add each step *)
40,                        (* starting number of terms *)
Richardson[#1, #2, 24] &]  (* extrapolation function *)
answer = res[[2]] + bit      (* add the little bit to get the complete integral *)
(*
-7.4084729050246807289001990592571553616689002705033883146404497941764\
1713462740784131023302538755726602927879849993920908117244826729334375\
32927540835628437772008149977875*10^-6
*)


So the integral is estimated to be

-7.40847290*^-6


with an absolute error estimated by SearchSumValue[] to be at most 2.*^-15.

Further remarks (updated): The amplitude of pr[x] - g[x] appears to be about (15 Sqrt[(π/2)] )/(32 x^5), which can be seen from Simplify@Normal@Series[pr[x] - g[x], {x, Infinity, 5}]. I did a similar analysis as above with pr[x] - g[x] + (15 Sqrt[(π/2)] )/(64 x^5) with a different node grid and got the same result for the integral (subtracting the integral of (15 Sqrt[(π/2)] )/(64 x^5) from the sum, of course). Here is a comparison with the "amplitude" (15 Sqrt[(π/2)] )/(32 x^5) (takes a while to run):

Plot[{RealExponent[g[x] - pr[x]], RealExponent[(15 Sqrt[(π/2)] )/(32 x^5)]},
{x, 10, 100}, MaxRecursion -> 2, PlotRange -> All, WorkingPrecision -> 5000,
PlotLegends -> {HoldForm@Log10@Abs[pr[x] - g[x]], HoldForm@Log10[(15 Sqrt[(π/2)] )/(32 x^5)]}]


Code dump:

Clear[Richardson]
Richardson[A_, n_, N_] :=
Total@Table[(A[n + k]*(n + k)^N*
If[OddQ[k + N], -1, 1])/(k! (N - k)!), {k, 0, N}];

Clear[SearchSumValue]
SearchSumValue[partialSumFunc_, accGoal_, step_: 10, startStep_: 40, methodFunc_: Shanks] :=
Block[{res, pf = 2},
res = NestWhile[{#[[1]] + step, #[[3]],
methodFunc[partialSumFunc, #[[1]] + step]} &, {1,
methodFunc[partialSumFunc, startStep],
methodFunc[partialSumFunc, startStep + step]},
Abs[N[#[[2]], pf*accGoal] - N[#[[3]], pf*accGoal]] > 10^-accGoal &];
(*n-steps,estimate,error*)
Append[res[[{1, 3}]],
Abs[N[res[[2]], pf*accGoal] - N[res[[3]], pf*accGoal]]]];


I'm adding this as a CW answer because (1) it shows how one can probably solve the problem within NIntegrate (if you allow me to split the integrand into two parts), (2) it's a fairly bad way to go about it, and (3) it was interesting to see that it could probably work, even if it is interesting only to me. The comments in my other answer about how bearish the integrand is with respect to numerical precision still apply and are probably the root of the difficulty of the integral.

The trick

tmp = pr[x] - g[x] // TrigExpand;
tmp = tmp /. Sin[x^2]^2 -> 1 - Cos[x^2]^2 // Expand;
iC = Cos[x^2] (Plus @@
Cases[tmp, f_ /; ! FreeQ[f, Cos[x^2]] :> f/Cos[x^2]]);
tmp = tmp - iC // Expand;
iS = Sin[x^2] (Plus @@
Cases[tmp, f_ /; ! FreeQ[f, Sin[x^2]] :> f/Sin[x^2]]);

OwnValues[iC]
OwnValues[iS]
HoldForm[pr[x] - g[x] - (iC + iS)] -> pr[x] - g[x] - (iC + iS) // Simplify
(*
{HoldPattern[iC] :>
Cos[x^2] (1/8 - 3/(16 x^4) + (π x^2)/16 - (Sqrt[π/2] Cos[x^2])/(8 x) - 1/8 Cos[x^2]^2 + (7 Cos[x^2]^2)/(32 x^4) - 1/4 π x^2 FresnelC[Sqrt[2/π] x]^2 - (3 Sqrt[π/2] Sin[x^2])/(16 x^3) + 1/4 Sqrt[π/2] x Sin[x^2] + (9 Cos[x^2] Sin[x^2])/(128 x^6) - (3 Cos[x^2] Sin[x^2])/(32 x^2))}

{HoldPattern[iS] :>
Sin[x^2] (3/(128 x^6) - 1/(32 x^2) - (3 Sin[x^2]^2)/(128 x^6) + Sin[x^2]^2/(32 x^2))}

pr[x] - g[x] - (iC + iS) -> 0
*)


Then the integral can be computed as the sum of two calls to NIntegrate for a chosen "precision goal" pg and working precision wp:

iCRes = NIntegrate[iC, {x, 10, Infinity}, PrecisionGoal -> pg,
WorkingPrecision -> wp, Method -> "ExtrapolatingOscillatory"];
iSRes = NIntegrate[iS, {x, 10, Infinity},
Method -> "ExtrapolatingOscillatory"]
iCRes + iSRes


Now iCRes, which contains FresnelC[], is very difficult to compute, and it seems impossible to achieve the precision goal pg. But iSRes is quite easy, and an accurate result is obtained with machine precision. Below is a table of results for increasing pg and wp. It seems that the result is, or might be, converging to the result -7.4 * 10^-6 in my other answer. It's hard to have any confidence in that without knowing the other answer though, since a precision goal in the low 100s produces a result with only one digit of accuracy.

Perhaps with a precision goal of 500 and 2-3 days of computation, two digits of accuracy might be achieved.