9. Numerical Routines: SciPy and NumPy¶

9.2.1. Single integrals¶

The function quad is the workhorse of SciPy’s integration functions. Numerical integration is sometimes called quadrature, hence the name. It is normally the default choice for performing single integrals of a function \(f(x)\) over a given fixed range from \(a\) to \(b\)

The general form of quad is scipy.integrate.quad(f, a, b), where f is the name of the function to be integrated and a and b are the lower and upper limits, respectively. The routine uses adaptive quadrature methods to numerically evaluate integrals, meaning it successively refines the subintervals (makes them smaller) until a desired level of numerical precision is achieved. For the quad routine, this is about \(10^{-8}\), although it usually does even better.

As an example, let’s integrate a Gaussian function over the range from 0 to 1

We first need to define the function \(f(x)=e^{-x^2}\), which we do using a lambda expression, and then we call the function quad to perform the integration.

In [1]: import scipy.integrate

In [2]: f = lambda x : exp(-x**2)

In [3]: scipy.integrate.quad(f, 0, 1)
Out[3]: (0.7468241328124271, 8.291413475940725e-15)

The function call scipy.integrate.quad(f, 0, 1) returns two numbers. The first is 0.7468..., which is the value of the integral, and the second is 8.29...e-15, which is an estimate of the absolute error in the value of the integral, which we see is quite small compared to 0.7468.

Because quad requires a function name as its first argument, we can’t simply use the expression exp(-x**2). On the other hand, we could use the usual def statement to create a normal function, and then use the name of that function in quad. However, it’s simpler here to use a lambda expression. In fact, we can just put the lambda expression directly into the first argument, as illustrated here

In [4]: scipy.integrate.quad(lambda x : exp(-x**2), 0, 1)
Out[4]: (0.7468241328124271, 8.291413475940725e-15)

That works too! Thus we see a lambda expression used as an anonymous function, a function with no name, as promised in the section Anonymous functions (lambda).

In [5]: scipy.integrate.quad(lambda x : exp(-x**2), 0, inf)
Out[5]: (0.8862269254527579, 7.101318390472462e-09)

In [6]: scipy.integrate.quad(lambda x : exp(-x**2), -inf, 1)
Out[6]: (1.6330510582651852, 3.669607414547701e-11)

The quad function handles infinite limits just fine. The absolute errors are somewhat larger but still well within acceptable bounds for practical work.

The quad function can integrate standard predefined NumPy functions of a single variable, like exp, sin, and cos.

In [7]: scipy.integrate.quad(exp, 0, 1)
Out[7]: (1.7182818284590453, 1.9076760487502457e-14)

In [8]: scipy.integrate.quad(sin, -0.5, 0.5)
Out[8]: (0.0, 2.707864644566304e-15)

In [9]: scipy.integrate.quad(cos, -0.5, 0.5)
Out[9]: (0.9588510772084061, 1.0645385431034061e-14)

Let’s integrate the first order Bessel function of the first kind, usually denoted \(J_1(x)\), over the interval from 0 to 5. Here is how we do it, using scipy.special.jn(v,x) where v is the (real) order of the Bessel function:

In [10]: import scipy.special

In [11]: scipy.integrate.quad(lambda x: scipy.special.jn(1,x),0,5)
Out[11]: (1.177596771314338, 1.8083362065765924e-14)

Because the SciPy function scipy.special.jn(v, x) is a function of two variables, v and x, we cannot use the function name scipy.special.jn in quad. So we use a lambda expression, which is a function of only one variable, x, because we have set the v argument equal to 1.

In [12]: p = np.poly1d([2, 5, 1])

In [13]: p
Out[13]: poly1d([2, 5, 1])

In [14]: p(1), p(2), p(3.5)
Out[14]: (8, 19, 43.0)

In [15]: P = polyint(p)

In [16]: P
Out[16]: poly1d([ 0.66666667,  2.5     ,  1.      ,  0.      ])

In [17]: q=P(5)-P(1)

In [18]: q
Out[18]: 146.66666666666666

9.2.2. Double integrals¶

