# 2. Linear Algebra and Simultaneous Equations#

This chapter was last revised on August 14, 2024

**References:**

Chapter 2

*Systems of Equations*of [Sauer, 2022].Chapter 6

*Direct Methods for Solving Linear Systems*of [Burden*et al.*, 2016] and some early sections of Chapter 7.Chapter 8

*Linear Systems*of [Kincaid and Chenney, 1990], and parts of Chapter 8.

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:

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)}.\)

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\).)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

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.When \(A\) is

*sparse*: most elements are zero, but not necessarily with all the non-zero elements near the main diagonal.

(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
```