Question 1. Let A = LU be a LU-factorisation of a NxN matrix A into a product of a lower-triangular matrix L and a upper-triangular matrix U. The solution x = U^-1L^-1b of Ax = b can be found efficiently by a "forward substitution" followed by a "backward substitution" using the following formulas derived in section §4.3.2 of the Lecture Notes

[L^-1b]_i = z_i = (b_i - _j=1∑^i-1l_ijz_j ) 1/l_ii, for i = 1, 2....., N,

[U^-1L^-1b] = x_i = (z_i - _j=i+1∑^N u_ijx_j ) 1/u_ii, for i = N, N-1, ....., 1,

Matlab code for forward substitution (1) is provided in Listing 1. Observe how elements of b that are no longer needed are replaced. Write a Matlab function b = backsolve(U,b) to perform backward substitution.

Listing 1: Matlab function implementing equation (1) for the solution of a lower-triangular system by forward substitution.

function b = for wardsolve (L, b)

[N , N ] = size(L);

b (1)= b (1)/ L (1 ,1);

for i = 2: N

b (i) = ( b ( i) - L (i ,1: i -1)* b (1: i -1))/ L (i , i );

end

Question 2. Let A = LU be a LU-factorisation of a given N × N matrix A into a product of a unit-lower-triangular matrix L with l_ii = 1 and a upper-triangular matrix U, to be determined. Show that the elements of U and L are given by

u_ij = (a_ij - _k=1∑^i-1l_iku_kj ), for i = 1, 2,..... j,     (3)

l_ij = (a_ij - _k=1∑^j-1l_iku_kj )1/ujj, for i = j + 1, j + 2,..... ,N,     (4)

for every j = 1, 2, ...., N .

[HINTS: Start from the explicit expression for matrix multiplication. Take into account the location of zeroes in L and U . May be useful to do Question 3 below before you do Question 2.]

Question 3. Illustrate the usage of equations (3) and (4) and find the LU-factorisation of the square matrix

[HINTS: Follow Example 4.29 from the Lecture Notes but proceed column-by-column rather than crosswise.]

Observe the following.

1. The u's and l's that occur in the right-hand side of (3) and (4) are already determined by the time they are needed. This makes the method possible and efficient.

2. The a's corresponding to u's and l's that are already determined are never needed again. This allows to overwrite the matrix A which is useful for pivoting, see below, and for reducing computer memory requirements.

3. Equation (3) in the case of i = j (its final application) is exactly the same as equation (4) except for the division in the latter equation; in both cases the upper limit of the sum is k = j - 1(= i - 1). This fact allows to incorporate partial pivoting (to select a large by modulus pivot (diagonal element ujj) for the division in (4)), which is essential for the method to work. Practically, compute l'_ij = (a_ij - j-1Σ_k=1 likukj), swap
the rows of A (which you have been overwriting with L and U ) corresponding to u_jj and m = max {|u_jj|,|l'_{j +1,j}|, |l'_{j +2j}|,....|l'_N,j| and only then divide by the pivot m to fully implement equation (4). Keep a permutation matrix P as you are now actually factorising PA = LU, instead.

Question 4. Write a Matlab function [P,L,U] = ColumnCroutLU(A) to perform the column- wise Crout method given by equations (3) and (4), overwriting the matrix A and performing partial pivoting as described in Question 3. Your function must take as a single argument the square matrix A and return a permutation matrix P, a unit-lower-triangular matrix L and upper-triangular matrix U of the same size. A Matlab function to perform one pivoting row swap is given in Listing2.

Use your ColumnCroutLU to factorise the matrices,

switching pivoting on and off. Verify that your factorisations satisfy PA = LU.

[HINTS: The code of Listing1can serve as an example. The inbuilt Matlab functions tril and triu can be used to extract L and U from the overwritten A.]

Listing 2: Matlab function to perform one pivoting row swap. Takes as arguments a matrix A to be pivoted and a pivot position (j,j) and returns a permutation matrix P and a modified A.

function[P , A ] = ppivot (A , j) [N , N ] =size( A );
P =eye( N );
[ amax , s] =max(abs( A ( j:N , j ))); s = s+ j -1;
A ([ j , s ] ,:) = A ([ s , j ] ,:);
P ([ j , s ] ,:) = P ([ s , j ] ,:);
end

Question 5. Study and run the Matlab function MysteriousCode provided in Listing3. What problem does it solve? Explain the related theory. Explain and comment the given code.

Listing 3: A Matlab function that performs a mysterious computation all the while using your own functions written for Questions 1 and 4.

function MysteriousCode
x = [0:1:10]';
y = cos( x );
A =vander( x );
[P ,L , U ]= ColumnCroutLU(A);
checkfact = L*U - P*A
z = forwardsolve (L, P * y );
a= backsolve (U, z );
xx = [0:0.1:10] ' ;
B = vander( xx );
BB = B (: ,end-10:end);
yy = BB * a; zz = cos( xx );
err =abs( zz - yy );
figure; plot(x ,y ,'-o');
hold on;plot( xx , yy ,'r');hold on;
figure;plot( xx , err );
end