# Evaluating numerically $\int_0^{\infty}e^{-t^2 /100} \sin \pi t$

What is an appropriate method to approximate $$I=\int_0^\infty e^{-t^2 /100} \sin \pi t \ dt?$$ This is for a Physics problem, but in fact I need this in general, as my professor and book taught us nothing about numerical methods for integrals. Also, I found on the internet various techniques for proper ones, but I'm having trouble with making $I$ manageable.

Finally, note that I have indeed come to know $$\int_0^\infty e^{-t^2 /\alpha} \sin \pi t \ dt =\sqrt{\alpha}F\left(\frac{\sqrt{\alpha} }{2}\pi\right),$$ where $F$ is Dawson's integral, however what I'm asking for is approximations.

• The simplest methods you can think of like $I \approx \frac{x_{\rm max}}{n}\sum_{i=1}^nf\left(\frac{x_{\rm max}i}{n}\right)$ is usually more than good enough. For example $x_{\rm max} = 50$ and $n=1000$ gives the result to $0.2\%$ accuracy. Commented Jul 24, 2015 at 19:45

Option 1:

Cut off the domain on $[0,M]$, for some $M$ chosen such that the integral over $[M,\infty)$ is at most half of the desired tolerance. Then use a standard quadrature routine like the trapezoidal rule or Simpson rule on $[0,M]$, with a small enough step size that the error on $[0,M]$ is at most half of the desired tolerance.

Option 2:

Use option 1, but with larger step sizes where the integrand is smaller. This is helpful in your problem, because your integrand is extremely small outside a moderately large interval around $0$ (like $[0,100]$).

Option 3:

Define $F(t)$ to be the even extension of $\sin(\pi t)$ from $[0,\infty)$ to $\mathbb{R}$. (So $F(t)=\sin(\pi t)$ for $t \geq 0$ and $F(t)=\sin(-\pi t)$ for $t<0$.)

Write $$I=\int_0^\infty e^{-t^2/100} \sin(\pi t) dt = \frac{1}{2} \int_{-\infty}^\infty e^{-t^2/100} F(t) dt =5 \int_{-\infty}^\infty e^{-s^2} F(10s) ds$$

Option 4:

Similar to option 3, write

$$I=\frac{\sqrt{50}}{2} \int_{-\infty}^\infty e^{-s^2/2} F(\sqrt{50}s) ds$$

and then recognize this as $\frac{\sqrt{50}}{2} E(F(\sqrt{50}X))$, where $X$ is a standard normal random variable. Then compute this expectation using a Monte Carlo procedure. (This is probably the worst of the four options I've given here.)

Option 5 (perhaps the best, suggested by @uranix):

Write

$$I=10 \int_0^\infty e^{-s^2} F(10s) ds$$

and then use Gauss quadrature directly. As far as I know, this particular form of Gauss quadrature does not have a special name and its parameters are not well-tabulated, but one can use the Golub-Welsch algorithm to generate the nodes and weights directly. The only thing the Golub-Welsch algorithm needs for this case that it can't generate by itself is $\int_0^\infty e^{-x^2} dx$ which is of course $\frac{\sqrt{\pi}}{2}$.

• Your $F(x)$ function has a discontinuous derivative, I expect convergence issues with Gauss-Hermite rule. For oscillatory integral Filon rule might be used instead of trapezoidal or Simpson. Commented Jul 24, 2015 at 20:33
• @uranix Gauss quadrature is convergent for any bounded, continuous integrand (the weight function not included in the "integrand"). The only question is how high of a degree of a polynomial is required in order to achieve a given estimate (in an appropriate metric depending on the weight function) of the integrand by a polynomial. I don't think that is actually so high in this case. (For instance, you can approximate $F$ by a smooth cutoff of $F$, which is equal to $F$ outside $(-2\delta,2\delta)$ but is zero inside $(-\delta,\delta)$, then apply the standard error estimate to the cutoff.)
– Ian
Commented Jul 24, 2015 at 20:44
• Sure, but it may be easier to construct a Gaussian rule for weight $w(x) = e^{-x^2}$ for $[0, \infty)$ interval. Didn't find the corresponding polynomial set, but the procedure is straightforward. Commented Jul 24, 2015 at 22:07
• @uranix That's a fair point (Golub-Welsch isn't that hard to implement, and the all-important normalization integral $\int_0^\infty e^{-x^2} dx$ is well-known.)
– Ian
Commented Jul 25, 2015 at 0:20

