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 |

2 | f(y) = ∫ f(x)δ(x − y) dx | Sampling f at y |

3 | p(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 |

(Tf)(y) = ∫ f(x)e

T(λx e

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)2}dx

= Σ(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.

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.