Linear Equations
Linear Equations
A system of linear equations (or linear system) is a collection of one or more linear equations involving the same variables. If there are $m$ equations with $n$ unknown variables to solve for
$$ \begin{align*} a_{1,1} x_1 + a_{1,2} x_2 + \ldots + a_{1,n} x_n & = b_1 \\ a_{2,1} x_1 + a_{2,2} x_2 + \cdots + a_{2,n} x_n & = b_2 \\ \vdots & \\ a_{m,1} x_1 + a_{m,2} x_2 + \cdots + a_{m,n} x_n & = b_m. \end{align*} $$The system of linear equations can be written in matrix form ${Ax=b}$, where
$$ \begin{align*} A = \left( \begin{array}{cccc} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{array} \right), \quad x = \left( \begin{array}{c} x_1 \\\ x_2 \\\ \vdots \\\ x_n \end{array} \right), \quad b = \left( \begin{array}{c} b_1 \\ b_2 \\ \vdots \\ b_m \end{array} \right), \end{align*} $$so that $A \in \mathbb{R}^{m \times n}$, $x \in \mathbb{R}^n$ and $b \in \mathbb{R}^m$.
A square matrix $A$ is symmetric if ${A=A^T}$, that is, $a_{i,j}=a_{j,i}$ for all indices $i$ and $j$.
A square matrix is said to be Hermitian if the matrix is equal to its conjugate transpose, i.e. $a_{i,j}=\overline{a_{j,i}}$ for all indices $i$ and $j$. A Hermitian matrix is written as $A=A^H$.
For a nonsingular matrix, the following all hold:
- Nonsingular matrix has full rank
- A square matrix is nonsingular if and only if the determinant of the matrix is non-zero.
- If a matrix is singular, both versions of Gaussian elimination will fail due to division by zero, yielding a floating exception error.
Direct Methods
Gaussian elimination is a method to solve systems of linear equations based on forward elimination (a series of row-wise operations) to convert the matrix, $A$, to upper triangular form (echelon form), and then back substitution to solve the system. The row operations are:
- row swapping
- row scaling, i.e. multiplying by a non-zero scalar
- row addition, i.e. adding a multiple of one row to another
Forward elimination is written as
1 For $k=1$ to $n-1$:
2 $\hspace{2cm}$ For $i = k + 1$ to $n$:
3 $\hspace{4cm}$ For $j = k$ to $n$:
4 $\hspace{6cm} a_{i,j} = a_{i,j} - \dfrac{a_{i,k}}{a_{k,k}} a_{k,j}$
5 $\hspace{4cm}$ End
6 $\hspace{2cm} b_{i} = b_{i} - \dfrac{a_{i,k}}{a_{k,k}} b_{k}$
7 End
Then back substitute, starting with the last unknown first:
8 Set $x_n = \dfrac{b_n}{a_{n,n}}$
9 For $i$ =$n-1$ to $1$:
10 $\hspace{2cm} y = b_i$
11 $\hspace{2cm}$ For $j = n$ to $i+1$:
12 $\hspace{4cm} y = y - a_{i,j} x_j$
13 $\hspace{2cm}$ End
14 $\hspace{2cm} x_i = \dfrac{y}{a_{i,i}}$
15 End
A pivot element is the element of a matrix, $A$, which is selected to do certain calculations first. Pivoting helps reduce errors due to rounding.
1 Find maximal absolute values vector $s$ with entries $s_i= \max_{j=1,\ldots,n}\left| a_{i,j}\right|$
2 For $k=1$ to $n-1$:
3 $\hspace{2cm}$ For $i = k$ to $n$:
4 $\hspace{4cm}$ Compute $\left| \dfrac{a_{i,k}}{s_i} \right|$
5 $\hspace{2cm}$ End
6 $\hspace{2cm}$ Find row with largest relative pivot element, denote this as row $j$
7 $\hspace{2cm}$ Swap rows $k$ and $j$
8 $\hspace{2cm}$ Swap entries $k$ and $j$ in vector $s$
9 $\hspace{2cm}$ Do forward elimination on row $k$
10 End
The matrix will be in row-echelon form, so that back substitution can now be performed.
Cholesky’s algorithm to compute the entries of the triangular matrix $L$ can be written as
1 For $i=1$ to $n$:
2 $\hspace{2cm}$ For $j = 1$ to $i-1$:
3 $\hspace{4cm}$ Set $y = a_{i,j}$
4 $\hspace{4cm}$ For $k = 1$ to $j-1$:
5 $\hspace{6cm}$ $y = y - l_{i,k} l_{j,k}$
6 $\hspace{4cm}$ End
7 $\hspace{4cm}$ $l_{i,j} = y / l_{j,j}$
8 $\hspace{2cm}$ End
9 $\hspace{2cm}$ Set $y = a_{i,i}$
10 $\hspace{2cm}$ For $k= 1$ to $i-1$:
11 $\hspace{4cm}$ $y = y - l_{i,k} l_{i,k}$
12 $\hspace{2cm}$ End
13 $\hspace{2cm}$ If $y \le 0$:
14 $\hspace{4cm}$ exit as there is no solution
15 $\hspace{2cm}$ Else:
16 $\hspace{4cm}$ $l_{i,i} = \sqrt{y}$
12 $\hspace{2cm}$ End
17 End
Indirect Methods
$$ Q x_{k +1}= \left( Q - A \right) x_{k} + b. $$This is equivalent to
$$ x_{k +1} = \left( I - Q^{-1} A \right) x_{k} + Q^{-1} b. $$where the $\lambda_i$ are the eigenvalues of the matrix.
This sequence may converge to $x = A^{-1} b$, depending on $A$.
$$ x_{k +1}= x_{k} + \omega \left( b - A x_{k} \right). $$where $D$ is the matrix comprised of the diagonal elements of $A$.
where $L$ is the lower triangular part of $A$, and $D$ is the diagonal matrix of $A$.
Fundamental Theorem of Numerical Analysis
A numerical method is said to be stable if and only if any initial error $e_0$ is damped during the iterations, i.e. ${\left\| e_k \right\| < \left\| e_0 \right\| }$.
Note that $\Vert x \Vert $ is a norm of a vector, such as ${\Vert x \Vert_2 = \sqrt{x_0^2 + x_1^2 + \ldots + x_n^2} }$.
If the iterative method for a linear system is stable then
$$ e_k = \left( I - Q^{-1} A \right)^k e_0, $$so then ${\Vert I - Q^{-1} A \Vert < 1}$ for ${\Vert e_k \Vert < \Vert e_0 \Vert}$, where ${\Vert I - Q^{-1} A \Vert}$ is a matrix norm.