The scipy.integrate function dblquad can be used to numerically evaluate double integrals of the form

The general form of dblquad is

scipy.integrate.dblquad(func, a, b, gfun, hfun)

where func if the name of the function to be integrated, a and b are the lower and upper limits of the x variable, respectively, and gfun and hfun are the names of the functions that define the lower and upper limits of the y variable.

As an example, let’s perform the double integral

We define the functions f, g, and h, using lambda expressions. Note that even if g, and h are constants, as they may be in many cases, they must be defined as functions, as we have done here for the lower limit.

In [19]: f = lambda x, y : 16*x*y

In [20]: g = lambda x : 0

In [21]: h = lambda y : sqrt(1-4*y**2)

In [22]: scipy.integrate.dblquad(f, 0, 0.5, g, h)
Out[22]: (0.5, 5.551115123125783e-15)

Once again, there are two outputs: the first is the value of the integral and the second is its absolute uncertainty.

Of course, the lower limit can also be a function of \(y\), as we demonstrate here by performing the integral

The code for this is given by

In [23]: g = lambda y : 1-2*y

In [24]: scipy.integrate.dblquad(f, 0, 0.5, g, h)
Out[24]: (0.33333333333333326, 3.700743415417188e-15)

9.4.1. Continuous and discrete Fourier transforms¶

The Fourier transform of a function \(g(t)\) is given by

where \(f\) is the Fourier transform variable; if \(t\) is time, then \(f\) is frequency. The inverse transform is given by

Here we define the Fourier transform in terms of the frequency \(f\) rather than the angular frequency \(\omega = 2\pi f\).

The conventional Fourier transform is defined for continuous functions, or at least for functions that are dense and thus have an infinite number of data points. When doing numerical analysis, however, you work with discrete data sets, that is, data sets defined for a finite number of points. The discrete Fourier transform (DFT) is defined for a function \(g_n\) consisting of a set of \(N\) discrete data points. Those \(N\) data points must be defined at equally-spaced times \(t_n=n\Delta t\) where \(\Delta t\) is the time between successive data points and \(n\) runs from 0 to \(N-1\). The discrete Fourier transform (DFT) of \(g_n\) is defined as

where \(l\) runs from 0 to \(N-1\). The inverse discrete Fourier transform (iDFT) is defined as

The DFT is usually implemented on computers using the well-known Fast Fourier Transform (FFT) algorithm, generally credited to Cooley and Tukey who developed it at AT&T Bell Laboratories during the 1960s. But their algorithm is essentially one of many independent rediscoveries of the basic algorithm dating back to Gauss who described it as early as 1805.

9.4.2. The SciPy FFT library¶

The SciPy library scipy.fftpack has routines that implement a souped-up version of the FFT algorithm along with many ancillary routines that support working with DFTs. The basic FFT routine in scipy.fftpack is appropriately named fft. The program below illustrates its use, along with the plots that follow.

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt

width = 2.0
freq = 0.5

t = np.linspace(-10, 10, 101)   # linearly space time array
g = np.exp(-np.abs(t)/width) * np.sin(2.0*np.pi*freq*t)

dt = t[1]-t[0]       # increment between times in time array

G = fftpack.fft(g)   # FFT of g
f = fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]
f = fftpack.fftshift(f)     # shift frequencies from min to max
G = fftpack.fftshift(G)     # shift G order to coorespond to f

fig = plt.figure(1, figsize=(8,6), frameon=False)
ax1 = fig.add_subplot(211)
ax1.plot(t, g)
ax1.set_xlabel('t')
ax1.set_ylabel('g(t)')

ax2 = fig.add_subplot(212)
ax2.plot(f, np.real(G), color='dodgerblue', label='real part')
ax2.plot(f, np.imag(G), color='coral', label='imaginary part')
ax2.legend()
ax2.set_xlabel('f')
ax2.set_ylabel('G(f)')

plt.show()

Function \(g(t)\) and its DFT \(G(f)\).

The DFT has real and imaginary parts, both of which are plotted in the figure.

