Improper Integrals I: Infinite bound

5.3.1. Improper Integrals I: Infinite bound#

When infinity is involved in either the integral interval or the integrand, we have an improper integral. In this section, we deal with the former and the latter type of improper integral is discussed in the next section.

If the upper and/or lower limit envolves \(\infty\), the methods of peicewise intetration won’t work since infinite segments must be added up.

For example, consider

\[ \int_1^\infty x^{-2}\, dx\,. \]

Despite the infinite interval this integral converges to \(1\). One might hope that the infinity can be approximated by some large number \(L\). Then, the approxmated integral is

\[ \int_1^L x^{-2}\, dx = 1 - \frac{1}{L} \]

If we want to get four digits of accuracy, \(L\) must be larger than \(10^{4}\). If step size \(h=10^{-2}\) is used, we have sum up \(10^6\) segments.

Exercise 3.2.1: Evaluate the integral using the trapezoidal rule.

The output of the above exercise suggests that this approach is feasible for this simple integral. However, if we need a higher accuracy, the accumulation of round-off error may become an issue. Furthermore, if the integrand is more complicated, it is difficult to find an appropriate upper bound.

There are better strategies. In the following, two approaches are introduced.

5.3.1.1. Method of variable transformation#

One way to overcome the above problem is to split the integral to two parts

\[ \int_a^\infty f(x)\, dx = \int_a^b f(x)\, dx + \int_b^\infty f(x)\,dx \]

where \(b\) is a constant to be chosen. The first term in the right hand side can be integrated by the trapezoidal or the simpson rule. The second term needs to be transformed to a numerically computable form by introducing a new variable \(y=\displaystyle\frac{1}{x}\). Then, the integral we need to compute is

\[ \int_b^\infty f(x)\, dx = \int_0^{1/b} \frac{1}{y^2}\, f\left (\frac{1}{y} \right )\, dy \]

Now the infinity disappeared. However, the new form is not necessarily trouble free since the integrand may not be defined at the lower bound (divided-by-zero error). If we can evaluate \(\displaystyle\lim_{y \rightarrow 0} \frac{1}{y^2} f\left (\frac{1}{y} \right )\) analytically by hand, then standard methods such as the Simpson method works if an appropriate \(b\) is chosen.

Example 3.2.1 Evcaluate \(\displaystyle\int_0^\infty \frac{1}{x^2+1} dx\). The exact answer is known to be \(\pi/2\).

\[ \int_0^\infty \frac{1}{x^2+1} dx = \int_0^b \frac{1}{x^2+1} dx + \int_0^{1/b} \frac{1}{y^2}\frac{1}{\frac{1}{y^2}+1} dy = \int_0^b \frac{1}{x^2+1} dx + \int_0^{1/b} \frac{1}{y^2+1} dy \]

Now we conveniently choose \(b=1\). Then, the two seprate integrals become ideitical. Hence,

\[ \int_0^\infty \frac{1}{x^2+1} dx = \int_0^1 \frac{2}{x^2+1} dx \]

Now, we integrate this using the Simson’s rule.

import numpy as np
import scipy.integrate as integrate
N=100
x = np.linspace(0,1,N+1)

f = 2./(x**2+1)
s=integrate.simpson(f,x=x)
error=s-np.pi/2
print("simpson =",s)
print("  exact =",np.pi/2.)
print("  error =",error)

simpson = 1.5707963267948768
  exact = 1.5707963267948966
  error = -1.9761969838327786e-14

We obtain 14 digits of precision.

There is another variable transformation for \(\int_0^\infty f(x) dx\). Using \(x = (1-y)/y\) and its inverse \(y = 1/(x+1)\)

\[ \int_0^\infty f(x) \, dx = \int_0^1 \frac{1}{y^2} f\left(\frac{1-y}{y}\right)\, dy \]

