Skip to article frontmatterSkip to article content

Improper integrals

When the interval of integration or the integrand itself is unbounded, we say an integral is improper. Improper integrals present particular challenges to numerical computation.

9.7.1Infinite interval

When the integration domain is (,)(-\infty,\infty), the integrand has to decay as x±x \to \pm \infty in order for the improper integral to be finite. This fact brings up the possibility of truncating the domain:

f(x)dxMMf(x)dx.\int_{-\infty}^{\infty} f(x)\, dx \approx \int_{-M}^{M} f(x)\, dx.

This integral can be discretized finitely by, say, the trapezoid formula, or an adaptive integrator. However, this approach can be inefficient.

9.7.2Double exponential transformation

In order to do better than direct truncation, we want to encourage the function to decay faster. In practice this means a change of variable, x(t)x(t). If x(t)|x(t)| grows rapidly as t|t| \to \infty, then f(x(t))f(x(t)) will decay more rapidly in tt than in xx.

One way to accomplish this feat is to use

x(t)=sinh(sinht). x(t) = \sinh\left( \sinh t \right).

Noting the asymptotic behavior as t±t \rightarrow \pm\infty that

sinh(t)12et, \left| \sinh(t) \right| \sim \frac{1}{2} e^{ |t| },

we find that in the same limits,

x(t)±12exp(12et).x(t) \approx \pm \frac{1}{2} \exp\left( \frac{1}{2} e^{ |t| } \right).

Thus, (9.7.4) is often referred to as a double exponential transformation.

By the chain rule,

f(x)dx=f(x(t))dxdtdt=f(x(t))cosh(sinht)coshtdt.\begin{split} \int_{-\infty}^\infty f(x)\, dx &= \int_{-\infty}^\infty f(x(t))\frac{dx}{dt}\, dt \\ &= \int_{-\infty}^\infty f(x(t))\, \cosh\left( \sinh t \right) \cosh t \, dt. \end{split}

The exponential terms introduced by the chain rule grow double exponentially, but the more rapid decay of ff in the new variable more than makes up for this.

Function 9.7.1 implements double exponential integration by applying the adaptive integrator Function 5.7.1 to (9.7.7). It truncates the interval to MtM-M\le t \le M by increasing MM until the integrand is too small to matter relative to the error tolerance.

9.7.3Integrand singularity

If ff asymptotically approaches infinity as xx approaches an integration endpoint, its exact integral may or may not be finite. If ff is integrable, then the part of the integration interval near the singularity needs to be more finely resolved than the rest of it.

Let’s consider

01f(x)dx,\int_0^1 f(x)\,dx,

where ff and/or a derivative of ff is unbounded at the left endpoint, zero. The change of variable

x(t)=21+exp(2sinht) x(t) = \frac{2}{1+\exp(2 \sinh t)}

satisfies x(0)=1x(0)=1 and x0+x\to 0^+ as tt\to \infty, thereby transforming the integration interval to t(0,)t\in(0,\infty) and placing the singularity at infinity. The chain rule implies

01f(x)dx=0f(x(t))dxdtdt=0f(x(t))coshtcosh(sinht)2dt.\begin{split} \int_{0}^1 f(x)\, dx &= \int_{0}^\infty f(x(t)) \frac{dx}{dt}\, dt \\ &= \int_{0}^\infty f(x(t)) \frac{\cosh t}{\cosh(\sinh t)^2} \, dt. \end{split}

Now the growth of ff and cosht\cosh t together are counteracted by the double exponential denominator, allowing easy truncation of (9.7.13). This variable transformation is paired with adaptive integration in Function 9.7.2.

Double exponential integration is an effective general-purpose technique for improper integrals that usually outperforms interval truncation in the original variable. There are specialized methods tailored to specific singularity types that can best it, but those require more analytical work to use properly.

9.7.4Exercises

  1. ⌨ Use Function 9.7.1 to estimate the given integral with error tolerances 103,106,109,101210^{-3},10^{-6},10^{-9},10^{-12}. For each result, show the actual error and the number of nodes used.

    (a) 11+x2+x4dx=π3\displaystyle\int_{-\infty}^\infty \dfrac{1}{1+x^2+x^4}\, dx = \dfrac{\pi}{\sqrt{3}}

    (b) ex2cos(x)dx=e1/4π\displaystyle\int_{-\infty}^\infty e^{-x^2}\cos(x)\, dx = e^{-1/4}\sqrt{\pi}

    (c) (1+x2)2/3dx=πΓ(1/6)Γ(2/3)\displaystyle\int_{-\infty}^\infty (1+x^2)^{-2/3}\, dx = \dfrac{\sqrt{\pi}\,\Gamma(1/6)}{\Gamma(2/3)} (use gamma() for Γ()\Gamma())

  2. ⌨ Use Function 9.7.2 to estimate the given integral, possibly after rewriting the integral into the form (9.7.11) with a left-endpoint singularity. Use error tolerances 103,106,109,101210^{-3},10^{-6},10^{-9},10^{-12}, and for each result, show the actual error and the number of nodes used.

    (a) 01(logx)2dx=2\displaystyle\int_{0}^1 (\log x)^2\, dx = 2

    (b) 0π/4tan(x)dx=π2\displaystyle\int_{0}^{\pi/4} \sqrt{\tan(x)}\, dx = \dfrac{\pi}{\sqrt{2}}

    (c) 0111x2dx=π2\displaystyle\int_{0}^1 \frac{1}{\sqrt{1-x^2}}\, dx = \dfrac{\pi}{2}

  3. For integration on a semi-infinite interval such as x[0,)x\in [0,\infty), another double exponential transformation is useful: x(t)=exp(sinht)x(t)=\exp\left( \sinh t \right).

    (a) ✍ Show that t(,)t\in(-\infty,\infty) is mapped to x(0,)x\in (0,\infty).

    (b) ✍ Derive an analog of (9.7.7) for the chain rule on 0f(x)dx\int_0^\infty f(x)\,dx.

    (c) ✍ Show that truncation of tt to [M,M][-M,M] will truncate xx to [1/μ,μ][1/\mu,\mu] for some positive μ.

    (d) ⌨ Write a function intsemi(f,tol) for the semi-infinite integration problem. Test it on the integral

    0exxdx=π.\displaystyle\int_0^\infty \frac{e^{-x}}{\sqrt{x}}\,dx = \sqrt{\pi}.