The fft function returns the \(N\) Fourier components of \(G_n\) starting with the zero-frequency component \(G_0\) and progressing to the maximum positive frequency component \(G_{(N/2)-1}\) (or \(G_{(N-1)/2}\) if \(N\) is odd). From there, fft returns the maximum negative component \(G_{N/2}\) (or \(G_{(N-1)/2}\) if \(N\) is odd) and continues upward in frequency until it reaches the minimum negative frequency component \(G_{N-1}\). This is the standard way that DFTs are ordered by most numerical DFT packages. The scipy.fftpack function fftfreq creates the array of frequencies in this non-intuitive order such that f[n] in the above routine is the correct frequency for the Fourier component G[n]. The arguments of fftfreq are the size of the the orignal array g and the keyword argument d that is the spacing between the (equally spaced) elements of the time array (d=1 if left unspecified). The package scipy.fftpack provides the convenience function fftshift that reorders the frequency array so that the zero-frequency occurs at the middle of the array, that is, so the frequencies proceed monotonically from smallest (most negative) to largest (most positive). Applying fftshift to both f and G puts the frequencies f in ascending order and shifts G so that the frequency of G[n] is given by the shifted f[n].

The scipy.fftpack module also contains routines for performing 2-dimensional and \(n\)-dimensional DFTs, named fft2 and fftn, respectively, using the FFT algorithm.

As for most FFT routines, the scipy.fftpack FFT routines are most efficient if \(N\) is a power of 2. Nevertheless, the FFT routines are able to handle data sets where \(N\) is not a power of 2.

scipy.fftpack also supplies an inverse DFT function ifft. It is written to act on the unshifted FFT so take care! Note also that ifft returns a complex array. Because of machine roundoff error, the imaginary part of the function returned by ifft will, in general, be very near zero but not exactly zero even when the original function is a purely real function.

9.5.1. Basic computations in linear algebra¶

SciPy has a number of routines for performing basic operations with matrices. The determinant of a matrix is computed using the scipy.linalg.det function:

In [1]: import scipy.linalg
In [2]: a = array([[-2, 3], [4, 5]])
In [3]: a
Out[4]: array([[-2,  3],
               [ 4,  5]])

In [5]: scipy.linalg.det(a)
Out[5]: -22.0

The inverse of a matrix is computed using the scipy.linalg.inv function, while the product of two matrices is calculated using the NumPy dot function:

In [6]: b = scipy.linalg.inv(a)

In [6]: b
Out[6]: array([[-0.22727273,  0.13636364],
               [ 0.18181818,  0.09090909]])

In [7]: dot(a,b)
Out[7]: array([[ 1.,  0.],
               [ 0.,  1.]])

9.5.2. Solving systems of linear equations¶

Solving systems of equations is nearly as simple as constructing a coefficient matrix and a column vector. Suppose you have the following system of linear equations to solve:

The first task is to recast this set of equations as a matrix equation of the form \(\mathsf{A}\, \mathbf{x} = \mathbf{b}\). In this case, we have:

Next we construct the array \(\mathsf{A}\) and vector \(\mathbf{b}\) as NumPy arrays:

In [8]: A = array([[2, 4, 6], [1, -3, -9], [8, 5, -7]])
In [9]: b = array([4, -11, 2])

Finally we use the SciPy function scipy.linalg.solve to find \(x_1\), \(x_2\), and \(x_3\).

In [10]: scipy.linalg.solve(A,b)
Out[10]: array([ -8.91304348,  10.2173913 ,  -3.17391304])

which gives the results: \(x_1=-8.91304348\), \(x_2= 10.2173913\), and \(x_3= -3.17391304\). Of course, you can get the same answer by noting that \(\mathbf{x}=\mathsf{A}^{-1}\mathbf{b}\). Following this approach, we can use the scipy.linalg.inv introduced in the previous section:

Ainv = scipy.linalg.inv(A)

In [10]: dot(Ainv, b)
Out[10]: array([ -8.91304348,  10.2173913 ,  -3.17391304])

