Approximating specific integrals via sums: Tight error predictions

Question

Dear math enthusiasts,

I am facing a particular problem where I am looking at integrals of the form $$I = \int_{-\infty}^\infty p(t) {\rm e}^{-t^2} {\rm d}t,$$ where $p(t)$ are certain polynomials. These are either even or odd symmetric, the even symmetric ones being the more interesting case since for odd ones, the integral is zero (in fact, they are auto- and cross-products of Hermite polynomials but I think this detail is not relevant). For now, we can just consider the simplest example $p(t)=t^2$.

The reason why I am looking at the integrals is that I actually have sums of the form $$S(t) = \sum_{n=-\infty}^\infty t_0 p(t-nt_0) {\rm e}^{-(t-nt_0)^2}$$ which I want to quantify. For $t_0$ small enough, these sums are very close to $I$ for any $t$. I need to quantify how close, i.e., I am interested in $$\max_t |S(t)-I|.$$

So, what I did was to interpret $S$ as a quadrature of the integral $I$ and use the residual formulas for quadratures. Since it is a linear quadrature, the standard results predict a residual of the order $|I-S|< {\rm const} \cdot t_0^2 \cdot \max|f''(t)|$, where $f(t) = p(t){\rm e}^{-t^2}$ (which gives $\max|f''(t)| = 2$ for $p(t)=t^2$). In other words, the error should decay quadratically with $t_0$. All this is not surprising and well within what I expected.

Until I tried it and realized that empirically, the error decays much much faster with $t_0$ than this pessimistic bound predicts. Here is an example:

In this example, I computed $S$ for $p(t)=t^2$, varying the grid spacing $t_0$. I plot $\max_t|I-S(t)|/I$ where the exact value $I$ is equal to $\sqrt{\pi}/2$.

Obviously, $t_0=1$ is too coarse but as I make it finer, the error goes to zero very rapidly (Note that the plot is doubly logarithmic!). In fact, around $t=0.35$ it reaches the numerical accuracy of my double floating point but I would expect the exponentially decaying trend to continue. The predicted upper bound is shown in the dashed line (it is a line with slope 2 due to the double logarithmic plot).

So here is my question: Can I make tighter predictions of $|S(t)-I|$ as $t_0\rightarrow 0$? I know it is easy to construct examples where the residual formula is basically tight, so it must have to do with the particular function I am integrating, especially the ${\rm e}^{-t^2}$ term. I keep getting reminded of Gaussian tails (erfc functions of some sort) but I cannot put my hands on how I could get there. Just before posting I stumbled upon Euler-Maclaurin, but it confuses me as well since what I read about it talks about finite sums (and derivatives being evaluated at the borders) while mine seems infinite (and everything becomes zero far enough from $t=0$).

It looks like a standard result and I would not be surprised about a very simple answer that I was just not seeing.

Any hint is appreciated, many thanks in advance!

edit: Thanks to user14717 I got what I needed. For someone stumbling across a similar problem, here is what worked: Theorem 5.1 in [*] says the following: Let $a>0$ such that the function $w(t)$ to be integrated is analytic in the string $|{\rm Im}(t)|<a$ and decays to zero uniformly as $|t|\rightarrow \infty$. Then: $$|I-S| \leq 2\sqrt{\pi} \frac{M}{{\rm e}^{2\pi a/t_0}-1},$$ where $M$ is a constant satisfying $\int |w(t+ib)| {\rm d}t \leq M$ for all $b \in (-a,a)$. If we apply this for $w(t) = {\rm e}^{-t^2}$ (i.e., my $p(t)=1$), we obtain $M={\rm e}^{a^2}$ as the best $M$ for a given $a$. This gives the family of bounds $$|I-S| \leq 2\sqrt{\pi} \frac{{\rm e}^{a^2}}{{\rm e}^{2\pi a/t_0}-1},$$ which is valid for any $a>0$. To obtain the tightest bound, we need to minimize over $a$. This is not possible analytically, however, for $t_0<1$, the value $a=\frac{\pi}{t_0}$ is very close to the optimum. Inserting it gives the bound $$|I-S| \leq \frac{2\sqrt{\pi}}{{\rm e}^{(\pi/t_0)^2}-{\rm e}^{-(\pi/t_0)^2}} \approx 2\sqrt{\pi} {\rm e}^{-(\pi/t_0)^2}.$$

And now, let us plot it:

Adapting this to any $p(t)$ should be a breeze now.

I could barely be happier! :)

Thanks so much!

[*] http://epubs.siam.org/doi/pdf/10.1137/130932132

