Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 3.	The Solution of Linear Systems  AX = B
	Algorithm 3.1	Back Substitution	
	Algorithm 3.2	Upper-Triangularization Followed by Back Substitution	
	Algorithm 3.3	PA = LU Factorization with Pivoting	
	Algorithm 3.4	Jacobi Iteration	
	Algorithm 3.5	Gauss-Seidel Iteration

Algorithm 3.1 (Back Substitution). To solve the upper-triangular system:
	a1,1 x1	+	a1,2 x2	+ ... +	a1,n-1 xn-1	+	a1,n xn	 =  b1
		           	a2,2 x2	+ ... +	a2,n-1 xn-1	+	a2,n xn	 =  b2
					                                    :		              :		        :
				                          an-1,n-1 xn-1  + an-1,n  xn	=	bn-1
						                                               an,n  xn	=  bn
Proceed with the method only if all the diagonal elements are non-zero.   
First compute  xn = bn/an,n  and then use the rule

xr   =  ÊÊ(br - ar,r+1 xr+1 - ... - a2,n-1 xn-1 - a2,n xn)/ar,rÊ      
 
for  r = n-1,n-2, ...,1.

Remark.  Additional steps have been placed in the algorithm which 
perform the side calculation of finding the determinant of  A. 
Section	3.3	Upper-Triangular Linear Systems   Page  145

Off[General::spell1];
Clear[BackSub,a0,b0];
BackSub[a0_,b0_] :=
Module[{a,b,x,k,sum,mna,nb,ma,na},
  a = a0;
  b = b0;
  mna = Dimensions[a];
  nb = Length[b];
  ma = mna[[1]];
  na = mna[[2]];
  det = 1;
  For[k=1,k<=ma,k++,
    det = det a[[k,k]]; ];
  If[ma==na && det!=0,
    Module[{i,j},
    x = Table[0,{i,ma}];
    x[[ma]] = b[[ma]]/a[[ma,ma]];
    Clear[i,j,k];
    For[i=ma-1,1<=i,i--,
      sum = 0;
      For[j=i+1,j<=ma,j++,
        sum = sum + a[[i,j]] x[[j]] ]; 
      x[[i]] = (b[[i]] - sum)/a[[i,i]]; ]; ],
    x = MatrixIsSingular];
  Return[x]; ];
On[General::spell1];

Chapter 3, Example 3.12, Page 144.

Use back substitution to solve the upper-triangular linear system  AX = B:

A = {{3,-1,2,3},{0,-2,7,-4},{0,0,6,5},{0,0,0,3}};
A = N[A,19];
TableForm[A]

B = {20,-7,4,6};
B = N[B,19];
TableForm[B]

X = BackSub[A,B];
TableForm[X]

Verify by matrix multiplication to see if  AX - B  =  0.

Y = A.X-B;
TableForm[Y]

Compare with Mathematica's built in routine.

Z = LinearSolve[A,B];
TableForm[Z]
TableForm[A.Z-B]

Algorithm 3.2 (Upper-Triangularization Followed by Back Sub.)
To construct the solution to  AX = B, by first reducing the augmented
matrix [A,B] to upper triangular form and then performing back
substitution.
Section	3.4	Gaussian Elimination and Pivoting   Page  156

Let the coefficients of A and the constants of B be stored in the augmented matrix  [A,B] = (ai,j)  which has n rows and n+1 columns; that is, the column vector B = (ai,n+1) is stored in column n+1 of the augmented matrix [A,B]. The system can be written

 a1,1 x1 +   a1,2 x2 +...+ a1,n-1 xn-1 + a1,n xn  =  a1,n+1

 a2,1 x1 + a2,2 x2 +...+ a2,n-1 xn-1 +  a2,n xn  =  a2,n+1
       .          .                   .            .      .

       :          :                   :            :      :
an-1,1 x1 + an-1,2 x2 +...+ an-1,n-1 xn-1 + an-1,n xn  =  an-1,n+1

