2
$\begingroup$

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$
  • $\begingroup$ Don't know, so this isn't an answer, but sure does look like a bug. For abs(a) big enough, the calculator ought to just approximate it as + or - infinity, and use the asymptotic values of the integrated gaussian, related to the well-known "erf" function. $\endgroup$ – DarenW Apr 8 '14 at 23:06
  • $\begingroup$ I tried this on my TI-84 plus and while I got a different value, 7.36107732x$10^{-14}$, the value is obviously not $\Pi$. This is a very interesting bug. $\endgroup$ – TheBluegrassMathematician Apr 8 '14 at 23:24
  • $\begingroup$ Related discussion. $\endgroup$ – user127096 Apr 13 '14 at 22:22
1
$\begingroup$

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$
  • $\begingroup$ The TI doesn't really know how to do integrals exactly, but instead probably uses this approach. You can try doing $\int_0^1\text{rand}\ dx$. It will spit out different values ~$0.5$. Sometimes it will take a while to think it converged to a value, and sometimes it will determine that it never converges. $\endgroup$ – user137794 Apr 8 '14 at 23:33
  • 1
    $\begingroup$ Search results suggest that TI-84+ uses Gauss-Kronrod quadrature, but I could not find the details such as the number of nodes. $\endgroup$ – user127096 Apr 13 '14 at 22:23
1
$\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
*)
$\endgroup$
  • $\begingroup$ If you have a reference for the algorithm used by TI-86, then please add it to your answer as this type of information can be hard to track down. Kind regards $\endgroup$ – Carl Christian Mar 15 at 9:31
  • 1
    $\begingroup$ @CarlChristian I added the reference for the TI-84. The Guidebooks can be found here, but the one for the TI-86 does not seem to indicate the method(s) used. $\endgroup$ – Michael E2 Mar 15 at 14:15
  • 1
    $\begingroup$ @CarlChristian I found this resource at TI.com (it downloads a PDF). It indicates that the TI-86 uses a 7-node rule, which would not give the same results as OP observes (or that I have observed) on the TI-84. I don't have a TI-86 to test. $\endgroup$ – Michael E2 Mar 15 at 14:29
  • $\begingroup$ Thank you! You have already given me more information than I could hope for. $\endgroup$ – Carl Christian Mar 15 at 15:16
0
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.