Computing Gaussian-Weighted Integrals for Small Lengthscales

Question

I am attempting to compute the following integral:

$$ I(l) = \int_{[-1,1]^2}dxdy \sqrt{1-x^2} \sqrt{1-y^2}e^{-\frac{(x-y)^2}{l^2}} $$

I am using Mathematica for symbolic calculations and translating the results into a Fortran code with the help of Mathematica's FortranForm function. This method will be invoked frequently, so it is crucial to have a formula that is both simple and executes very quickly. Numerical integration methods are not suitable for this purpose.

If $ l \gg 1 $ or$ l \sim 1 $, I can easily find a solution using series expansion, ChebyshevT, or Pade approximation. However, none of these approaches works effectively when $l \ll 1 $ because the Gaussian becomes very thin, unlike the nearly constant behavior observed when $ l \gg 1 $. Do you have any suggestions on how to address this issue?

It's also worth mentioning that I will need to compute integrals of the form:

$$ J(x; l) = \int_{-1}^{1}dy \sqrt{1-y^2} e^{-\frac{(x-y)^2}{l^2}} $$

So, numerically computing the integral and fitting some curve on it might not be the best idea.

Thanks in advance, L.

Answer 1

Update: It appears that the first integral has a solution. Taylor expanding the Gaussian leads to an infinite series, each term is easy to integrate and, by luck, the infinite sum can be rewritten into a simple expression using HyperGeometricPFQ functions. Since those functions have existing implementations that are numerically stable, I am totally happy with that.

However, I am still struggling with the second integral $J(x;l)$, and the series does not seem to be converging to a simple expression in that case...

Answer 2

$$I(x;l)=\int_{-1}^{1}\sqrt{(1-y^{2})}e^{-\frac{(x-y)^{2}}{l^{2}}}dy $$

For small $l$ , the integral is dominated by the largest value of the exponential. Thus, when $-1<x<1$ we can ignore all $y$ that are not very near $x$ and write $y=x+t$

$$I(x;l)=\int\sqrt{(1-(x+t)^{2}}e^{-\frac{t^{2}}{l^{2}}}dt\approx\int\left(1-\frac{1}{2}x^{2}+O(t^{2})\right)e^{-\frac{t^{2}}{l^{2}}}dt\approx(1-\frac{1}{2}x^{2})\int_{-\infty}^{\infty}e^{-\frac{t^{2}}{l^{2}}}dt$$

At precisely $x=-1$ we may still use this method, but the peak is is at the edge, so we should use half the value.

When $|x|>1$ we can expand near the larger domain edge, e.g. when $x<-1$, expand near $y=-1$ setting $y\approx-1+t$:

$$I(x;l)=\int\sqrt{(1-(-1+t))(1+(-1+t)}e^{-\frac{(x+1-t)^{2}}{l^{2}}}=\int\sqrt{(2-t)t}e^{-\frac{(x+1)^{2}}{l^{2}}+\frac{(x+1)t}{l^{2}}-\frac{t^{2}}{l^{2}}}dt$$

For small $l$ we can ignore the $\left(\frac{t}{l}\right)^{2}$ term in the exponential, since it kicks in after $|\frac{(x+1)t}{l^{2}}|$ becomes large and the rest of the integral is negligible, so that

$$I(x;l)\approx2e^{-\frac{(x+1)^{2}}{l^{2}}}\int_{0}^{\infty}\sqrt{t}e^{\frac{(x+1)}{l^{2}}t}dt$$ (recall that x+1<0 so the integral converges).

Edit As @ian pointed out, when $x$ is close to the integration endpoints, special care needs to be taken, since the Gaussian integral can't be assumed to extend from $(-\infty,\infty)$. I have ignored this edge case in my answer.