Engineering Analysis ENG 3420 Fall 2009Lecture 14PivotingTridiagonal systems of linear equationsTridiagonal system solverEngineering Analysis ENG 3420 Fall 2009Dan C. MarinescuOffice: HEC 439 BOffice hours: Tu-Th 11:00-12:0022Lecture 14Lecture 14 Last time: Solving systems of linear equations (Chapter 9) Graphical methods Cramer’s rule Gauss elimination Today: Discussion of pivoting Tri-diagonal system solver Examples Next Time LU Factorization (Chapter 10)function x=GaussNaive(A,b)ExA=[A b];[m,n]=size(A);q=size(b);if (m~=n) fprintf ('Error: input matrix is not square; n = %3.0f, m=%3.0f \n', n,m);Endif (n~=q) fprintf ('Error: vector b has a different dimension than n; q = %2.0f \n', q);endn1=n+1;for k=1:n-1for i=k+1:n factor=ExA(i,k)/ExA(k,k); ExA(i,k:n1)= ExA(i,k:n1)-factor*ExA(k,k:n1); EndEndx=zeros(n,1);x(n)=ExA(n,n1)/ExA(n,n);for i=n-1:-1:1 x(i) = (ExA(i,n1)-ExA(i,i+1:n)*x(i+1:n))/ExA(i,i);end>> A=[1 1 1 0 0 0 ; 0 -1 0 1 -1 0 ; 0 0 -1 0 0 1 ; 0 0 0 0 1 -1 ; 0 10 -10 0 -15 -5 ; 5 -10 0 -20 0 0]A = 1 1 1 0 0 0 0 -1 0 1 -1 0 0 0 -1 0 0 -1 0 0 0 0 1 -1 0 10 -10 0 -15 -5 5 -10 0 -20 0 0b =[0 0 0 0 0 200]>> b=b'b =00000200>> x = GaussNaive(A,b) x =NaNNaNNaNNaNNaNNaNPivoting If a coefficient along the diagonal is 0 (problem: division by 0) or close to 0 (problem: round-off error) then the Gauss elimination causes problems. Partial pivoting Æ determine the coefficient with the largest absolute value in the column below the pivot element. The rows can then be switched so that the largest element is the pivot element. Complete pivoting Æ check also the rows to the right of the pivot element are also checked and switch columns.function x=GaussPartialPivot(A,b)ExtendedA=[A b];[m,n]=size(A);q=size(b);if (m~=n) fprintf ('Error: input matrix is not square; n = %3.0f, m=%3.0f \n', n,m);Endif (n~=q) fprintf ('Error: vector b has a different dimension than n; q = %2.0f \n', q);endn1=n+1;for k=1:n-1[largest,i]=max(abs(ExtendedA(k:n,k)));nrow = i +k -1;if nrow ~=kExtendedA([k,nrow],:) = ExtendedA([nrow,k],:);endendfor k=1:n-1for i=k+1:n factor=ExtendedA(i,k)/ExtendedA(k,k); ExtendedA(i,k:n1)= ExtendedA(i,k:n1)-factor*ExtendedA(k,k:n1); EndEndx=zeros(n,1);x(n)=ExtendedA(n,n1)/ExtendedA(n,n);for i=n-1:-1:1 x(i) = (ExtendedA(i,n1)-ExtendedA(i,i+1:n)*x(i+1:n))/ExtendedA(i,i);endA = 1 1 1 0 0 0 0 -1 0 1 -1 0 0 0 -1 0 0 -1 0 0 0 0 1 -1 0 10 -10 0 -15 -5 5 -10 0 -20 0 0>> k=4; n=6; A(k:n,k)ans =00-20>> k=4; n=6; A(2,k:6)ans =1 -1 0>> k=4; n=6; [largest,i]=max(abs(A(k:n,k))); nrow = i +k -1, largest, inrow =6largest =20i =3Tridiagonal systems of linear equations A tridiagonal system of linear equations Æ a banded system with a bandwidth of 3: Can be solved using the same method as Gauss elimination, but with much less effort because most of the matrix elements are already 0.f1g1e2f2g2e3f3g3⋅⋅⋅⋅⋅ ⋅⋅⋅ ⋅en−1fn−1gn−1enfn⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ x1x2x3⋅⋅⋅xn−1xn⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎫ ⎬ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪ =r1r2r3⋅⋅⋅rn−1rn⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎫ ⎬ ⎪ ⎪ ⎪ ⎭ ⎪ ⎪ ⎪Tridiagonal system
View Full Document