edit2: For an even symmetric polynomial $p(t)$ it is then very easy to show that $$|S-I|\leq \frac{2{\rm e}^{a^2} h(a)}{{\rm e}^{2\pi a/t_0}-1},$$ where $h(a)$ is a polynomial of same degree as $p(t)$. For $a=\pi/t_0$, this gives $$|S-I|\leq \frac{2{\rm e}^{\pi^2/t_0^2}}{{\rm e}^{2\pi^2/t_0^2}-1}h(\pi/t_0) \approx 2{\rm e}^{-\pi^2/t_0^2}h(\pi/t_0),$$

i.e., still exponential convergence, as expected.

In fact, for a degree $2k$ polynomial $p(t) = \sum_n \alpha_n t^{2(k-n)}$, an explicit form of $h(t)$ is given by $$h(t) = \sqrt{\pi} \sum_n \sum_\ell |\alpha_n|\frac{k! (2\ell)!}{4^\ell (\ell!)^2 (k-\ell)!} b^{2(k-n-\ell)},$$ though it doesn't matter much since the relevant part is the exponential convergence.

Answer 1

It is a standard result. What you are observing is a well-known property of the Euler-Maclaurin expansion \begin{align*} h\frac{f(a) + f(b)}{2} + h\sum_{k = 1}^{n-1} f(a+kh) = \int_{a}^{b} f(x) \, \mathrm{d}x + \sum_{k = 1}^{\infty} \frac{ h^{2k} }{(2k)!} B_{2k}(f^{(2k-1)}(b) - f^{(2k-1)}(a)) \end{align*} Now, you are worried that this doesn't apply since your integrand is infinite. But in finite precision, $\exp(-t^2) = \epsilon$ which means that in practice your integration can only be performed over the range $[-\sqrt{-\log(\epsilon)}, \sqrt{-\log(\epsilon)}]$. In double precision, this is about $[-6, 6]$.

Now, all the derivatives of $f$ are also bell-shaped, so as $a\to -\infty$, $f^{(2k-1)}(a) \to 0$ very fast. Assuming that $B_{2k}f^{(2k-1)}(a)/(2k)!$ is tiny, we can state that your trapezoidal sum converges faster than any power of $h$. Let's examine this claim in more detail: Using the representation of the Bernoulli numbers \begin{align*} B_{2k} = (-1)^{k+1}\frac{2(2k)!}{(2\pi)^{2k}}\zeta(2k) \end{align*} we can write the error as \begin{align*} E(h) := \sum_{k = 1}^{\infty} (-1)^{k+1}h^{2k}\frac{2}{(2\pi)^{2k}}\zeta(2k)(f^{(2k-1)}(b) - f^{(2k-1)}(a)) \end{align*} In your case, it makes sense to make $a = -b$, and since we assume that $f$ is even then all odd-order derivatives of $f$ are odd. Then \begin{align*} E(h) = 4\sum_{k = 1}^{\infty} (-1)^{k+1}\left(\frac{h}{2\pi}\right)^{2k}\zeta(2k)f^{(2k-1)}(b) \end{align*}

Now, $1 < \zeta(2k) < 2 \forall k \in \mathbb{N}$ so \begin{align*} |E(h)| \le 8\sum_{k = 1}^{\infty} \left(\frac{h}{2\pi}\right)^{2k}|f^{(2k-1)}(b)| \end{align*} All that we need is for the error to be smaller than the machine epsilon to get the result, as $f(0) \sim \mathcal{O}(1)$ and hence any corrections smaller than the unit roundoff are not observable.

I have an idea for getting this argument a bit more rigorous, but it needs some more TLC (or a counterexample):

Assume $f(x) = (x^{2n} + \cdots)\exp(-x^2)$ where $\cdots$ are terms with degree lower than $2n$. Then $f^{(2k-1)}(x) = (x^{2n+2k-1} +\cdots)\exp(-x^2)$ and hence for large $x$, $|f^{(2k-1)}(x)| \le C|x^{2k-1}|x^{2n}\exp(-x^2)$, so that \begin{align*} |E(h)| \le 8Cb^{2n-1}\exp(-b^2)\sum_{k = 1}^{\infty} \left(\frac{hb}{2\pi}\right)^{2k} = 8Cb^{2n-1}\exp(-b^2) \frac{(hb/2\pi)^2}{1-(hb/2\pi)^2} \end{align*} assuming $bh/2\pi < 1$. Now just choose $b$ large enough that $|E(h)| < \epsilon$ and the result is proved.

To prove the result without recourse to arguments about finite precision requires complex analysis. This is discussed by Trefenthen here.