This is not necessary a better approach than the previous one. One advantage is that we don’t have to devide the interval. In fact, many canned routines such as scipy.integrate.quad uses this method. So, you may want to try this approach first. If not satisfactory, use the previous one.

Exercise 3.2.2: Solve Exercise 3.2.1 using this method. Is there any improvment in the accuracy?

5.3.1.2. Gaussian Quadratures#

There amazing methods, known as Gaussian quadrature. Some of them are quite useful for physics as shown in Appication section. In the following, we introduce two of them Gauss-Laggeurre quadrature and Gauss–Hermite quadrature.

5.3.1.2.1. Gauss-Lagguerre quadrature#

the Gaussian-Laguerre Quadrature approximate the following form of integral as a simple weighted sum:

\[ \int_0^\infty f(x) e^{-x} \, dx \approx \sum_{i=1}^{n} w_i f(x_i) \]

where the bscissa \(x_i\) is the \(i\)-th root of Laguerre polynomial \(L_n(x)\) and the weight \(w_i\) is given by

\[ w_i = \frac{x_i}{(n+1)^2[L_{n+1}(x_i)]^2} \]

Amazingly, \(n\) can be as small as 4 to get a resonable approximation despite that the integral interval is infinite. It is almost a miracle that the sum of several terms approximate the improper integral. Now, a question is if this method can be used in the absence of the exponential function. Suppose that we want to evalute \(\int_0^\infty g(x)\, dx\). Can we do the following?

\[ \int_0^\infty g(x)\, dx = \int_0^\infty g(x) e^{x} e^{-x}\, dx \approx \sum_{i=1}^{n} w_i g(x_i) e^{x_i} \]

It looks cheating. It is known to work in some cases but fails in many other cases. There is no guarntee.

The following code generates the abscissas and the weights for any ginve \(n\) with 20 digits of precision. See Chapter xxx for root finding methods and Chapter yyy for special functions.

# We use sympy package
from sympy import symbols, poly, laguerre

# this function computes abscissas and weights for ginve n
def laguerre_weights_roots(n):
    x = symbols("x")
    roots = poly(laguerre(n, x)).all_roots()
    x_i = [rt.evalf(20) for rt in roots]
    w_i = [(rt / ((n + 1) * laguerre(n + 1, rt)) ** 2).evalf(20) for rt in roots]
    return x_i, w_i

# set the number of abscissas
n=4

# compute all abscissas and weights
[x,w]=laguerre_weights_roots(n)

# print out the results
print(" {0:^3} {1:^25} {2:^25}".format("i","x_i","w_i"))
for i in range(n):
    print("{0:3d} {1:25.19e} {2:25.19e}".format(i,x[i],w[i]))

  i             x_i                       w_i           
3.2254768961939231180e-1  6.0315410434163360164e-1
1.7457611011583465757e+0  3.5741869243779968664e-1
4.5366202969211279833e+0  3.8887908515005384272e-2
9.3950709123011331292e+0  5.3929470556132745010e-4

print(x)
import numpy as np
y=np.array([x[0], x[1], x[2]],dtype=float)
print(y)

print(np.sin(y))

[0.32254768961939231180, 1.7457611011583465757, 4.5366202969211279833, 9.3950709123011331292]
[0.32254769 1.7457611  4.5366203 ]
[ 0.31698389  0.98473267 -0.98459241]

Example 3.2.2 Let us evaluate \(\int_0^\infty \sin(b x) e^{-a x}\, dx\). The exact anser is \(\frac{b}{a^2 + b^2}\). We try \(a=b=1\). Then, the answer should be \(0.5\).

import numpy as np

# Here we use the table of the abscissas and weights
# instread of calculating them.
# Smetime this is mor convenient.
x=np.array([1.7027963230510100e-1, 9.0370177679937991e-1,
            2.2510866298661307,    4.2667001702876588,
            7.0459054023934657,    1.0758516010180995e+1,
            1.5740678641278005e+1, 2.2863131736889264e+1])

