Is this caused by round-off errors?

Question

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?

Answer 1

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.

Answer 2

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