A Numerical Curiosity

Some years ago, with my HP 35 hand calculator and too much time on my hands, I decided to numerically compute ∫ e−x2dx from −∞ to ∞. I couldn’t remember if it was √π or half that. I first tried a Δx = 0.1 expecting that the answer would make it clear despite the large integration error which would pessimistically be about Δx2. To my astonishment I got 10 digits of accuracy. The answer is √π. This code yields nearly 15 digits of accuracy—within rounding error. This Algol 68 program provides 44 digits. Here is a compiler that supports these precisions. What can account for this astonishing accuracy?

Trying Δx = 1.0 yields only 3 digits—there is no lurking identity. Initializing x to −10.3 yields 44 digits again. Increasing the integration limits to [−35, 35] yields 428 correct digits! This is not what you learn in either theoretical or practical numerical analysis!

This gnawed on me for perhaps a year until I thought of an explanation of sorts. In the following Σ denotes summation with n running over all integers and ∫ denotes integration with real x from −∞ to ∞. δ(x) is that illegitimate function which integrates to 1 over any interval including 0 and to 0 for all other intervals.

Fix some positive value for Δx. In the following “a ≃ b” means that lim(a − b) → 0 as Δx → 0.
1∫ f(x)dx ≃ Σf(Δx • n)Δx Numerical integration
2f(y) = ∫ f(x)δ(x − y) dxSampling f at y
3p(x) = ΔxΣδ(x − n • Δx) Picket fence function (or Dirac comb)
4∫ f(x)dx ≃ ∫ (p(x)•f(x))dx Numerical integration of f as inner product of p and f
For a complex valued function f on the real line let Tf be the Fourier transform of f.
(Tf)(y) = ∫ f(x)e2πixydx and f(y) = ∫ (Tf)(x)e−2πixydx
T(λx e−x2) = λx (√π)e−(πx)2 Authority, proof

Now we consider some simple Hilbert space lore. These complex functions on the real line form a vector space over the complex numbers. The expression (f, g) = ∫ f(x)•(g(x))*dx serves as an inner product on this vector space. (x* = complex conjugate of x.) A significant theorem is that the Fourier transform preserves this inner product:
(f, g) = (Tf, Tg).
As a linear transformation on the vector space, T preserves vector lengths and angles between vectors; it is a rigid transformation. Our original numerical integration amounts to (p, λx e−x2). According to the theorem this is (Tp, T(λx e−x2)) = (Tp, λx (√π)e−(πx)2). But what might Tp be? Recall that in equation 3 above p depends on Δx. Tp is another picket fence function whose spacing is 1/Δx.
(Tp)(x) = (1/Δx)Σδ(x − n/Δx).
Our expression thus becomes
(Tp, λx (√π)e−(πx)2) = (λx(1/Δx)Σδ(x − n/Δx), λx (√π)e−(πx)2)
= Σ(1/Δx)∫ δ(x − n/Δx) • (√π)e−(πx)2dx
= Σ(1/Δx) (√π)e−(π(n/Δx))2

This last sum addresses the mystery. The term for n = 0 is: (1/Δx) (√π)e−02. Ignoring (1/Δx) for a bit we can see that this is the correct answer and that the remaining terms are the truncation error in the original numerical integration.
The next terms, for n = ± 1 are each:
(1/Δx) (√π)e−(π(1/Δx))2.
When Δx = 0.1, (π(1/Δx))2 is about 987 and e−987 is so small that the first term dominates.

You may have noticed that pesky (1/Δx). I think that shows up because p is doubly illegitimate. The limit of p as Δx goes to 0 is the constant function λx 1. Neither p nor its limit belongs to Hilbert space.

Nonetheless we have a numerical theory.
Lets try Δx = 1/2. The error term is 2(√π)e−(2π)2 = 2.5*10−17. The numerical integration yields an error of 9*10−17. I am sort of satisfied with this explanation.

Worrying the (1/Δx)

We could say that ∫ f(x)dx = (λx 1, f(x)) except that λx 1 isn’t in Hilbert space since is it not square integrable. No one has succeeded in stretching Hilbert enough to include λx 1. We are skating on thin ice here and the ice gives away before we get to these functions. Heaviside was adept, however, at bulling thru such calculations and guessing how to fix the formulas so that they were useful for engineering. They were a means of discovery.

I think that the function δ could be replaced by a series of ever narrower functions in some limiting process. Then the function p would have to be replaced by some double limit process.

They know about this.

This scheme does not extend in the obvious way to integration limits other than −∞ to ∞. The routine di(x, dt) here computes the integral from x to ∞ using dt. A small improvement in dt changes the answer drastically.