Introduction to Python for Science, Release 0.9.23

The polynomial functions shown have a special syntax that uses NumPy’s polyval function for generating polynomials. If p is a list or array of N numbers and x is an array, then

polyval(p, x) = p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]

For example, if p = [2.0, 5.0, 1.0], polyval generates the following quadratic polynomial:

polyval(p, x) = 2.0*x**2 + 5.0*x + 1.0

SciPy’s special.legendre(n) and special.laguerre(n) functions output the coefficients p needed in polyval to produce the nth-order Legendre and Laguerre polynomials, respectively. The special library of SciPy has functions that specify many other polynomial functions in this same way.

9.2 Linear algebra

Python’s mathematical libraries, NumPy and SciPy, have extensive tools for numerically solving problems in linear algebra. Here we focus on two problems that arise commonly in scientific and engineering settings: (1) solving a system of linear equations and (2) eigenvalue problems. In addition, we also show how to perform a number of other basic computations, such as finding the determinant of a matrix, matrix inversion, and LU decomposition. The SciPy package for linear algebra is called scipy.linalg.

9.2.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 [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:

Introduction to Python for Science, Release 0.9.23

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

[0., 1.]])

9.2.2Solving 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:

2x1 + 4x2 + 6x3 = 4 x1 3x2 9x3 = 11 8x1 + 5x2 7x3 = 1

The first task is to recast this set of equations as a matrix equation of the form A x = b.

Next we construct the array A and vector 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 x1, x2, and x3.

In [10]: scipy.linalg.solve(A,b)

Out[10]: array([ -8.91304348, 10.2173913 , -3.17391304])

which gives the results: x1 = 8:91304348, x2 = 10:2173913, and x3 = 3:17391304. Of course, you can get the same answer by noting that x = A 1b. 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])

Introduction to Python for Science, Release 0.9.23

which is the same answer we obtained using scipy.linalg.solve. Using scipy.linalg.solve is numerically more stable and a faster than using x = A 1b, 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 A is given by

01

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 A and then try to solve the system as we did before.

In [11]: A[2] = array([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.2.3 Eigenvalue problems

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

Ax = x

where A is a square matrix, x is a column vector, and is a scalar (number). Given the matrix A, the problem is to find the set of eigenvectors x and their corresponding eigenvalues 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

Introduction to Python for Science, Release 0.9.23

eigenvectors. Let’s return to the matrix we were using previously and find its eigenvalues and eigenvectors.

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 Ax = 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 Ax = x : are equal. In general, the eigenvalues can be complex, so their values are reported as complex numbers.