$\int_{-1}^u\text{exp}(\frac{1}{x^2-1})dx$

Question

I'm trying to compute $\int_{-1}^u\text{exp}(\frac{1}{x^2-1})dx$ where $u\in[-1,1]$.

This is a crucial element of this paper and I need to be able to compute it quickly in Mathematica thousands of times over for a simulation.

I am aware that this post may give some clue as to the answer (i.e. Meijer's $G$-function or Whittaker's $W$ function may come into play) however, because $u$ is not fixed, I cannot use it directly.

If we could even compute $\int_{0}^u\text{exp}(\frac{1}{x^2-1})dx$ where $u\in[0,1]$, because of the symmetry of the bump function, we could figure all the rest out.

Any ideas as to where to start to begin expressing this in terms of another function? Or should I resign myself to numerical integration each time?

Answer 1

This integral can be calculated when $u = 1$. Call it $I$. First let $x = \sin(u)$. The integral becomes: \begin{equation} I = 2\int_{0}^{\frac{\pi}{2}}e^{-(1+\tan^{2}u)}\cos(u) du \end{equation} Changing to $y = \tan(u)$ and introducing the factor $t$ in the exponent yields: \begin{equation} I(t) = 2\int_{0}^{\infty}\frac{e^{-t(1+y^2)}}{(1+y^2)^{3/2}} dy \end{equation} We are seeking $I(1)$. Upon differentiating with respect to $t$ and letting $y^2 = u$ yields: \begin{equation} \frac{dI}{dt} = e^{-t} \int_{0}^{\infty} \frac{e^{-tu}}{\sqrt{u(1+u)}}du \end{equation} This integral can be evaluated in terms of the Modified Bessel Function $K_{0}$ (Gradshteyn 3.364.3), which after integration over $t$ is expressed in terms of both $K_{0}$ and $K_{1}$, giving for $I(1)$: \begin{equation} I = 2+\frac{2}{e}(K_{0}(1)-K_{1}(1)) \end{equation}

Answer 2

Since you need to reference lots of values of this, you want to be able to approximate the function efficiently which points to some form of interpolation. Broadly speaking:

1. Calculate the integral numerically for some fixed number of points (e.g. $u = -1, -0.9, \ldots, 0.9, 1$).

2. Use those points to fit a smooth curve through them, such as a polynomial.

3. Create a function that returns the value of the fitted curve for any given input.

In Mathematica, you should be able to use the InterpolatingPolynomial function to achieve step 2. Because the function is so smooth, you probably don't need to use too many points to get a decent interpolation.

Answer 3

If your main concern is accurate numerical calculation (rather than finding a closed form expression), then I think the ideal tool for this problem would be to use Chebyshev function approximation.

First, you calculate the coefficients for the Chebyshev series of the integrand. For this particular function, you need a few hundred of them to achieve double precision accuracy. Then there's a simple formula to calculate the Chebyshev coefficients of the integral. With those coefficients, you can quickly calculate the integral to machine precision at any point in the domain $[-1,1]$.

In python, the numpy library can do all this in a single line of code:

f = np.polynomial.Chebyshev.interpolate(lambda x: np.exp(1/(x**2 -1)), 256).integ()

For details on the theory, refer to the work of Nick Trefethen and collaborators (the Chebfun package). But, briefly, the Chebyshev series is a way of representing a function on a finite interval by expanding it as a sum of orthogonal polynomials: $$ e^{\frac{1}{x^2-1}} = \sum_{n=0}^\infty a_n T_n(x) \quad \forall x \in [-1,1] $$

$$T_n(\cos(x)) := \cos(n x)$$

$$ a_n := \frac{\int_{-1}^{1} e^{\frac{1}{x^2-1}} T_n(x) \frac{dx}{\sqrt{1-x^2}}} {\int_{-1}^{1} T_n(x)^2 \frac{dx}{\sqrt{1-x^2}}} $$ Then with those coefficients you can use the properties of the Chebyshev polynomials to integrate it: $$ \int_0^x e^{\frac{1}{t^2-1}} dt = \sum_{n=1}^\infty \frac{a_{n-1} - a_{n+1}}{2n} T_n(x) \quad \forall x \in [-1,1] $$ In practice, we compute enough coefficients to achieve a desired accuracy. There are fast algorithms for evaluating the series at a given point.

Answer 4

After graphing this function, it looks like some sort of logistic function or error function. So a good approximation is, by trial and error, $\frac{\mathrm{erf}(x)}3+\frac 29=I(x)$, where $I(x)$ is the function in question. The $\frac 29$ comes from the fact that $I(0)\approx0.222$. Another simpler but worse approximation is $\frac1{1+e^{-x}}-\frac5{18}$, where the fraction here is $-\frac12+\frac29$ so that the y-intercept of the logistic function is close to $I(0)$.