MATH 417 Lecture 18

« previous | Thursday, March 27, 2014 | next »

Quiz on Tuesday over Chapter 5 (diff eqs)

Matrices

Component-wise addition, subtraction, and scalar multiplication

Gaussian-Jordan Elimination

(Sections 6.1 and 6.2 Row echelon form:

Find solution to $A\,{\vec {x}}={\vec {b}}$

${\begin{array}{c|c}A&b\end{array}}\to {\begin{array}{c|c}I&x\end{array}}$

1. If $a_{1,1}=0$, find another row $j$ such that $a_{j1}\neq 0$ and switch
2. Divide row $1$ by $a_{11}$
3. Add row 1 scaled by $-a_{i,1}$ to remaining rows $i$
4. Recur on submatrices to get upper triangular matrirx
5. Jordan steps perform similar operations to upper half of matrix (starting from bottom) to get identity matrix.

Gaussian and Jordan's steps each take $O\left({\frac {2}{3}}\,n^{3}\right)$ operations. Hence solving system of equations is $O\left({\frac {4}{3}}\,n^{3}\right)$

In Matlab, x = A \ b is a fast matrix equation solver. Don't use x = A^(-1) * b

Instability

Example

Let $\epsilon <{\frac {1}{2}}$

${\begin{bmatrix}\epsilon &1\\1&1\end{bmatrix}}\,{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}}$

$\left[{\begin{array}{cc|c}\epsilon &1&1\\1&1&2\end{array}}\right]$

Performing gaussian elimination gives solutions

$x_{1}={\frac {1}{\epsilon }}\,\left(1-x_{2}\right)$ $x_{2}={\frac {2-{\frac {1}{\epsilon }}}{1-{\frac {1}{\epsilon }}}}$ for very small $\epsilon$, we get $x_{2}\approx {\frac {-{\frac {1}{\epsilon }}}{-{\frac {1}{\epsilon }}}}=1$, which is actually a pretty good estimation, but substituting it back into $x_{1}gives{\frac {1}{\epsilon }}\cdot 0=0$, whereas the exact solution for $\epsilon =0$ is $x_{1}=x_{2}=1$.

In reality, we are not solving for $A\,{\vec {x}}={\vec {b}}$, we are actually solving for $\left(A+\Delta A\right)\,{\hat {x}}=\left({\vec {b}}+\Delta {\vec {b}}\right)$. We really care how close ${\hat {x}}$ is to ${\vec {x}}$ given $\left|\Delta A\right|$ and $\left|\Delta b\right|$ are small values.

Matlab Exercise

${\vec {x}}={\begin{bmatrix}1\\1\\\vdots \\1\end{bmatrix}}$

Generate random matrix $A$ and compute $A\,{\vec {x}}={\vec {b}}$

Solve for $A\,{\hat {x}}=b$.

This usually works only 1 out of 10 times, so Gaussian elimination is unstable

Gaussian Elimination with Partial Pivoting

Added stability

At step $i$, we are operating on row $i$ with current diagonal entry $a_{ii}$.

If $a_{ii}=0$, find $j$ such that $\left|a_{ji}\right|=\max _{i\leq k\leq n}\left|a_{ki}\right|$ and swap rows $i$ and $j$.

Now whenever we multiply remaining rows $k$ by $-\left({\frac {a_{ki}}{a_{ii}}}\right)$, we have $\left|{\frac {a_{ki}}{a_{ii}}}\right|\leq 1$

Example

Applying this to our example above gives

${\begin{bmatrix}\epsilon &1\\1&1\end{bmatrix}}\,{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}}$

The solutions match with exact solutions: $x_{2}={\frac {1-2\epsilon }{1-\epsilon }}\approx {\frac {1}{1}}=1$ and $x_{1}=2-x_{2}\approx 1$

Gaussian Elimination with Full Pivoting

When looking at $a_{11}$,

1. find $a_{ij}$ such that $\left|a_{ij}\right|=\max _{1\leq k,s\leq n}\left|a_{ks}\right|$ (largest entry in entire matrix)
2. Interchange rows $1$ and $i$ and columns $1$ and $j$
3. Change $x_{1}$ to $x_{j}$ (bookkeeping headache)

This is better, but still not perfect: there are still matrices that might cause it to fail.

(this method is what Matlab uses)

Eigenvalues and Eigenvectors

$\lambda$ is an eigenvalue of $A$ if there exists ${\vec {x}}\neq {\vec {0}}$ such that $A\,{\vec {x}}=\lambda \,{\vec {x}}$.

This implies $0=A\,{\vec {x}}-\lambda \,{\vec {x}}=\left(A-\lambda \,I\right)\,{\vec {x}}$, which further implies $\left|A-\lambda \,I\right|=0$.

Computing this determinant gives a polynomial of $\lambda$ of degree $n$ (for $n$ \times $n$ matrix) that has roots $\lambda$ (i.e. $p(\lambda )=0$).

Example

${\begin{vmatrix}\epsilon -\lambda &1\\1&1-\lambda \end{vmatrix}}=0$

Gives equation $\lambda ^{2}-\left(1+\epsilon \right)\,\lambda +\epsilon -1=0$

and solutions $\lambda ={\frac {1+\epsilon \pm {\sqrt {\left(1+\epsilon \right)^{2}-4(\epsilon -1)}}}{2}}$

As $\epsilon \to 0$, we get $\left|{\frac {\lambda _{2}}{\lambda _{1}}}\right|\approx 3$. Therefore Gaussian elimination will be stable.