which is the same answer we obtained using scipy.linalg.solve. Using scipy.linalg.solve is numerically more stable and a faster than using \(\mathbf{x}=\mathsf{A}^{-1}\mathbf{b}\), so it is the preferred method for solving systems of equations.

You might wonder what happens if the system of equations are not all linearly independent. For example if the matrix \(\mathsf{A}\) is given by

where the third row is a multiple of the first row. Let’s try it out and see what happens. First we change the bottom row of the matrix \(\mathsf{A}\) and then try to solve the system as we did before.

In [11]: A[2] = array([1, 2, 3])

In [12]: A
Out[12]: array([[ 2,  4,  6],
                [ 1, -3, -9],
                [ 1,  2,  3]])

In [13]: scipy.linalg.solve(A,b)
LinAlgError: Singular matrix

In [14]: Ainv = scipy.linalg.inv(A)
LinAlgError: Singular matrix

Whether we use scipy.linalg.solve or scipy.linalg.inv, SciPy raises an error because the matrix is singular.

9.5.3. Eigenvalue problems¶

One of the most common problems in science and engineering is the eigenvalue problem, which in matrix form is written as

where \(\mathsf{A}\) is a square matrix, \(\mathbf{x}\) is a column vector, and \(\lambda\) is a scalar (number). Given the matrix \(\mathsf{A}\), the problem is to find the set of eigenvectors \(\mathbf{x}\) and their corresponding eigenvalues \(\lambda\) that solve this equation.

We can solve eigenvalue equations like this using scipy.linalg.eig. the outputs of this function is an array whose entries are the eigenvalues and a matrix whose rows are the eigenvectors. Let’s return to the matrix we were using previously and find its eigenvalues and eigenvectors.

A = array([[2, 4, 6],[1, -3, -9],[8, 5, -7]])

In [15]: A
Out[15]: array([[ 2,  4,  6],
                [ 1, -3, -9],
                [ 8,  5, -7]])

In [16]: lam, evec = scipy.linalg.eig(A)

In [17]: lam
Out[17]: array([ 2.40995356+0.j, -8.03416016+0.j,
                -2.37579340+0.j])

In [18]: evec
Out[18]: array([[-0.77167559, -0.52633654,  0.57513303],
                [ 0.50360249,  0.76565448, -0.80920669],
                [-0.38846018,  0.36978786,  0.12002724]])

The first eigenvalue and its corresponding eigenvector are given by

In [19]: lam[0]
Out[19]: (2.4099535647625494+0j)

In [20]: evec[:,0]
Out[20]: array([-0.77167559,  0.50360249, -0.38846018])

We can check that they satisfy the \(\mathsf{A}\mathbf{x} = \lambda \mathbf{x}:\)

In [21]: dot(A,evec[:,0])
Out[21]: array([-1.85970234,  1.21365861, -0.93617101])

In [22]: lam[0]*evec[:,0]
Out[22]: array([-1.85970234+0.j,  1.21365861+0.j,
                -0.93617101+0.j])

Thus we see by direct substitution that the left and right sides of \(\mathsf{A}\mathbf{x} = \lambda \mathbf{x}:\) are equal. In general, the eigenvalues can be complex, so their values are reported as complex numbers.

In [22]: A = array([[2, 4, 6], [1, -3, -9], [8, 5, -7]])
Out[22]: B = array([[5, 9, 1], [-3, 1, 6], [4, 2, 8]])

In [23]: lam, evec = scipy.linalg.eig(A,B)

In [24]: lam
Out[24]: array([-1.36087907+0.j,  0.83252442+0.j,
                -0.10099858+0.j])

In [25]: evec
Out[25]: array([[-0.0419907 , -1.        ,  0.93037493],
                [-0.43028153,  0.17751302, -1.        ],
                [ 1.        , -0.29852465,  0.4226201 ]])

9.6.1. Single equations of a single variable¶

Solving a single nonlinear equation is enormously simpler than solving a system of nonlinear equations, so that is where we start. A word of caution: solving non-linear equations can be a tricky business so it is important that you have a good sense of the behavior of the function you are trying to solve. The best way to do this is to plot the function over the domain of interest before trying to find the solutions. This will greatly assist you in finding the solutions you seek and avoiding spurious solutions.