Besides numerical integration methods, you can use the asymptotic series for Dawson's integral to get very good approximations. Taking $7$ terms of the series: $$I \approx \dfrac{1}{\pi} + \dfrac{2}{100 \pi^3} + \dfrac{12}{10^4 \pi^5} + \dfrac{120}{10^6 \pi^7} + \dfrac{1680}{10^8 \pi^9} + \dfrac{30240}{10^{10} \pi^{11}} + \dfrac{665280}{10^{12} \pi^{13}}$$ with error approximately $6.2 \times 10^{-15}$.

Since $e^{-t^2/100}$ gets much smaller as $t$ gets larger, you can just reduce the upper bound to a reasonably small number and apply one of the various techniques you've found on the internet for proper integrals.

The fact that $e^{-t^2/100}$ becomes small enough to ignore is strengthened by the fact that the $\sin$ function will make sign changes of the function, therefore some parts of the graph will cancel out.

• If you want to approximate , it suffice to compute $$\int_{0}^{40}e^{-\frac{t^2}{100}}sin(\pi t) dt$$ because $t>35 ,e^{-\frac{t^2}{100}}< 0.00001$ Commented Jul 24, 2015 at 19:46
• @daryakhosrotash It's not quite so simple. You need to do something along the lines of this calculation: $\int_x^\infty e^{-t^2} dt = \int_{x^2}^\infty e^{-t} \frac{1}{2t} dt \leq \frac{1}{2x^2} \int_{x^2}^\infty e^{-t} dt = \frac{1}{2x^2} e^{-x^2}$. Otherwise you can mistakenly get extremely incorrect results; for instance, you can "conclude" that $\int_1^\infty \frac{1}{x} dx$ is finite with a naive argument like that.
– Ian
Commented Jul 24, 2015 at 20:50

One way is to simply use whatever numerical integration technique you like (you can find many online, popular common ones include trapezoid rule and Simpson's rule), with a fixed step size that decreases as $\alpha$ increases (e.g. using derivative-based estimates on error to choose your step-size), and then stop after $t$ becomes large enough that the rest of the integral is negligible compared to your current estimate. For this, a numerical way is to use approximations of the Erf($x$) function, since the tail error in your integral, i.e. the integral from $x$ to $\infty$ is bounded by the Erf($x$) function after suitable variable transformation to bring the exponential into standard form (i.e., removing $\alpha$). The reason why this is a tail error bound is that $|\sin y| \leq 1$.

Dawson's integral can be approximated by using a continued fraction, since it is a hypergeometric function fulfilling the differential equation $F'+x F = 1$. For large values of the parameter $x$ (say $x\geq 10$), a good approximation of the Dawson function is provided by: $$F(x)\approx\frac{x}{1+2x^2}.$$

This looks, to me, like a theta function, which has a functional equation that is very useful for this type of calculation. I first saw this in Bellman's little book on theta functions, which is currently available here: http://www.amazon.com/Brief-Introduction-Theta-Functions-Mathematics/dp/0486492958

After a little searching, I found this: http://press.princeton.edu/chapters/s10_7563.pdf

Here is the relevant stuff (hope I have transcribed it correctly):

Let $$\Theta(z \mid \tau) =\sum\limits_{n=-\infty}^{\infty} e^{\pi i n^2 \tau}e^{2\pi i n z}.$$

Then $$\Theta(z \mid -1/\tau) =\sqrt{\tau \over i}e^{\pi i \tau z^2} \Theta(z\tau \mid \tau).$$ With $x$ real and $\tau = it$, this becomes $$\sum\limits_{n=-\infty}^{\infty} e^{-\pi n^2 /t}e^{2\pi i n x} =t^{1/2}e^{-\pi t x^2}\sum\limits_{n=-\infty}^{\infty} e^{-\pi n^2 t}e^{-2\pi n xt}.$$

Note that for small $t$ the left side converges quickly, and that for large $t$ the right side converges quickly.

In the same spirit as Robert Israel's answer, an asymptotic expansion of Dawson's integral is $$F(x)=\frac{1}{2 x}+\frac{1}{4 x^3}+\frac{3}{8 x^5}+\frac{15}{16 x^7}+\frac{105}{32 x^9}+\frac{945}{64 x^{11}}+\frac{10395}{128 x^{13}}+\frac{135135}{256 x^{15}}+O\left(\frac{1}{x^{16}}\right)$$ For $x=10$, the absolute error is $4.33\times 10^{-14}$ and the relative error is $8.62\times 10^{-13}$.

The general form is $$F(x)=\sum_{k=0}^\infty \frac{(2k-1)!!}{2^{k+1}x^{2k+1}}$$