Problem with numerical evaluation of analytically solved integral, solution way off

Question

The following command in Version 9.0.1:

N[Integrate[x^50*Sin[x], {x, 0, 1}]]

gives $1.4615\times 10^{48}$ which is way off from the correct solution which is between $0$ and $0.5$.

Interestingly the indefinite integral comes out correct and if you increase the precision Mathematica gives an error and the following number: $0.0163$. The same result is given by NIntegrate.

What exactly is the problem here, is this a known issue and is there a workaround?

Answer 1

This is standard cancellation error (please Google for it). See below for more details.

---

Take a look at the exact result first:

result = Integrate[x^50*Sin[x], {x, 0, 1}]

(* ==>
    16432804687774250383441481995831940788236063969597816674967907249 Cos[1] + 
      50 (-608281864034267560872252163321295376887552831379210240000000000 + 
           511851539169698127179811081351899421435935914235737810293853873 Sin[1])

 *)

Expand[result]

(* ==>
    -30414093201713378043612608166064768844377641568960512000000000000 + 
     16432804687774250383441481995831940788236063969597816674967907249 Cos[1] +
     25592576958484906358990554067594971071796795711786890514692693650 Sin[1]
 *)

This is approximately 0.016, a very small number, which arises as a difference of two very large numbers. So now it's clear that we need to work with a very high precision (lots of digits) when carrying out the subtraction, otherwise we'll get what's called catastrophic cancellation and thus an incorrect result. This is a phenomenon you need to learn about and be aware of if you do any kind of numerical work on a computer. Mathematica can do computations with lots of digits to avoid this problem, but systems which always use the CPU's built in floating point arithmetic (machine precision) won't be able to work around the problem simply by increasing precision.

By default, Mathematica uses machine precision when you apply N[...], i.e. approximately 16 digits (strictly, it uses 64-bit IEEE floating point, which implies a 53-bit significand, so the number of decimal digits is $\log_{10} 2^{53}$). In this example you have numbers that have 64 digits to the left of the decimal point, so it's impossible to get a precise result if the result would have the first significant digit two places to the right of the decimal point. That is, even a value as large as $10^{48}$ is as near to zero as makes no difference! Crucially, however, precision tracking is disabled when working with machine numbers (for performance reasons), so that even if catastrophic precision loss occurs, Mathematica has no way of knowing it.

The solution is to use arbitrary precision arithmetic, which you can do by explicitly requesting a precision in N[...]. Mathematica will then automatically increase the number of digits used internally in order to get a result that is precise to at least the number of digits requested.

Let's try it:

N[Expand[result], 10]

(* ==> 0.01628978362 *)

This works fine. But if we use the unexpanded form, it doesn't:

N[result, 10]

(* ==> N::meprec: Internal precision limit $MaxExtraPrecision = 50.`
       reached while evaluating ... *)

Mathematica warns you that after internally increasing the precision to the requested precision plus 50 digits (i.e. a total of 60 digits in this case), it still can't get a precise result (so you know that the result it prints might not be correct). We can allow more than 50 digits of precision increase:

Block[{$MaxExtraPrecision = 100}, N[result, 10]]

(* ==> 0.01628978362 *)

Now it works fine.

Note that this precision adjustment is part of the functionality of N, and it will not be done automatically e.g. when carrying out simple arithmetic operations. But, since precision is always tracked for arbitrary-precision values, one can easily check to see if precision loss has occurred by using Precision.

If you use NIntegrate, the system never actually computes an exact closed form solution. It uses numerical methods which in this case will give a precise result even without needing to use a high working precision:

NIntegrate[x^50*Sin[x], {x, 0, 1}]

(* ==>  0.0162898 *)

Answer 2

NIntegrate[x^50*Sin[x], {x, 0, 1}] works fine and results to 0.0162898.

Also putting N[Integrate[x^50*Sin[x], {x, 0, 1}], 20] you get a correct result.

If you try to calculate the anti-derivative using :

anti[x_] := Module[{t}, Integrate[t^50*Sin[t], t] /. t -> x] then you will see that you get the exact integral form.

Now

anti[1]==16432804687774250383441481995831940788236063969597816674967907249 Cos[1] + 25592576958484906358990554067594971071796795711786890514692693650 Sin[1]

and

anti[0]==30414093201713378043612608166064768844377641568960512000000000000

N[anti[1] - anti[0], 20]==0.016289783620195801683 but N[anti[1]-anti[0]]==0. which is caused by the round to only 8 digits of Sin[1] and Cos[1].

Try Trace[N[Integrate[x^50*Sin[x], {x, 0, 1}]]] to see where exactly it goes way off!