an,1 x1 +   an,2 x2 + ... +   an,n-1 xn-1 +   an,n xn  =  an,n+1

Row operations will be used to eliminate  xp  in column  p.
REMARKS. Rather than actually interchange rows, row interchanges are recorded in the pointer vector Row(p). If either the matrix A or vector B is to be used later, then a copy must be stored someplace else besides the augmented matrix.

Off[General::spell1];
Clear[GaussJordan,a0];
GaussJordan[a0_] :=
Module[{a,i,j,mna,n,n1,p,row,x,
        cond,Nonsingular,Singular},
  a = a0;
  mna = Dimensions[a];
  n  = mna[[1]];
  n1 = mna[[2]];
  x = Table[0,{i,n}];
  cond = Nonsingular;
  Clear[row];
  For[j=1,j<=n,j++, row[j]=j];
  For[p=1,p<=n-1,p++,
  Block[{k},
    For[k=p+1,k<=n,k++,
      If[Abs[a[[row[k],p]]] > Abs[a[[row[p],p]]],
        t=row[p]; row[p]=row[k]; row[k]=t ]; ];
    If[a[[row[p],p]] == 0,
      Print["A(",row[p],",",p,") = ",a[[row[p],p]]];
      cond = Singular;
      Print["The matrix is singular."];
      Goto[getout]; ];
    Clear[k];
    For[k=p+1,k<=n,k++,
      Module[{c},
      m = a[[row[k],p]]/a[[row[p],p]];
      For[c=p+1,c<=n+1,c++,
        a[[row[k],c]]=a[[row[k],c]]-
        m a[[row[p],c]] ]; ]; ]; ]; ];
 If[a[[row[n],n]] == 0,
     Print["A(",row[n],",",n,") = ",a[[row[n],n]]];
     cond = Singular;
     Print["The matrix is singular."];
     Goto[getout]; ]; 
 x[[n]] = a[[row[n],n+1]]/a[[row[n],n]];
 Clear[k];
 For[k=n-1,1<=k,k--,
   sum:= 0;
   For[c=k+1,c<=n,c++,
     sum = sum + a[[row[k],c]] x[[c]];
     x[[k]] = (a[[row[k],n+1]]-sum)/a[[row[k],k]]; ]; ];
 Label[getout];
 Return[x]; ];
On[General::spell1];

Chapter 3, Exercise 8, Page 158.

Use upper-triangularization followed by back-substitution to solve the linear system  AX = B:

A = {{4,8,4,0},{1,5,4,-3},
     {1,4,7,2},{1,3,0,-2}};
A = N[A,19];
TableForm[A]

B = {8,-4,10,-4};
B = N[B,19];
TableForm[B]

AB = {{4,8,4,0,8},{1,5,4,-3,-4},
     {1,4,7,2,10},{1,3,0,-2,-4}};
AB = N[AB,19];
TableForm[AB]

X = GaussJordan[AB];
TableForm[X]

Verify by matrix multiplication to see if  AX - B  =  0.

Y = A.X - B;
TableForm[Y]

Compare with Mathematica's built in routine.

Z = LinearSolve[A,B];
TableForm[Z]
TableForm[A.Z-B]

Algorithm 3.3 (PA = LU Factorization with Pivoting). 
To construct the solution to the linear system AX = B  by performing
the steps:
1.	Find a permutation matrix  P , lower-triangular matrix  L   and 
		 upper-triangular matrix  U  which satisfy   PA = LU.
2.	Compute  PB  and form the equivalent linear system   LUX = PB.
3.	Solve the lower-triangular system	LY = PB  for  Y.
4.	Solve the upper-triangular system	UX = Y  for  X.
Remarks: Since the diagonal elements of L are all 1's these values do not need to be stored.  The coefficients of L below the main diagonal and the nonzero coefficients of U overwrite the matrix A. The permutation matrix P is not actually needed.  The pointer vector row serves in its place to keep track of row interchanges.  If matrix A is to be used later, then a copy must be stored someplace else.  The value  det  is the determinant 
of A.
Section	3.6	Triangular Factorizaton   Page  175

