Homework Set #10, Math 575A

Part I: This part is the analytical part.

1. Solve the following system twice. First, use the basic Gaussian elimination using the standard algorithm. Second, determine the factorization of the form $P A = L U$. Third, use the Gaussian elimination with scaled pivoting using the scaled pivoting algorithm.
   1. 
      
      $$\left[ \begin{array}{rrr} -1 & 1 & -4 \\ 2 & 2 & 0 \\ 3 ... ...\left[ \begin{array}{r} 0\ 1\ \frac{1}{2}\end{array}\right]$$
      
      

   2. 
      
      $$\left[ \begin{array}{rrrr} 6 & -2 & 2 & 4 \\ 12 & -8 & 4 &... ...left[ \begin{array}{r} 0\ -10\ -39\ -16 \end{array}\right]$$
      
      

2. Count how many short operations (additions and subtractons) and long operations (multiplications and divisions) are involved in solving a tridiagonal system following following the Crout algorithm. Count how many short operations and long operations are involved in the Gaussian elimination with scaled row pivoting. Give the detailed process of your counting. To check your counting it is suggested that you put counters in your code and compare your estimate to the one given in the code for many different values of $n$, the matrix size.

3. Look up ILU, the incomplete LU decomposition. Also look up the ILU with thresholding, in which the smaller elements of $L$ and $U$ are discarded.

   For the matrix
   

   $$\left[ \begin{array}{rrrr} 10 & 0.1 & 0.1 & 0.1 \\ 0.1 & 0... ... 0.1 &0 & 0.1 & 0 \\ 0.1 & 0 & 0 &0.1 \end{array} \right]$$
   
   
   1. Compute the condition number for the above matrix.
   2. Compute the LU decomposition
   3. Compute an ILU preconditioned matrix and its condition number.
   4. What can you conclude regarding the condition number of the original matrix and the ILU decomposition?
   5. What about the cost of the ILU?

Part II: This part should be done using MATLAB.

1. The following is a MATLAB implementation of the basic Gaussian elimination algorithm for solving the linear system of equations $A x = b$. The matlab script can be downloaded by clicking here. Modify this code to implement the algorithm of the Gaussian elimination with scaled row pivoting. Note that since you are only to use the original matrix for storage, you will destroy the original matrix entries in the process. Test your code with the the linear systems in problem #1 of Part I to verify the correctness of your code. You should output the permutation vector, the LU factorization of the matrix in addition to the solution. Hand in the code and output. Both the basic Gaussian elimination and and the Gaussian elimination with scaled row pivoting are important subjects for this course. Please be advised that you must show that you understand every detail of the algorithm.

   % Gauss_np.m
   % Gaussian Elimination without pivoting.
   
     clear;
     format short;
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Step 0: Assign the matrix A and the vector b
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     n = 4;
     A = [ 6,  -2, 2,   4;
          12,  -8, 6,  10;
           3, -13, 9,   3;
          -6,   4, 1, -18 ];
     b = [ 0.23; 0.32; 0.33; 0.31 ];
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Step 1: Basic Gaussian Elimination. 
   % The matrix A is overwritten to store its LU factorization.
   %   L is unit lower triangular:
   %      L(i,j) = 0 for j>i; L(i,i) = 1;  L(i,j) = A(i,j) for i>j.
   %   U is upper triangular:
   %      U(i,j) = A(i,j) for j>=i; U(i,j) = 0 for i>j.
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     for k = 1:(n-1)
       for i = (k+1):n
           A(i,k) = A(i,k) / A(k,k);
         for j = (k+1):n
           A(i,j) = A(i,j) - A(i,k) * A(k,j);
         end
       end
     end
   
   % output the LU factorization of the original matrix A.
     A
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Step 2: Forward substitution to solve  L * y = b where 
   % L(i,j) = 0 for j>i; L(i,i) = 1;  L(i,j) = A(i,j) for i>j.
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     y = zeros(n,1);                % Initialize y to be a column vector
                                    % of length n with zeros entries.
     y(1) = b(1);
     for i = 2:n
         y(i) = b(i);
       for j = 1:(i-1)
         y(i) = y(i) - A(i,j)*y(j);
       end
     end
   
   % output y (the intermediate solution):
     y
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Step 3: Back substitution to solve  U * x = y
   % U(i,j) = A(i,j) for j>=i; U(i,j) = 0 for i>j.
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     x = zeros(n,1);                % Initialize x to be a column vector
                                    % of length n with zeros entries.
     x(n) = y(n) / A(n,n);
     for i = (n-1):-1:1
         x(i) = y(i);
       for j = (i+1):n
         x(i) = x(i) - A(i,j)*x(j);
       end
         x(i) = x(i) / A(i,i);
     end
   
   % output x (the final solution):
     x