w=np.array([3.6918858934163753e-1, 4.1878678081434296e-1,
            1.7579498663717181e-1, 3.3343492261215652e-2,
            2.7945362352256725e-3, 9.0765087733582131e-5,
            8.4857467162725315e-7, 1.0480011748715104e-9])

gauss=(w*np.sin(x)).sum()
exact=1/2.

print(" Exact={0:18.12e}\n Gauss={1:18.12e}\n Error={2:18.12e}".format(exact,  gauss, abs(exact-gauss)))

 Exact=5.000000000000e-01
 Gauss=4.999877537353e-01
 Error=1.224626470031e-05

Example 3.2.3

Considering we used only 8 points, five digits agreement is quite good.

Exercise: Try \(n=4\) and find how good or bad it is.

5.3.1.2.2. Gauss-Hermite quadrature#

Another popular Gaussian quadrature is Gauss-Hermite quadrature which approximates the following improper integral:

\[ \int_{-\infty}^\infty f(x) e^{-x^2}\, dx \approx \sum_{i=1}^n w_i \left[f(x_i)+f(-x_i) \right] \]

where the bscissa \(x_i\) is the \(i\)-th root of Hermite polynomial \(H_n(x)\) and the weight \(w_i\) is given by

\[ w_i = \frac{2^{n-1} n! \sqrt{\pi}}{n^2[H_{n-1}(x_i)]^2} \]

The following code generates the abscissas and the weights for any ginve with 20 digits of precision. See Chapter xxx for root finding methods and Chapter yyy for special functions.

from sympy import *

def hermite_weights_roots(n):
    x = symbols('x')
    roots = poly(hermite(n, x)).all_roots()
    x_i = [rt.evalf(20) for rt in roots]
    w_i = [(2**(n-1) * factorial(n) * sqrt(pi) / (n * hermite(n-1, rt))**2).evalf(20) for rt in roots]
    return x_i, w_i

n=8
[x,w]=hermite_weights_roots(n)

print(" {0:^3} {1:^25} {2:^25}".format("i","x_i","w_i"))
for i in range(n):
    print("{0:3d} {1:25.19e} {2:25.19e}".format(i,x[i],w[i]))

  i             x_i                       w_i           
-2.9306374202572440192e+0  1.9960407221136761921e-4
-1.9816567566958429259e+0  1.7077983007413475456e-2
-1.1571937124467801947e+0  2.0780232581489187954e-1
-3.8118699020732211685e-1  6.6114701255824129103e-1
3.8118699020732211685e-1  6.6114701255824129103e-1
1.1571937124467801947e+0  2.0780232581489187954e-1
1.9816567566958429259e+0  1.7077983007413475456e-2
2.9306374202572440192e+0  1.9960407221136761921e-4

Example 3.2.2 Let us evaluate \(\int_{-\infty}^\infty \cos(x) e^{-x^2}\, dx\). The exact answer is \(\sqrt{\pi}e^{-1/4}\).

import numpy as np

# We use 8-point method.
x=np.array([-2.9306374202572440192, -1.9816567566958429259,
            -1.1571937124467801947, -0.3811869902073221168,
             0.3811869902073221168,  1.1571937124467801947,
             1.9816567566958429259,  2.9306374202572440192])

w=np.array([1.9960407221136761921e-4, 1.7077983007413475456e-2,
            2.0780232581489187954e-1, 6.6114701255824129103e-1,
            6.6114701255824129103e-1, 2.0780232581489187954e-1,
            1.7077983007413475456e-2, 1.9960407221136761921e-4])

gauss=(w*np.cos(x)).sum()
exact=np.sqrt(np.pi)* np.exp(-1./4.)

print(" Exact={0:18.12e}\n Gauss={1:18.12e}\n Error={2:18.12e}".format(exact,  gauss, abs(exact-gauss)))

 Exact=1.380388447043e+00
 Gauss=1.380388447031e+00
 Error=1.184230491447e-11

The agreement is quite good.