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

Example of the rapid decay of the approximation error for p(t)=t^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:

New upper bound

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.

up vote 3 down vote accepted

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.

  • Great, this is helpful! Numerical issues let aside I could simple choose any fixed $[a,b]$ big enough, in a second step even looking at the limit $[a,b]\rightarrow[-\infty,\infty]$. Then I have $f^{(2k-1)}(t) = O(t^{2k+1}e^{-t^2})$ and I can bound $B_{2k}/(2k)!$ via Stirling, right (since it is constant, the $e^{-t^2}$ should take care of it all I guess)? Looks good, tomorrow I'll have a look at what I can do about the error term. I'd really like to get the the correct scaling law (without killing myself on the way). – Florian Apr 10 '17 at 17:03
  • Great, the addition helps a lot. Just one more clarification how to get from the last inequality to "faster than any power of $h$": I would assume the argument is that if I were to truncate the sum to any finite number of terms $k=1...K$ then the limit would surely be zero (since $e^{-b^2}h^{2k}$ always goes to zero) - but I cannot say this for the infinite sum (since even if all terms are zero the infinite sum could have any value)? – Florian Apr 10 '17 at 18:54
  • Well I guess my $C_k$ are not - if we keep forming derivatives we essentially get Hermite polynomials in front of the exp and their coefficients grow with $k$, if I am not mistaken. But that's fine, I think I can work with this answer. Thank you very much for taking the effort to reply in such detail, greatly appreciated! – Florian Apr 10 '17 at 19:27
  • Thanks for the additional edit! The only problem I see is that I cannot choose $b$ "large enough" since we need $bh/2\pi<1$ and $h$ is not necessarily very small. What I'm really after is a point $h$ below which $|R(t)-I|$ is practically zero, so I need a fast decay in $h$ ("a little smaller than one" seems to work, I'd like to quantify that). It seems it's almost there, the only problem being that the $b^{2k}$ is still inside the sum which means I cannot let $b$ grow too much withouth constraining $h$ (and, rigorously speaking I really want to say something about $\lim b\rightarrow \infty$). – Florian Apr 11 '17 at 8:25
  • In double precision, $b \sim 6$, so we only need $h < 1$. I don't see why that is a big restriction. As to the limit $b \to \infty$, this mode of analysis isn't going to help much, you'll need another strategy. – user14717 Apr 11 '17 at 13:24

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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