First, execute the following group of cells:

Off[General::spell1];
Clear[LUfactor,a0,row];
LUfactor[a0_] :=
Module[{a,i,j,p,mna,n,m,rowk,rowp},
  a = a0;
  mna = Dimensions[a];
  n  = mna[[1]];
  m = mna[[2]];
  det = 1;
  For[j=1,j<=n,j++, row[j]=j];
  For[p=1,p<=n-1,p++,
  Block[{c,k},
    For[k=p+1,k<=n,k++,
      If[Abs[a[[row[k],p]]]>Abs[a[[row[p],p]]],
        t=row[p]; row[p]=row[k]; row[k]=t; det=-det]; ];
    det = det a[[row[p],p]];
    If[det == 0,
      Print["The matrix is singular."];
      Goto[getout]; ];
    Clear[k];
    For[k=p+1,k<=n,k++,
      rowk = row[k];
      rowp = row[p];
      a[[rowk,p]] = a[[rowk,p]]/a[[rowp,p]];
      For[c=p+1,c<=n,c++,
         a[[rowk,c]] = a[[rowk,c]] - 
                       a[[rowk,p]] a[[rowp,c]];];];];];
  Label[getout];
  det = det a[[row[n],n]];
  Return[a]; ];

Clear[SolveLU,lu0,b0];
SolveLU[lu0_,b0_] :=
Module[{a,c,i,j,k,n,p,x,rowk,sum},
  a = lu0;
  b = b0;
  n = Length[a];
  x = Table[0,{i,n}];
  For[k=1,k<=n,k++,
    If[a[[row[k],k]]==0, Goto[getout] ]; ];
  x[[1]] = b[[row[1]]];
  Clear[k];
  For[k=2,k<=n,k++,
    sum = 0;
    rowk = row[k];
    For[c=1,c<=k-1,c++,
      sum = sum + a[[rowk,c]] x[[c]] ];
    x[[k]] = b[[rowk]] - sum; ];
  x[[n]] = x[[n]]/a[[row[n],n]];
  Clear[k];
  For[k=n-1,1<=k,k--,
    sum = 0;
    rowk = row[k];
    For[c=k+1,c<=n,c++,
      sum = sum + a[[rowk,c]] x[[c]] ];
    x[[k]] = (x[[k]] - sum)/a[[rowk,k]]; ];
  Label[getout];
  Return[x]; ];
On[General::spell1];

Chapter 3, Exercise 7, Page 178.

Use LU factorization with pivoting to solve the linear system  AX = B:

A = {{2,4,-4,0},
     {1,5,-5,-3},
     {2,3,1,3},
     {1,4,-2,2}};
A = N[A,19];
TableForm[A]

LU = LUfactor[A];
TableForm[LU]

B = {12,18,8,8};
B = N[B,19];
TableForm[B]

X = SolveLU[LU,B];
TableForm[X]

Verify by matrix multiplication to see if  AX - B  =  0.

Y = A.X - B;
TableForm[Y]

Compare with Mathematica's built in routine.

Z = LinearSolve[A,B];
TableForm[Z]
TableForm[A.Z-B]

Algorithm 3.4 (Jacobi Iteration). To solve the linear system AX = B  by starting with  P0 = 0, and generating a sequence {Pk} which converges to the solution  P,  i.e.  AP = B.   A sufficient condition for the method to be applicable is that  A  is diagonally  dominant. 
Section	3.7	Iterative Methods for Linear Systems   Page  186

