Today in Calculus class I was bored so I decided to try and approximate $\pi$ by evaluating $ \left( \displaystyle \int_{-a}^a e^{-x^2} dx \right)^2$ on my calculator for larger and larger values of $a$.
However, in doing so, I noticed something peculiar. When I plugged in $a = 100$ it gave me a value of approximately $\pi$, as expected. However, when I plugged in $a = 1000$ I got an answer of about $2.7 \times 10^{-7}$. In fact, I was able to narrow it down to figure out that $a = 892.26$ and below (with reasonable assumption) gave me a correct value while $a = 892.27$ and above (with reasonable assumption) gave me some very small value like I found with $a = 1000$
What's going on here?
edit: I realize that this is a problem with the calculator's integration method. I am aware that larger values of $a$ should converge closer to $\pi$.
$\endgroup$ 33 Answers
$\begingroup$TI-84 integration algorithm is an adaptive Gauss-Kronrod 15-point rule, with apparently an absolute error threshold of $10^{-5}$ (TI-84 Plus and TI-84 Plus Silver Edition Guidebook, 2010, p. 41). When the absolute error estimate from the rule is greater than the threshold, the interval is subdivided and the integral recalculated over the intervals. The intervals are recursively subdivided until the error estimate is within the threshold.
For $a \le 892.26$, the initial error estimate is greater than the threshold, and recursive subdivision is applied. For $a \ge 892.27$, the initial error estimate is less than the threshold, and no subdivision is applied.
Below is the equivalent computation done with Mathematica. For the OP's integral, to reproduce the integral, I had, surprisingly, to apply a common strategy used to cancel out odd functions when integrating over symmetric intervals:$$\int_{-a}^{a} f(x) \; dx = \int_0^a {\left(f(x) + f(-x)\right)} \; dx$$Here is the code:
{nodes, (* vector of nodes for the interval {0, 1} *) wts, (* weights vector *) errwts} = (* error weights vector *) NIntegrate`GaussKronrodRuleData[7, MachinePrecision];
Length@nodes
(* 15 *)The following shows that $10^{-5}$ is the threshold:
fvec = Exp[-Rescale[nodes, {0, 1}, {0, 892.26}]^2]; (* vector of function values *)
fvec.wts (2 * 892.26) (* Dot product of function values and nodes scaled by interval length yields the integral *)
fvec.errwts (2 * 892.26) (* Scaled dot product with the error weights estimates the error *)
(* 0.0000100015 integral estimate 0.0000100015 error estimate
*)
fvec = Exp[-Rescale[nodes, {0, 1}, {0, 892.27}]^2];
fvec.wts (2 * 892.27)
fvec.errwts (2 * 892.27)
(* 9.99831*10^-6 integral estimate 9.99831*10^-6 error estimate
*)The following calculates the integral and error estimate. The integral agrees with the OP's:
fvec = Exp[-Rescale[nodes, {0, 1}, {0, 1000}]^2];
fvec.wts (2000)
fvec.errwts (2000)
(* 2.71313*10^-7 integral estimate 2.71313*10^-7 error estimate
*)First, it doesn't recognize your input as the error function, but just as a function to be integrated numerically. When it does so, probably what is happening is that the calculator starts by evaluating the integrand at a series of points between the upper and lower limits of integration, using something like the trapezoid rule to approximate the integral. Then it cuts the interval in half and computes it again. If they agree (closely enough) it believes it has converged and reports the result. If they disagree, it keeps working on smaller and smaller pieces (maybe changing the spacing to put lots of points where the function is changing rapidly) until it converges. When the interval becomes long enough, the spacing between the points gets wide enough that it misses the hump around $x=0$ completely, which is why it reports such a small number. You might find that if you set $a$ even larger, it returns exactly zero because the value of the integrand is so small at all the points it samples.
$\endgroup$ 2 $\begingroup$We know that $$\lim_{a\to\infty}\left(\int^a_{-a}e^{-x^2}dx\right)^2=\lim_{a\to\infty}\dfrac{1}{4}\pi\left(\operatorname{erf}(a)^2+\operatorname{erf}(-a)^2-2\operatorname{erf}(a)\operatorname{erf}(-a)\right)=\\ \dfrac{\pi}{4}(1+1-(-2))=\pi$$ Hence, I am led to believe that this is a bug in the calculator itself, since the value should be getting closer to $\pi$.
$\endgroup$
