Linear Algebra and Simultaneous Equations

2. Linear Algebra and Simultaneous Equations#

This chapter was last revised on August 14, 2024

References:

The main topic of this chapter is solving simultaneous linear equations, with a few sections on the related topics of computing eigenvalues and eigenvectors, and on solving system of simultaneous systems of non-linear equations by merging ideas from the previous chapter on Root-finding.

The problem of solving a system of \(n\) simultaneous linear equations in \(n\) unknowns, with matrix-vector form \(A x = b\), is quite thoroughly understood as far as having a good general-purpose methods usable with any \(n \times n\) matrix \(A\): essentially, Gaussian elimination (or row-reduction) as seen in most linear algebra courses, combined with some modifications to stay well away from division by zero: partial pivoting. Also, good robust software for this general case is readily available, for example in the Python packages NumPy and SciPy.

Nevertheless, one has to be careful about possibly dividing by zero oand the effects of rounding error, and this basic algorithm can be very slow when \(n\) is large – as it often is when dealing with differential equations (even more so with partial differential equations). We will see that it requires about \(n^3/3\) arithmetic operations.

Thus I will summarise the basic method of row reduction or Gaussian elimination, and then build on it with methods for doing things more robustly, and then with methods for doing it faster in some important special cases:

  1. When one has to solve many systems \(A x^{(m)} = b^{(m)}\) with the same matrix \(A\) but different right-hand side vectors \(b^{(m)}.\)

  2. When \(A\) is banded: most elements are zero, and all the non-zero elements \(a_{i,j}\) are near the main diagonal: \(|i - j|\) is far less than \(n\). (Aside on notation: “far less than” is sometimes denoted \(\ll\), as in \(|i-j| \ll n\).)

  3. When \(A\) is strictly diagonally dominant: each diagonal element \(a_{i,i}\) is larger in magnitude that the sum of the magnitudes of all other elements in the same row.

Other cases not (yet) discussed in this text are

  1. When \(A\) is positive definite: symmetric (\(a_{i,j} = a_{j,i}\)) and with all eigenvalues positive. This last condition would seem hard to verify, since computing all the eigenvalues of \(A\) is harder that solving \(Ax = b\), but there are important situations where this property is automatically guaranteed, such as with Galerkin and finite-element methods for solving boundary value problems for differential equations.

  2. When \(A\) is sparse: most elements are zero, but not necessarily with all the non-zero elements near the main diagonal.

Remark 2.1 (Module numpy.linalg with standard nickname la)

Python package Numpy provides a lot of useful tools for numerical linear algebra through a module numpy.linalg.

Just as package numpy is used so often that there is a conventional nickname np, so numpy.linalg is usually nicknamed la.

Thus we will often start a section or Jupyter notebook which uses linear algebra with the following import commands.

import numpy as np
import numpy.linalg as la

# Import some items individually, so they can be used by "first name only".
from numpy import array, inf, zeros_like, empty