2
$\begingroup$

Let's consider this integration

Integrate[ E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}]

It returns Gamma[4 n μ] Gamma[4 n ν] Hypergeometric1F1Regularized[ 4 n ν, 4 n (μ + ν), 4 n s]

Now let's replace those letter by numbers

n = 20000; s = 0.01; μ = 10^-3; ν = 10^-3;

And let's run both the result of the integral and the integral itself

In[28]:= Integrate[ E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}] Out[28]= -1.786676093655969*10^352 In[29]:= Gamma[4 n μ] Gamma[ 4 n ν] Hypergeometric1F1Regularized[4 n ν, 4 n (μ + ν), 4 n s] Out[29]= 5.20048*10^228

The results differ. While I agree both look like big numbers, one is 124 orders of magnitude higher than the other one.

Given that all calculations were analytic (no numerical approximations) we should find the same results. What causes these different results? Is this difference caused by round-off errors? How can I solve such issues?

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ Yes, I believe it is round-off error. Try s = 0.01`20 and then s = 0.01`200. $\endgroup$ Commented Jan 21, 2015 at 0:58
  • $\begingroup$ Btw, what version are you using? I get different results (V10.0.2, Mac OSX 10.10.1). $\endgroup$ Commented Jan 21, 2015 at 1:00
  • $\begingroup$ with s = 0.0120` I get 0.*10^350 (in red) for the first formula and 5.2004811325650457130*10^228 for the second. I am using Version 10.0.1.0 in Mac OSX. I think I had a mistake first. Sorry about that. I corrected my post. Do we still get different results? $\endgroup$ Commented Jan 21, 2015 at 1:14
  • $\begingroup$ with s = 0.01200` I get in both cases some long number times 10 to the power 228. Seems much better. Was it indeed a round-off error? And the solution is to make to always write the good number of decimal in inputs so to tell Mathematica what degree of precision it should use? $\endgroup$ Commented Jan 21, 2015 at 1:17
  • $\begingroup$ Now I get the same numbers as in the update. The red in 0.*10^350 means that the precision loss was so great that M can only estimate that the result is closer to zero than 10^350 -- a significant loss of precision. Precision tracking does not occur with MachinePrecision numbers such as s = 0.01. Using arbitrary-precision numbers such as 0.01`20 etc, which turns on precision tracking, is a way to check calculations that seem off. I usually increase precision until the result stabilizes. $\endgroup$ Commented Jan 21, 2015 at 1:47

2 Answers 2

Reset to default
2
$\begingroup$

The problem is with numerical error in the procedure used by Integrate. Sufficiently high arbitrary precision, or exact input, is needed to ensure an accurate result from Integrate. On the other hand, MachinePrecision is sufficient for NIntegrate. There is nothing particularly numerically challenging about the integrand, so the NIntegrate result is not surprising.

Finally, applying N to the exact result to any precision yields an accurate result to the requested precision, as promised in the documentation. A high precision is not required.

Insufficient input precision

Block[{n = 20000, s = 0.01, μ = 10^-3, ν = 10^-3}, Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}] ] (* -1.786676093655969*10^352 *) Block[{n = 20000, s = 0.01`20, μ = 10^-3, ν = 10^-3}, Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}] ] (* 0.*10^350 *)

Given that the answer is about 10^228, this result suggests that the input needs more than 350 - 228 == 122 digits of precision to get at least one digit of accuracy in the result.

Sufficient input precision

122 more digits would be 142; 145 yields 3 digits.

Block[{n = 20000, s = 0.01`145, μ = 10^-3, ν = 10^-3}, Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}] ] (* 5.20*10^228 *)

Exact input yields an exact answer; applying N gives an accurate approximation.

Block[{n = 20000, s = 1/100, μ = 10^-3, ν = 10^-3}, N[Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}], 6] ] (* 5.20048*10^228 *)

NIntegrate gives an accurate answer quickly. It should be the preferred method for an approximate machine real result. 

Block[{n = 20000, s = 0.01, μ = 10^-3, ν = 10^-3}, NIntegrate[E^(4 n x s) (1 - x)^(-1 + 4 n μ) x^(-1 + 4 n ν), {x, 0, 1}] ] (* 5.20048*10^228 *)

Summary

NIntegrate should also be the first method used to check a result from Integrate for accuracy. Stability of the result under increasing precision can be used as a further check. These are fairly common checks one can find in answers and comments to integration problems on the site.

Improve this answer
$\endgroup$
1
$\begingroup$

As stated n the comments, the problem lies with the defnition s = 0.01 which causes all calculations to be done with machine precision. Higher precision is required, say 20.

Analytic solution

int1 = Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n mu) x^(-1 + 4 n nu), {x, 0, 1}]

ConditionalExpression[ Gamma[4 mu n] Gamma[4 n nu] Hypergeometric1F1Regularized[4 n nu, 4 n (mu + nu), 4 n s], Re[n nu] > 0 && Re[mu n] > 0]

Analytic solution evaluated with exact numbers then convered to an approximate number with precision of 20.

N[int1 /. {n -> 20000, s -> 1/100, mu -> 10^-3, nu -> 10^-3}, 20]

5.2004811325650457130*10^228

Integrated with exact numbers then convered to an approximate number with precision of 20

Integrate[E^(4 n x s) (1 - x)^(-1 + 4 n mu) x^(-1 + 4 n nu) /. {n -> 20000, s -> 1/100, mu -> 10^-3, nu -> 10^-3}, {x, 0, 1}] // N[#, 20] &

5.2004811325650457130*10^228

Using exact numbers in integrand then doing numerical ntegation wth working precision of 20

NIntegrate[ E^(4 n x s) (1 - x)^(-1 + 4 n mu) x^(-1 + 4 n nu) /. {n -> 20000, s -> 1/100, mu -> 10^-3, nu -> 10^-3}, {x, 0, 1}, WorkingPrecision -> 20]

5.2004811325650457130*10^228

Analytic solution evaluated with high precision numbers (minimum precision of 20).

int1 /. {n -> 20000, s -> 0.01`20, mu -> 10^-3, nu -> 10^-3}

5.2004811325650457130*10^228

Numeric integration of high precision integrand (minimum precision of 20).

NIntegrate[ E^(4 n x s) (1 - x)^(-1 + 4 n mu) x^(-1 + 4 n nu) /. {n -> 20000, s -> 0.01`20, mu -> 10^-3, nu -> 10^-3}, {x, 0, 1}]

5.20048*10^228

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.