Gaussian Elimination with Partial Pivoting

Elimination with partial pivoting and LU decomposition

The Gaussian elimination can be looked at as a two-step algorithm: the first step is an  LU decomposition,  QA=LU, where Q is a permutation matrix,  L is a unit lower triangular matrix that contains the multipliers, and U is an upper triangular matrix. Following the book, we consider Gaussian elimination with scaled row pivoting. In our algorithm, we keep the permutation vecrot in P (the matrix Q = ((Subscript[Q, i,j]))_(i,j=1)^n is defined by (Subscript[Q, i,j]) = (Subscript[δ, (Subscript[p, i]),j])). We  use "digit" in the algorithm so that every step is computed up to the given digits; the default value is set as digit =10.

Algorithm

GaussLU[a0_] := 
Module[{a,i,j,k,n,x,z, digit}, digit=10; 
  a = a0; n = Length[a]; x=Table[a[[i,j]],{i,1,n},{j,1,n}];
    For[i=1,i<=n,i++, p[i]=i];
    For[i=1,i<=n,i++, 
       s[i]=Max[Table[Abs[a[[i,j]]],{j,1,n}]] ];
     For[k=1,k<=n-1,k++, Block[{l,i,j},
      For[j=k+1,j<=n,j++,
         If[Abs[x[[p[j],k]]]/s[p[j]]   > Abs[x[[p[k],k]]]/s[p[k]], t=p[k]; p[k]=p[j]; p[j]=t ]; ];
      For[i=k+1,i<=n,i++,
         z = N[x[[p[i],k]]/x[[p[k], k]],digit]; x[[p[i],k]]=z; 
          For[l=k+1,l<=n,l++,
           x[[p[i],l]]=N[x[[p[i],l]]- z x[[p[k],l]],digit];    ];];];];
Clear[P,Q,L,U,T];     
    P=Table[p[i],{i,1,n}]; 
    Q=Table[0,{i,1,n},{j,1,n}];For[i=1,i<=n,i++, Q[[i,p[i]]]=1]; T=Q.x; 
    L=Table[If[i !=j,  If[i > j, T[[i,j]], 0],1],{i,n},{j,n}];  
    U=Table[If[i <=j, T[[i,j]], 0],{i,n},{j,n}];     
Return[x]; ];

The function of the algorithm is of the form  GaussLU[matrix],  the output is a square matrix, the lower trianglular is the L matrix without diagonal and the upper triangle, together with the diagonal, is the U matrix of the LU decomposition of the input matrix.

Example

Examle 1:  Here is a matrix of 4 by 4

A = {{2,4,-4,0},{1,5,-5,-3},{2,3,1,3},{1,4,-2,2}}; 
MatrixForm[A];

LU = GaussLU[A];
MatrixForm[LU]

The permutation verctor in the parital pivoting is

P

You can also see the matrix L and U explicitly by typing

MatrixForm[L]  
MatrixForm[U]

Solving Equations with Gaussian Elimination

After QLU decomposition, the elimination step has been carried out and we can solve the linear system of equations Ax = b
rather easily. We again use "digit" in the algorithm so that every step is computed up to the given digits; the default value is set as digit =10.

Algorithm

GaussSol[a0_,b0_] := Module[{a,b,k,n,x, digit}, digit=10; 
  a = a0;  b = b0; n = Length[b]; x = Table[0,{i, n}];
  For[k=1,k<=n-1,k++,
    For[i=k+1,i<=n,i++,
      b[[p[i]]] =N[ b[[p[i]]] - a[[p[i],k]] b[[p[k]]],digit] ];];
  For[i=n,1<=i,i--,
     x[[i]] = N[(b[[p[i]]] -  Sum[a[[p[i],j]] x[[j]],{j, i+1, n}])/a[[p[i],i]],digit]; ];
  Return[x]; ];

The function of the algorithm is of the form  GaussSol[A, B], where A is the matrix and B is the vector in the equation Ax =B. The output is the solution vecote x.

Example

Examle 1:  Here is a system of 4 by 4

A = {{2,4,-4,0},{1,5,-5,-3},{2,3,1,3},{1,4,-2,2}};
B = {12,18,8,8};
StringForm["A= ``    B = ``", MatrixForm[A], MatrixForm[B]]

LU=GaussLU[A];
GaussSol[LU,B]

Solving Tridiagonal System

If the matrix of the linear system of equations Ax = b is of the tridiagonal form 
                                         (({{({{(Subscript[d, 0]), (Subscript[c, 0])}, {(Subscript[a, 1]), (Subscript[d, 1])}}), ({{ ,  }, {(Subscript[c, 1]),  }})}, {({{ ,  (Subscript[a, 1])}, { ,  }}),  ({{(Subscript[d, 2]), ...}, {..., ...}})}}))
then the linear system of equtions can be solved easily with a forword substitution followed by a backward substitution. We again use "digit" in the algorithm so that every step is computed up to the given digits; the default value is set as digit =10.

Algorithm

TriDiagonal[a0_,d0_,c0_,b0_] :=Module[{i,xmult, digit}, digit =10;
  a = a0; d=d0; c=c0; b=b0;   n = Length[b];  x=Table[0,{i,n}];
  For[i=2,i<=n,i++,    xmult = a[[i-1]]/d[[i-1]];  d[[i]] = d[[i]]-xmult*c[[i-1]];
    b[[i]] = b[[i]]-xmult*b[[i-1]];];  x[[n]] = N[b[[n]]/d[[n]],digit];
  For[i=n-1,1<=i,i--,     x[[i]] = N[(b[[i]] - c[[i]]*  x[[i+1]])/d[[i]],digit];];
 Return[x]; ];

The function of the algorithm is of the form  TriDiagonal[a, d, c,  B], where {a}, {b}, {c},  are elements of the tridiagonal matrix and B is the right hand side of Ax = b.

Example

Examle 1:  Here is a system of 4 by 4

a = {4,-4,0}; d={1,5,-5,-3}; c={2,3,1}; b={1,4,-2,2};
A= Table[Switch[i-j,-1,c[[i]],0, d[[i]],1,a[[i-1]],_,0], {i,4},{j,4}];
MatrixForm[A]

TriDiagonal[a ,d ,c ,b ]

Examle 2:  Here is a system of 6 by 6

a = {-1,-1,-1,-1,-1}; d={2,2,2,2,2,2}; c={-1,-1,-1,-1,-1}; b={1,2,3,4,5,6} ;
A= Table[Switch[i-j,-1,c[[i]],0,d[[i]],1,a[[i-1]],_,0], {i,6},{j,6}];
MatrixForm[A]

TriDiagonal[a ,d ,c ,b ]