MATH2071: LAB 6: Solving linear systemsIntroduction Exercise 1Test matrices Exercise 2The linear system “problem” Exercise 3Gauß factorization Exercise 4Permutation matrices Exercise 5PLU factorization Exercise 6PLU solution Exercise 7A system of ODEs Exercise 8Extra Credit: Pivoting Exercise 9Exercise 10Exercise 11Extra Credit1 IntroductionIn many numerical applications, notably the ordinary and partial differential equations this semester, thesingle most time-consuming step is the solution of a system of linear equations. We now begin a majortopic in Numerical Analysis, Numerical Linear Algebra. We assume you are familiar with the concepts ofvectors and matrices, matrix multiplication, vector dot products, and the theory of the solution of systemsof linear equations, the matrix inverse, eigenvalues and so on. We now concentrate on the practical aspectsof carrying out the computations required to do linear algebra in a numerical setting.This lab will take three sessions. If you print this lab, you may prefer to use the pdf version.We begin by noting several special test matrices that we will refer to when trying out algorithms. Thenwe state the linear system problem and consider three methods of solution, using the determinant, the inversematrix, or Gauß factorization. We will find that factorization into three simply-solved factors is the bestway to go, and we will write a m-file to perform the factorization and another to solve systems given thefactorization. We will end up with an example using our m-files as part of the numerical solution of a partialdifferential equation.2 Test matricesWe will now consider a few simple test matrices that we can use to study the behavior of the algorithmsfor solution of a linear system. We can use such problems as a quick test that we haven’t made a seriousprogramming error, but more importantly, we also use them to study the accuracy of the algorithms. Wewant test problems which can be made any size, and for which the exact values of the determinant, inverse,or eigenvalues can be determined. Matlab itself offers a “rogues’ gallery” of special test matrices throughthe Matlab function gallery, and it offers several special matrices. We will use our own functions for somematrices and Matlab’s functions for others.Some test matrices we will be using are:• The second difference matrix:• The Frank matrix,• The Hilbert matrix.• The Pascal matrix.1• The magic square matrix.2.1 The second difference matrixThe second difference matrix arises in approximating the second derivative on equally spaced data. Youhave seen this matrix and calculated its eigenvalues, eigenvectors, and determinant in Lab 5. The formulais:Ai,j=2 for j = i−1 for |j − i| = 10 for |j − i| > 1So that the matrix looks like:2 −1−1 2 −1−1 2 −1...−1 2 −1−1 2where blanks represent zero values. As you have seen, the matrix is positive definite, symmetric, andtridiagonal. The determinant is ±(n + 1) where n denotes the size (A is n × n). Matlab provides a functionto generate this matrix called gallery(’tridiag’,n), but we will use the m-file dif2.m.Do not use the gallery(’tridiag’,n) form in these labs except when instructed to use it. Matlabhas two forms of storage for matrices, “sparse” and “full.” Everything we have done so far has used the“full” form, but gallery(’tridiag’,n) returns the “sparse” form. The “sparse” matrix form requires somespecial handling.Download the m-file, dif2.m now.2.2 The Frank matrixRow i of the n × n Frank matrix has the formula:Ai,j=0 for j < i − 2n + 1 − i for j = i − 1n + 1 − j for j ≥ iThe Frank matrix for n = 5 looks like:5 4 3 2 14 4 3 2 10 3 3 2 10 0 2 2 10 0 0 1 1The determinant of the Frank matrix is 1, but is difficult to compute numerically, because of roundoff errors.This matrix has a special form called Hessenberg form wherein all elements below the first subdiagonal arezero. Matlab provides the Frank matrix in its “gallery” of matrices, gallery(’frank’,n), but we will usean m-file frank.m. It turns out that the inverse of the Frank matrix also has a simple formula.Download the m-file frank.m now.22.3 The Hilbert matrixThe Hilbert matrix is related to the interpolation problem on the interval [0,1]. The matrix is given by theformula Ai,j= 1/(i + j − 1). For example, with n = 5 the Hilbert matrix is:11121314151213141516131415161714151617181516171819The Hilbert matrix arises in interpolation and approximation contexts because it happens that Ai,j=∫10(xi−1)(xj−1)dx. The Hilbert matrix is at once nice because its inverse has integer elements and also notnice because it is extremely difficult to compute the inverse using the usual formulæ for matrix inverses.Matlab has special functions for the Hilbert matrix and its inverse, called hilb(n) and invhilb(n), andwe will be using these functions in this lab.2.4 The Pascal matrixThe Pascal matrix is generated by a form of Pascal’s triangle. It is defined as Pi,j=(i+j−2j−1)(the binomialcoefficient). For example, with n = 5, the Pascal matrix isP =1 1 1 1 11 2 3 4 51 3 6 10 151 4 10 20 351 5 15 35 70As you can tell, the colored portion of this matrix is a Pascal triangle, written on its side.A related lower triangular matrix can be written (for n = 5) asL =1 0 0 0 01 −1 0 0 01 −2 1 0 01 −3 3 −1 01 −4 6 −4 1Except for signs, this matrix is a Pascal triangle written asymmetrically. The matrixL1=1 0 0 0 01 1 0 0 01 2 1 0 01 3 3 1 01 4 6 4 1is truly a Pascal triangle written asymmetrically. Matlab provides a function pascal so that P=pascal(5),L=pascal(5,1), and L1=abs(L).The Pascal matrix, like the Frank matrix, is very ill-conditioned for large values of n. In addition, thePascal matrix satisfies the conditions3• det(P ) = det(L) = det(L1) = 1.• L2= I (I is the identity matrix)• P = LLT= L1LT1More information is available in a Wikipedia article http://en.wikipedia.org/wiki/Pascal matrix3 The linear system “problem”The linear system “problem” can be posed in the following way. Find an n-vector x that satisfies the matrixequationAx = b (1)where