Why do I get number with Precision larger than error estimate?

Question

I'm trying to integrate a large function like this

func[ξ_] = 1/(1167598657575 π Sqrt[1 - ξ]) 4 Sqrt[2/143] ξ Sqrt[(1 - ξ)/(1 + ξ)] (1 + 72 (-2 + 2 ξ^2) + 1197 (-2 + 2 ξ^2)^2 + 7980 (-2 + 2 ξ^2)^3 + 209475/8 (-2 + 2 ξ^2)^4 + 92169/2 (-2 + 2 ξ^2)^5 + 706629/16 (-2 + 2 ξ^2)^6 + 43263/2 (-2 + 2 ξ^2)^7 + 1081575/256 (-2 + 2 ξ^2)^8) ((8957659658839916544 - 38598804293820612608 ξ^2 + 69110375109906923520 ξ^4 - 66520053962590126080 ξ^6 + 37182509594538147840 ξ^8 - 12158541247606161408 ξ^10 + 2222643299988013056 ξ^12 - 202385520787169280 ξ^14 + 7047164723274720 ξ^16 - 40061762318655 ξ^18) ((1 + ξ) EllipticE[(2 ξ)/(1 + ξ)] - EllipticK[(2 ξ)/(1 + ξ)]) + 8 ξ^2 (279926864338747392 - 1127483203586621440 ξ^2 + 1863097027060039680 ξ^4 - 1626763404855214080 ξ^6 + 805078969740165120 ξ^8 - 224757138041339904 ξ^10 + 33055899621212160 ξ^12 - 2165655547777920 ξ^14 + 40988123520795 ξ^16) EllipticK[(2 ξ)/(1 + ξ)]);

in the following modes (error estimate extraction is from this answer):

arb16 = NIntegrate[func[x], {x, 0, 1}, WorkingPrecision -> 16, 
  IntegrationMonitor :> ((errors = Through[#@"Error"]) &)]
arb200 = NIntegrate[func[x], {x, 0, 1}, WorkingPrecision -> 200]

I get the numbers, and now try to get the precision of arb16. Here's what I have:

Row[{"arb16 precision: ", -Log10@Abs[arb200 - arb16],
  ", precision by error estimator: ", -Log10@Total@errors,
  "\nreturned precision ", Precision@arb16,
  ", accuracy ", Accuracy@arb16}]

arb16 precision: 10.3694098, precision by error estimator: 8.9544353783245526

returned precision 16., accuracy 17.7955803201903

So, the error estimate by NIntegrate appears correct, but Precision and Accuracy of the number are still way too optimistic. Why is it so? Is there a way to force the correct precision numbers to be given out by NIntegrate?

Answer 1

Consider any numerical integration method $I^*(f,a,b)$ that approximates the exact integral $I$ of a function $f$ over an interval $[a,b]$. It will be implemented by a computation represented by, say, I[f[x], {x, a, b}]. Conceiving an NIntegrate[] command in this way breaks the error into two independent components. The first is the error of the method $I^* - I$, where $I^* = I^*(f,a,b)$ is assumed to be computed exactly; this error is often called the "truncation error." The second is the difference between $I^*$ and I[f[x], {x, a, b}]; this error is usually thought of as "rounding error."

The error estimator of an NIntegrate[] method endeavors to yield an upper bound on the truncation error. The options PrecisionGoal and AccuracyGoal set the goal for this bound. The error estimate is often done by comparing the method with a less accurate method, e.g. in the "GaussKronrodRule" and the "ClenshawCurtisRule". For a well-behaved function, using either of these rules, the error can be greatly overestimated, as noted in this answer.

The option WorkingPrecision, which controls the Precision[] and Accuracy[] of the result, is used to control rounding error. (It also affects the default PrecisionGoal, which affects the bound on the truncation error discussed above and which @ilian also notes in his answer.) Precision@arb16 and Accuracy@arb16 reflect this rounding error. Normally this, too, is an upper bound.** So Precision[] has little to do with the truncation error in approximating the integral.

The general approach in numerical analysis has tended to focus on bounding error. Mathematica's error estimators and arbitrary-precision numerics follow this approach. I don't believe there is a way to get a better estimate of the error, without knowing the exact value of the integral or computing it to a higher precision. (One could perhaps use a rule with a good error estimator; the approach in the rule chebRule in the linked answer above does somewhat better on a function that is analytic on {0, 1} than "GaussKronrodRule", but the OP's function has a singularity at 1.)

---

**Technical note: The precision-tracking rules used to keep track of the precision of a computation are based on linear approximation of the error. They do not produce an upper bound in themselves but aim to be accurate. However, Mathematica keeps extra guard bits in the internal representation of arbitrary-precision numbers, so that the actual internal number computed usually has a rounding error of much less magnitude than the Precision[] and Accuracy[].

Answer 2

[Too long to type as a comment, but here is why PrecisionGoal is relevant]

Any number has a precision, and Precision[arb16] reflects only the uncertainty in the value of arb16 as a number. It is determined by both the precision of the input, and all the arithmetic operations subsequently done with it (see rounding error in Michael E2's excellent answer).

What it isn't necessarily related to is the error of the integration method, that is Precision[arb16] has nothing to do with how correct the obtained integral value is.

Since you only specified WorkingPrecision to be 16, but left the precision goal at the default setting (Automatic, which in this case is effectively half of the working precision), you should only expect the result to approximate the value of the integral up to (at least) 8 digits of precision. Since arb16 certainly satisfies this, NIntegrate has fulfilled its contract.

The way to tell NIntegrate you want a result with more correct digits is to specify PrecisionGoal, which, as documented, sets convergence criteria by putting an upper bound on the estimated integration error (aka truncation error in Michael E2's answer). For an example,

arb16a = NIntegrate[func[x], {x, 0, 1}, WorkingPrecision -> 16, PrecisionGoal -> 16]

(* -0.01601104510443394 *)

does indeed have (no less than) 16 correct digits.