Off[General::spell1];
Clear[Jacobi,A,B,V,P0,Tol,max];
Jacobi[A_,B_,P0_,Tol_,max_] := 
Block[{r,c,k,n,cond,stat,X,Xold,row,
       DiagonallyDominant,NotDominant},
  norm[V_] := Sqrt[V^2/.List->Plus]//N;
  sep = 1;
  k = 0;
  cond = 1;
  stat = DiagonallyDominant;
  X = P0 //N;
  Xold = P0;
  n = Length[X];
  For[r=1,r<=n,r++,
    row = 0;
    For[c=1,c<=n,c++,
      row = row + Abs[A[[r,c]]]//N];
    If[row >= 2 Abs[A[[r,r]]], cond = 0] ];
  If[cond == 0, stat=NotDominant; Goto[getout]];
  While[k<max && sep > Tol,
    For[r=1,r<=n,r++,
      sum = B[[r]]//N;
      For[c=1,c<=n,c++,
        If[c!=r,sum = sum-A[[r,c]]Xold[[c]]//N];];
      X[[r]] = sum/A[[r,r]]//N ];
    sep = norm[X - Xold];
    Print[N[X,14]];
    Xold = X;
    k = k + 1];
  Label[getout];
  Return[N[X,14]]; ];
On[General::spell1];

Chapter 3, Example 3.25, Page 180.

Use Jacobi iteration to solve the linear system  AX = B:

A = {{4,-1,1},
     {4,-8,1},
     {-2,1,5}};
A = N[A,19];
TableForm[A]

B = {7,-21,15};
B = N[B,19];
TableForm[B]

P0 = {1,2,2};
TableForm[P0]

X = Jacobi[A,B,P0,0.00000001,100]
TableForm[X]

Verify by matrix multiplication to see if  AX - B  =  0.

Y = A.X - B

Compare with Mathematica's built in routine.

Z = LinearSolve[A,B];
TableForm[Z]
TableForm[A.Z-B]

Algorithm 3.5 (Gauss-Seidel Iteration). To solve the linear system AX = B by starting with P0 = 0, and generating a sequence  {Pk}  which converges to the solution  P,  i.e.  AP  = B.   A sufficient condition for the method to be applicable is that  A  is diagonally  dominant. 
Section	3.7	Iterative Methods for Linear Systems   Page  187

Off[General::spell1];
Clear[GaussSeidel,A,B,V,P0,Tol,max];
GaussSeidel[A_,B_,P0_,Tol_,max_] := 
Block[{r,c,k,n,X,Xold,cond,row,stat,
       DiagonallyDominant,NotDominant},
  norm[V_] := Sqrt[V^2/.List->Plus]//N;
  sep = 1;
  k = 0;
  cond = 1;
  stat = DiagonallyDominant;
  X = P0 //N;
  Xold = P0;
  n = Length[X];
  For[r=1,r<=n,r++,
    row = 0;
    For[c=1,c<=n,c++,
          row = row + Abs[A[[r,c]]]//N];
    If[row >= 2 Abs[A[[r,r]]], cond = 0] ];
  If[cond == 0, stat=NotDominant; Goto[getout]];
  While[k<max && sep > Tol,
    For[r=1,r<=n,r++,
      sum = B[[r]]//N;
      For[c=1,c<=n,c++,
        If[c!=r, sum = sum-A[[r,c]]X[[c]]//N];];
      X[[r]] = sum/A[[r,r]]//N ];
    sep = norm[X - Xold];
    Print[N[X,14]];
    Xold = X;
    k = k + 1];
  Label[getout];
  Return[N[X,14]]; ];
On[General::spell1];

Chapter 3, Example 3.25, Page 180.

Use Gauss-Seidel iteration to solve the linear system:

A = {{4,-1,1},
     {4,-8,1},
     {-2,1,5}};
A = N[A,19];
TableForm[A]

B = {7,-21,15};
B = N[B,19];
TableForm[B]

P0 = {1,2,2};
TableForm[P0]

X = GaussSeidel[A,B,P0,0.00000001,100]
TableForm[X]

Verify by matrix multiplication to see if  AX - B  =  0.

Y = A.X - B;
TableForm[Y]

Compare with Mathematica's built in routine.

Z = LinearSolve[A,B];
TableForm[Z]
TableForm[A.Z-B]