Transcription of Solving a tridiagonal linear system
1 Solving a tridiagonal linear systemMATH2071: Numerical Methods in Scientific Computing jburkardt/classes/math20712020/ trididiagonal linear system has three nonzero stripes . tridiagonal SolutionEfficiently store and solve a tridiagonal system of linear A tridiagonal linear systemA linear systemAx=bis calledtridiagonalif, in thei-th equation, only the coefficientsai,i 1,ai,iandai,i+1are nonzero. The first and last equations will actually only have two nonzero coefficients. The name comesfrom the fact that, if we display the matrix, the nonzero entries fall along three diagonals:a1,1a1,2a2,1a2,2a2,3a3,2a3,3a3 , an 1,n 2an 1,n 1an 1, ,n 1an,nThis special form of a linear system is of interest to us for two reasons: the nonzero matrix entries can be stored using 3nspace rather thann2; the linear system can be solved usingO(n) operations rather thanO(n3).
2 Thus, if we recognize when we are dealing with a tridiagonal system , we can greatly reduce the necessarystorage and computational effort needed to obtain a tridiagonal linear system is one of the simplest examples of a sparse matrix. We will be able to usethis approach to solve the linear system associated with a discretized version of the one-dimensional How to solve a tridiagonal systemConsider the following tridiagonal linear system :3x1+4x2= 116x1+10x2+1x3= 2910x2+9x3= 47Of course, Gauss elimination will work on this problem, but we should notice that some of our work hasbeen done for us already.
3 In column 1, there is already a zero in the (3,1) position. For a larger tridiagonalsystem, we would see a much larger savings, since all but one of the subdiagonal entries are already solve this system , we begin by noticing that if we subtract 2 times equation #1 from equation #2, weeliminatex1from equation #2:3x1+4x2= 11+2x2+1x3= 710x2+9x3= 47 Now subtract 5 times equation #2 from equation #3, to eliminatex2from equation #3:3x1+4x2= 11+2x2+1x3= 7+4x3= 12 The original tridiagonal equations now have a format known asupper triangular. Upper triangular systemsare easy to solve. Now we do the back solve 12/4x2= (7 1x3)/2x1= (11 4x2)/3 Notice that the forward elimination required justn 1 simple operations, and the backward substitutionrequirednslightly more complicated for special cases where we encounter a zero pivot, any tridiagonal linear system can be solved thisway.
4 We sweep down the equations, eliminating variableifrom equationi+ 1. Then we sweep upwards, Solving for variablen, thenn 1, .., until we reach variable 1, and the system has been tridiagonal solution algorithmHere is an algorithm to solveAx=bwhenAis tridiagonal :Algorithm 1 tridiagonal linear system solverfor1 j n 1doifAj,j= , the algorithm has ifs Aj+1,j/Aj, +1,j+1 Aj+1,j+1 sAj,j+ +1 bj+1 forxn bn/An,nforn 1 j 1doxj (bj A(j, j+ 1)xj+1)/Aj,jend for2 Look at this algorithm, and consider how we are justified in saying that the solution of a tridiagonal systemrequiresO(n) arithmetic reason for this reduction in work is not just that a tridiagonal matrix has a special structure, but thatour algorithm realizes this structure, and takes advantage of it by only doing the work that it knows has tobe MATLAB tridiagonal solver1functionx = t r i d i a gs o l v e ( A, b )2n =s i z e( A, 1 ) ;3f o rj = 1 : n 14s = A( j +1, j ) / A( j , j ) ;5A( j +1, j +1) = A( j +1, j +1) s A( j , j +1) ;6b( j +1) = b( j +1) s b( j ) ;7end89x =zeros( n , 1 ) ;10x(n) = b(n) / A(n , n)11f o rj = n 1 : 1 : 112x( j ) = ( b( j ) A( j , j +1) x( j +1) ) / A( j , j ).
5 13end1415return16endListing 1: A compact tridiagonal storage schemeSince a tridiagonal matrix only has 3n 2 nonzero entries, it may be worth abandoning the fulln nstoragescheme we are used to. Instead, we can create 3 vectors,a, b, c, each of lengthn, and then, for each rowi,storingai=Ai,i 1,bi=Ai,i,ci=Ai,i+1. We will leavea1andcnto zero, since they don t correspond tolegal entries we write our solution algorithm, we just have to remember to replace each mention of an entry ofAby the corresponding element ofa,borc. Also, the right hand side will get renamed tof!ab c3x1+4x2= 1103 46x1+10x2+1x3= 29 6 10 110x2+9x3= 47109 0If we are very careful, we can figure out how to replace references to the full matrixAby references tothe subdiagonal, diagonal, and superdiagonal vectorsa, b, c.
6 This will allow us to create a new functiontridiagsparsesolve()which carries out Gauss elimination on the sparse version of our linear MATLAB sparse tridiagonal solverWe can use the following table of replacements to create a sparse storage version of our solver:Aj,j 1 ajAj,j bjAj,j+1 cjAj+1,j aj+1Aj+1,j+1 bj+1Aj+1,j+2 cj+11functionx = t r i d i a gs p a r s es o l v e ( a , b , c , f )23n =length( a ) ;45f o rj = 1 : n 16s = a ( j +1) / b( j ) ;7b( j +1) = b( j +1) s c ( j ) ;8f ( j +1) = f ( j +1) s f ( j ) ;9end1011x =zeros( n , 1 ) ;12x(n) = f (n) / b(n) ;13f o rj = n 1 : 1 : 114x( j ) = ( f ( j ) c ( j ) x( j +1) ) / b( j ) ;15end1617return18endListing 2: tridiagonal solver now takes advantage of both reduced space requirements,O(3n) and reduced workrequirementsO(n).
7 7 Exercise: a tridiagonal matrix associated with splinesIf we have a grid ofnpointsx, then we can define ann 1 vectorgsuch thatgi=xi+1 xi. These spacingscan be used to define ann ntridiagonal matrix associated with interpolation by cubic splines. ThefunctionsA=splinematrix(g)and[a,b,c]= splinesparsematrix(g)return dense and sparse versionsof this order to manufacture a problem, we choose a simple solution vectorxand computef=A x. However,for the sparse matrix version, we can callf=tridiagsparsemv(a,b,c,x)to compute this , to solve the linear system , we can callx=tridiagsolve(A,f)orx=tridiagsparse solve(a,b,c,f).
8 1. Decide on a sizen2. Set the vectorgton 1 random positive Compute the spline matrixAor spline sparse matrix [a, b, c].4. Decide on a solution Compute the right hand sidefby multiplying the matrix Call tridiagonal solver or tridiagonal sparse solver to Verifyx1 andx2 are equal .. actually, that||x1 x2||is Challenge: Computing a tridiagonal inverseStarting with one of the tridiagonal solver codes, make a new copy calledX = tridiaginverse(A)orX= tridiagsparseinverse(a,b,c). Instead of inputting a right hand side vectorf, set up a densen nmatrixF, which is initialized to the identity matrix.
9 Now solve for a densen nmatrixX, which will bethe inverse of the tridiagonal the original code referenced the vectorsf(j) orx(j), you will need to replace this by a referenceto an entire row,F(j,:) orX(j,:). This should be enough to get your code to your code to compute the inverse of the spline matrixAfrom the previous would be nice if the inverse of a tridiagonal matrix was also a tridiagonal matrix. It isn t!Verify that your output matrixXis the correct inverse of A by computingA X. If you worked with thesparse form, compute the product byAX = tridiagsparsemv(a,b,c,X).5
