Solve the following problem by hand using Gaussian Elimination (row reduction): \begin{equation*} \begin{pmatrix}1 \amp 2 \amp 3 \\ 4 \amp 5 \amp 6 \\ 7 \amp 8 \amp 0 \end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix}1 \\ 0 \\ 2 \end{pmatrix} \end{equation*}

Open a new MATLAB script and enter the left-hand matrix of the previous problem into matrix \(A\). Observe what happens when we do the following in sequence:

• left multiply \(A\) by the matrix \begin{equation*} L_1 = \begin{pmatrix}1 \amp 0 \amp 0 \\ -4 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1 \end{pmatrix} \end{equation*} This gives the matrix \(L_1 A\).

• left multiply \(L_1 A\) by the matrix \begin{equation*} L_2 = \begin{pmatrix}1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ -7 \amp 0 \amp 1 \end{pmatrix} \end{equation*} This gives the matrix \(L_2 L_1 A\).

• left multiply \(L_2 L_1 A\) by the matrix \begin{equation*} L_3 = \begin{pmatrix}1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 0 \amp -2 \amp 1 \end{pmatrix} \end{equation*}

Make a conjecture: If you wanted to multiply row \(j\) of n \(n\times n\) matrix by \(c\) and add it to row \(k\), that is the same as multiplying by what lower triangular matrix?

After the process from the previous problem you should notice that you now have an upper triangular matrix. Hence, in general, we have done this: \begin{equation*} L_3 L_2 L_1 A = U, \end{equation*} so if you solve for \(A\) we see that \(A\) can be written as \begin{equation*} A = L_1^{-1} L_2^{-1} L_3^{-1} U. \end{equation*}

Therefore, we could rewrite the problem \(A \bx = \bb\) as the problem \(LU\bx = \bb\) where \(L\) is which matrix?

If you have a lower triangular matrix of the forms discussed in the previous problem, how do you find it's inverse? (Experiment with MATLAB and work with your neighbors)

Now for the punch line: If we can write \(A \bx = \bb\) then if we can write it as \(LU \bx = \bb\) we can

1. Solve \(L \by = \bb\) for \(\by\) using forward substitution. Then,

2. solve \(U \bx = \by\) for \(\bx\) using backward substitution.

For our running example, write down \(L\) and \(U\) and solve the system.

Try the process again on the \(3\times 3\) system of equations \begin{equation*} \begin{pmatrix}3 \amp 6 \amp 8\\ 2 \amp 7 \amp -1 \\ 5 \amp 2 \amp 2 \end{pmatrix} \begin{pmatrix}x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix}-13 \\ 4 \\ 1 \end{pmatrix} \end{equation*}

That is: Find matrices \(L\) and \(U\) such that \(A \bx = \bb\) can be written as \(LU\bx = \bb\). The do two triangular solves to determine \(\bx\).

Notice that this time there isn't a “\(1\)” in the top left corner to begin with. Be careful.

Write three MATLAB functions:

1. Write a function called “Lsolve.m” that takes inputs a unit lower triangular matrix \(L\) and a vector \(\bb\) and outputs a vector \(\by\) that solves the equation \(L \by = \bb\).

2. Write a function called “Usolve.m” that takes inputs an upper triangular matrix \(U\) and a vector \(\by\) and outputs a vector \(\bx\) that solves the equation \(U \bx = \by\).

3. Write a function called “MyLU.m” that takes a square matrix \(A\) and outputs the matrices \(L\) and \(U\) such that \(A = LU\).

Test each of your functions individually on the problems that we've done thus far. (You can use the codes on pages 140 and 141 to get you going, but I'm asking for slightly different outputs on some of these so be careful). Then test your problem on a system that you don't know the solution to. Discuss where your code will fail. Use the following partial code to test your functions.

A = [...; ...; ... ];

b= [...; ...; ... ];

[L,U] = MyLU(A);

y = Lsolve(L,b);

x = Usolve(U,y)

MATLABExactAnswer = A\b % uses the MATLAB system solver

MyError = norm(x-MATLABExactAnswer)

For this problem we are going to run a numerical experiment to see how the process of solving the the equation \(A \bx = \bb\) using the \(LU\) factorization. Here is what you need to do:

Create a loop that does the following

• Build a random matrix of size \(n\times n\). You can do this with the code: A=rand(n,n)
• Build a random vector in \(\mathbb{R}^n\). You can do this with the code: b = rand(n,1)
• Find MATLAB's exact answer to the problem \(A\bx = \bb\) using the backslash: Xexact = A\b
• Write code that uses your three \(LU\) functions (MyLU, Lsolve, Usolve) to fina a solution to the equation \(A\bx = \bb\).
• Find the error between your answer and the exact answer using the code: error = norm(x - Xexact)
• Make a plot that shows how the error behaves as the size of the problem changes. You should run this for matrices of larger and larger size but be warned that the loop will run for quite a long time if you go above \(300 \times 300\) matrices. Just be patient.