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 |
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.
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.
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.