We begin with a concrete example. Suppose we want to find the solutions to the equation

Plots of \(\tan x\) and \(\sqrt{(8/x)^2-1}\) vs \(x\) are shown in the top plot in the figure Crossing functions, albeit with \(x\) replaced by \(\theta\). The solutions to this equation are those \(x\) values where the two curves \(\tan x\) and \(\sqrt{(8/x)^2-1}\) cross each other. The first step towards obtaining a numerical solution is to rewrite the equation to be solved in the form \(f(x)=0\). Doing so, the above equation becomes

Obviously the two equations above have the same solutions for \(x\). Parenthetically we mention that the problem of finding the solutions to equations of the form \(f(x)=0\) is often referred to as finding the roots of \(f(x)\).

Next, we plot \(f(x)\) over the domain of interest, in this case from \(x=0\) to 8. We are only interested in positive solutions and for \(x>8\), the equation has no real solutions as the argument of the square root becomes negative. The solutions, the points where \(f(x)=0\) are indicated by green circles; there are three of them. Another notable feature of the function is that it diverges to \(\pm\infty\) at \(x = \{0, \pi/2, 3\pi/2, 5\pi/2\}\).

Roots of a nonlinear function

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt

def tdl(x):
    y = 8./x
    return np.tan(x) - np.sqrt(y*y-1.0)

# Find true roots

rx1 = scipy.optimize.brentq(tdl, 0.5, 0.49*np.pi)
rx2 = scipy.optimize.brentq(tdl, 0.51*np.pi, 1.49*np.pi)
rx3 = scipy.optimize.brentq(tdl, 1.51*np.pi, 2.49*np.pi)
rx = np.array([rx1, rx2, rx3])
ry = np.zeros(3)
# print using a list comprehension
print('\nTrue roots:')
print('\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rx))

# Find false roots

rx1f = scipy.optimize.brentq(tdl, 0.49*np.pi, 0.51*np.pi)
rx2f = scipy.optimize.brentq(tdl, 1.49*np.pi, 1.51*np.pi)
rx3f = scipy.optimize.brentq(tdl, 2.49*np.pi, 2.51*np.pi)
rxf = np.array([rx1f, rx2f, rx3f])
# print using a list comprehension
print('\nFalse roots:')
print('\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rxf))

# Plot function and various roots

x = np.linspace(0.7, 8, 128)
y = tdl(x)
# Create masked array for plotting
ymask = np.ma.masked_where(np.abs(y)>20., y)

plt.figure(figsize=(6, 4))
plt.plot(x, ymask)
plt.axhline(color='black')
plt.axvline(x=np.pi/2., color="gray", linestyle='--', zorder=-1)
plt.axvline(x=3.*np.pi/2., color="gray", linestyle='--', zorder=-1)
plt.axvline(x=5.*np.pi/2., color="gray", linestyle='--', zorder=-1)
plt.xlabel(r'$x$')
plt.ylabel(r'$\tan x - \sqrt{(8/x)^2-1}$')
plt.ylim(-8, 8)

plt.plot(rx, ry, 'og', ms=5, label='true roots')

plt.plot(rxf, ry, 'xr', ms=5, label='false roots')
plt.legend(numpoints=1, fontsize='small', loc = 'upper right',
           bbox_to_anchor = (0.92, 0.97))
plt.tight_layout()
plt.show()

In [1]: run rootbrentq.py

True roots:
f(1.39547) = -6.39e-14
f(4.16483) = -7.95e-14
f(6.83067) = -1.11e-15

False roots:
f(1.57080) = -1.61e+12
f(4.71239) = -1.56e+12
f(7.85398) = 1.16e+12

9.6.2. Solving systems of nonlinear equations¶

Solving systems of nonlinear equations is not for the faint of heart. It is a difficult problem that lacks any general purpose solutions. Nevertheless, SciPy provides quite an assortment of numerical solvers for nonlinear systems of equations. However, because of the complexity and subtleties of this class of problems, we do not discuss their use here.