Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016
Lecture topic: Linear Algebra

Creating Matrices in Mathematica

Here is a quick way to create matrices in Mathematica when you know the form of the elements:

Clear[A,a,B,b,Id]
A=Array[a,{10,10}];
MatrixForm[A]
Table[a[i,j]=i/j^2,{i,1,10},{j,1,10}];
MatrixForm[A]

Here is the identity matrix:

Id=IdentityMatrix[10];
MatrixForm[Id]

Here is a vector:

B=Array[b,{10}];
MatrixForm[B]
Table[b[i]=w^i,{i,1,10}];
MatrixForm[B]

Here is how to join matrices to make an augmented matrix:

aug1=Join[A,Id,2];
MatrixForm[aug1]

CC=Array[c,{10,11}];
Table[c[i,j]=i/j^2,{i,1,10},{j,1,10}];
Table[c[i,11]=b[i],{i,1,10}];
MatrixForm[CC]

Clear[A,B,CC,a,b,c,Id,aug1]

Guassian Elimination

System of Equations

Let's define a matrix to represent a system of equations:

Clear[x,x1,x2,x3,x4,m,b]

x={x1,x2,x3,x4};
m={{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}};
b={4,1,-3,4};

TableForm[Table[Sum[m[[j,i]]x[[i]],{i,1,4}]==b[[j]],{j,1,4}]]

Here is the solution to the system of equations.

Solve[Table[Sum[m[[j,i]]x[[i]],{i,1,4}]==b[[j]],{j,1,4}],{x1,x2,x3,x4}]

Here is the system in matrix (tableau) form, which basically supresses the variables.

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

RowReduce

We can use Mathematica's built-in commands to row reduce this matrix, and of course we will have learned nothing.

MatrixForm[m]
RowReduce[m];
MatrixForm[%]

The above shows the solution to the system of equations, since the corresponding system of equations to the matrix above is

x1=RowReduce[m][[1,5]]
x2=RowReduce[m][[2,5]]
x3=RowReduce[m][[3,5]]
x4=RowReduce[m][[4,5]]

which is the same solution we found above when using NSolve.

The row reductions have transformed a matrix from its original form into an echelon form. This echelon form is an equivalent expression of an underlying system of equations.

RowReduce also works symbolically:

msymb={{a11,a12,a13,b1},{a21,a22,23,b2},{a31,a32,a33,b3}};
RowReduce[msymb];
MatrixForm[msymb]
MatrixForm[RowReduce[msymb]]

Echelon Form of a Matrix

A rectangular matrix is in row echelon form if it has the following properties:

All nonzero rows are above any rows of all zeros

Each leading entry of a row is in a column to the right of the leading entry of the row above it

All entries in a column below a leading entry are zeros

If a matrix in echelon form satisfies the following additional condition, then it is in reduced row echelon form:

Each leading 1 is the only nonzero entry in its column

Row Reduction Step by Step to get echelon form

Let's define each row, and then do the row reductions ourselves so we can see the process. This just defines the same system we used above, by assigning the rows from the matrix m we used earlier.

m1={R1,R2,R3,R4}=m;
MatrixForm[%]

This will do thee following row reductions:
(Subscript[R, 2])-> (Subscript[R, 2])-2 (Subscript[R, 1]) 
(Subscript[R, 3])-> (Subscript[R, 3])-3 (Subscript[R, 1]) 
(Subscript[R, 4])-> (Subscript[R, 4])+ (Subscript[R, 1])

m2={R1,R2-2R1,R3-3R1,R4+R1}; (* row reduction *)
{R1,R2,R3,R4}=%; (* assign new rows *)
MatrixForm[%]

This will do thee following row reductions:
(Subscript[R, 3])-> (Subscript[R, 3])-34(Subscript[R, 2]) 
(Subscript[R, 4])-> (Subscript[R, 4])+3 (Subscript[R, 2])

m3={R1,R2,R3-4R2,R4+3R2};
{R1,R2,R3,R4}=%;
MatrixForm[%]

Notice we didn’t get to reduced row echelon form, just row echelon form.

We have now finished the procedure known as Gaussian Elimination with Backward Substitution.

Back substitution will provide the solutions to the uderlying system of equations, since the above matrix represents a system of equations where starting from the bottom up we can determine each variable using backsubstitution.
R4: -13(Subscript[x, 4])=-13 -> (Subscript[x, 4])=1
R3: 3(Subscript[x, 3])+13(Subscript[x, 4])=13-> (Subscript[x, 3])=13-13(1)=0
R2: -(Subscript[x, 2])-(Subscript[x, 3])-5(Subscript[x, 4])=-7->(Subscript[x, 2])=7-0-5(1)=2
R1: (Subscript[x, 1])+(Subscript[x, 2])+3(Subscript[x, 4])=4->(Subscript[x, 1])=4-(2)-3(1)=-1

These are the same values we found earlier.

Row Reduction With Pivoting ((Subscript[(m^((k))), kk])=0)

m={{1,-1,2,-1,-8},{2,-2,3,-3,-20},{1,1,1,0,-2},{1,-1,4,3,4}};
MatrixForm[m]

m1={R1,R2,R3,R4}=m;
MatrixForm[m]

m2={R1,R2-2R1,R3-R1,R4-R1};
{R1,R2,R3,R4}=%;
MatrixForm[%]

At this point, we need to pivot the second row out, since the procedure cannot continue with a zero in the position (Subscript[m, 22]). The pivoting is just interchanging rows 2 and 3, (Subscript[R, 2])<->(Subscript[R, 3])

m3={R1,R3,R2,R4};
{R1,R2,R3,R4}=%;
MatrixForm[%]

m3={R1,R2,R3,R4+2R3};
{R1,R2,R3,R4}=%;
MatrixForm[%]

We have now finished the procedure.

Row Reduction With Partial Pivoting

Partial Pivoting acts to place the largest element on the diagional, to reduce the effect of rounding errors.

The pivoting procedure will help reduce rounding errors that occur when division by a small number occurs.

(1/0.0007)

(1/0.00071)

Youch! That error can bite you big time if you aren't careful.On a computer where significant number of digits are kept, you can reduce this type of error, but never eliminate it entirely. A good routine will use partial pivoting.

Example with No Pivoting

Here, I am setting the computer accuracy to be p, since I want to track round off errors.

p=8;

m={{SetAccuracy[3/10000,p],SetAccuracy[3,p],SetAccuracy[20001/10000,p]},{SetAccuracy[1,p],SetAccuracy[1,p],SetAccuracy[1,p]}};
MatrixForm[%]

m1={R1,R2}=m;
MatrixForm[%]

m2={R1,R2-R2[[1]]*R1/ R1[[1]]};
{R1,R2}=%;
MatrixForm[%]

Backsubstitution yields:

x2= R2[[3]]/R2[[2]]
x1= (R1[[3]]-R1[[2]]*x2)/R1[[1]]

Let's check this against Mathematica's result.

RowReduce[m]//MatrixForm

What happened? Well, the rounding errors bit us. If we do this procedure again, keeping more digits of accuracy,we will get something closer to the truth.

Example with Partial Pivoting

p=8;

m={{SetAccuracy[3/10000,p],SetAccuracy[3,p],SetAccuracy[20001/10000,p]},{SetAccuracy[1,p],SetAccuracy[1,p],SetAccuracy[1,p]}}

m1={R1,R2}=m;
MatrixForm[%]

The pivoting step is to always reorder the rows so that the largest coefficent in the column appears in the diagonal, and then use that row to make the other rows have zero elements.

m2={R2,R1};
{R1,R2}=%;
MatrixForm[%]

m3={R1,R2-R2[[1]]*R1/ R1[[1]]};
{R1,R2}=%;
MatrixForm[%]

Backsubstitution yields:

x2= R2[[3]]/R2[[2]]
x1= (R1[[3]]-R1[[2]]*x2)/R1[[1]]

Let's check this against Mathematica's result.

RowReduce[m]//MatrixForm

With partial pivoting we are able to reduce the effect of roundoff errors. Remember, without pivoting, we were NOT able to reproduce this result when we had eight digits of accuracy.

Examples:  No Pivoting and Pivoting as a function of Accuracy

Here we investigate the solution to the previous system of equations as we increase the accuracy we keep at each step.

res =1/3.;(*since we know the result, use this to compare error*)

$Post=. (*rests numerical evaluation, just a bit of Mathematica*)

data={};
Do[
m={{SetAccuracy[3/10000,p],SetAccuracy[3,p],SetAccuracy[20001/10000,p]},{SetAccuracy[1,p],SetAccuracy[1,p],SetAccuracy[1,p]}};
m1={R1,R2}=m;
m2={R1,R2-R2[[1]]*R1/ R1[[1]]};
{R1,R2}=m2;
x2= R2[[3]]/R2[[2]];
x1= (R1[[3]]-R1[[2]]*x2)/R1[[1]];
data=Append[data,{p,x1,x2,Abs[(x1-res)]}],{p,1,15,1}]
TableForm[data,TableHeadings->{None,{"Accuracy","x1","x2","Error"}}]

OK, so we are able to get better results by increasing the accuracy.

data={};
Do[
m={{SetAccuracy[3/10000,p],SetAccuracy[3,p],SetAccuracy[20001/10000,p]},{SetAccuracy[1,p],SetAccuracy[1,p],SetAccuracy[1,p]}};
m1={R1,R2}=m;
m2={R2,R1};
{R1,R2}=m2; (*the pivoting*)
m3={R1,R2-R2[[1]]*R1/ R1[[1]]};
{R1,R2}=m3;
x2= R2[[3]]/R2[[2]];
x1= (R1[[3]]-R1[[2]]*x2)/R1[[1]];
data=Append[data,{p,x1,x2,Abs[(x1-res)]}],{p,1,15,1}]
TableForm[data,TableHeadings->{None,{"Accuracy","x1","x2","Error"}}]

At some point, if you are working with decimals a number must be rounded. Partial pivoting reduced the effects of this round off error.

Complete Pivoting

Complete pivoting searches all remaining elements, not just down the column, to find the largest element in absolute value, and then moves that element to the pivot position.

Complete pivoting is numerically less desirable than partial pivoting.

There is also the concept of scaled partial pivoting.

Neither of these techniques sees as much use as partial pivoting.

Complete Gaussian Elimination with Partial Pivoting: A Module

Here is a module that does partial pivoting. Let's see what some of the syntax means.

The matrix we are working with:

m={{1,1,0,3,4},{2,1,-1,1,1},{3,-1,-1,2,-3},{-1,2,3,-1,4}}

The rows can be accessed as follows:

m[[1]]

m[[{1,3}]]

The elements can be accessed as follows:

m[[1,3]]

Row interchanges can be implemented by:

m
m[[{1,3}]] = m[[{3,1}]];
m

That's all we need to know to construct a module to do Gaussian Elimination! The module does complete Gaussian elimination with partial pivoting. The bit in red is the partial pivoting step.

GaussElim[A0_,n_]:=Module[ {A=A0}, 
Print[ MatrixForm[A] ]; 
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
Print[MatrixForm[A]]; ]; 
Return[A];  ]

GaussElim[m,4]

This is what the built-in Mathematica command RowReduce does:

RowReduce[m]//MatrixForm

Matrix Inversion

Let's get even more trickier, and use some built in Mathematica commands to manipulate matrices.

Here is a matrix we wish to invert:

m={{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}}
MatrixForm[m]

We form the augmented matrix, and use Gaussian Elimination on it:

maug=Join[m,IdentityMatrix[4],2]

MatrixForm[maug]

mat=GaussElim[maug,4]

The inverse of the matrix is extracted by

minv=Take[mat,{1,4},{5,8}]
MatrixForm[minv]

Let's check this out:

MatrixForm[m.minv]

What we have just done is equivalent to the Mathematica command Inverse:

Inverse[m]//MatrixForm

Matrix Determinant

I've made a new module that only reduces to an upper triangular matrix, so it does Gaussian Elimination with partial pivoting and back substitution. You can also keep track automatically of the effect the manipulations have on the determinant, since row exchanges change the determinant by -1, and if you multiply a row by λ then the deteminant is multiplied by λ.

DetMat[A0_,n_]:=Module[ {A=A0}, 
Print[ MatrixForm[A] ]; 
dete=1;
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ],{ A[[{p,k}]]  = A[[{k,p}]] ,dete=-dete};   ]; ];
dete=A[[p,p]] dete; 
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=p+1, i≤n, i=i+1, A[[i]] = A[[i]] - A[[i,p]] A[[p]]   ]; 
Print[MatrixForm[A],"       ",dete]; ]; 
Return[dete];  ]

DetMat[m,4]

This is the equivalent of Mathematica's intrinsic command Det:

Det[m]

Gaussian Elimination Modules

Test Matrix

m={{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}};
MatrixForm[%]

Id = IdentityMatrix[4];
MatrixForm[Id]

augment1=Join[m,Id,2];
MatrixForm[%]

m2={{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{3,-1,-1,2}};
MatrixForm[%]

augment2=Join[m2,Id,2];
MatrixForm[%]

Gaussian Elimination Partial Pivoting with Intermediate Output Module

GaussElim1[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]}, 
Print[ MatrixForm[A] ]; 
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]==0, {Print["Gaussian Elimination Failed"], Goto[end]}]];
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]];];];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
Print[MatrixForm[A]]; ]; 
Label[end];Return[A];  ]

GaussElim1[augment1];
MatrixForm[%]

GaussElim1[augment2];
MatrixForm[%]

Gaussian Elimination Partial Pivoting Backward Substitution with Intermediate Output Module

GaussElim2[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]}, 
Print[ MatrixForm[A] ]; 
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]==0, {Print["Gaussian Elimination Failed"], Goto[end]}]];
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=p+1, i≤n, i=i+1, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; 
Print[MatrixForm[A]]; ]; 
Label[end];Return[A];  ]

GaussElim2[augment1];
MatrixForm[%]

GaussElim2[augment2];
MatrixForm[%]

Gaussian Elimination Partial Pivoting with No Output Module

Here is a routine with no intermediate output:

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]==0, {Print["Gaussian Elimination Failed"], Goto[end]}]];
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Label[end];Return[A];  ] ;

GaussElim3[augment1];
MatrixForm[%]

GaussElim3[augment2];
MatrixForm[%]

Matrix Inversion Module

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];A3=Take[A2,{1,n},{n+1,2n}];Return[A3]];

m//MatrixForm

MatInv[m];
MatrixForm[%]
MatrixForm[%%.m]

Matrix Determinant Module Partial Pivoting

MatDet[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]}, 
dete=1;
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]==0, {Print["Gaussian Elimination Failed"], Goto[end]}]];
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ],{ A[[{p,k}]]  = A[[{k,p}]] ,dete=-dete};   ]; ];
dete=A[[p,p]] dete; 
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=p+1, i≤n, i=i+1, A[[i]] = A[[i]] - A[[i,p]] A[[p]]   ];  ]; 
Label[end];Return[dete];  ]

MatDet[m];
MatrixForm[%]
Det[m]

Solving Linear Systems of Equations

We can solve a linear system of equations in a variety of ways.

A x = b
A is an n × n matrix,
x, b  are n vectors.

The way we have just seen is through Gaussian elimination with backward substitution. There is also Cramer's rule, and through computation of (A^(-1)). Another way is through decomposing the matrix A into the product of a lower triangular matrix and an upper triangular matrix. Once this has been accomplished, we can write

A x = b
L U x = b
U x = y, L y=b.

We can solve L y =b by simple backsubstitution, and then solve U x = y again by back substitution. All the difficulties of solving linear systems of equations is moved into the problem of factoring the matrix A.

LU Decomposition: Doolittle’s Method is the LU in our text

The process of LU decomposition attempts to turn a matrix into the product of a lower triangular matrix and an upper triangular matrix. The procedure is outlined here.

Brute Force Approach Leads to (n^2) Equations in (n^2)+n Unknowns

Here is a brute force attempt to decompose a matrix into upper and lower parts:

Clear[i,j]
n=3;

A=Array[a,{n,n}];
MatrixForm[A]

I hate this new command for defining a lower diagional matrix, but it works:

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]]
MatrixForm[L]

U =Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]]
MatrixForm[U]

LU=L.U
MatrixForm[%]

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2)+n unknowns. If we impose different additional n conditions to get (n^2)+n equations in (n^2)+n unknowns, we will get different methods of determining an LU decomposition. Choosing wisely will lead to simple systems that are solvable by backsubstitution, which is somewhat surprising since the system of equations is nonlinear.

Solve[eqlist,varlist]

One's on Diagonal of Lower Triangular Matrix (Doolittle's Method)

Here is Doolittle's method to LU decompose a Hilbert Matrix:

n=6;

A=HilbertMatrix[n];
MatrixForm[A]

Here we have added the additional n=6 conditions that the lower diagonal matrix has 1's on the diagonal.

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]]
Do[L[[i,i]]=1,{i,1,n}]
MatrixForm[L]

U = Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]]
MatrixForm[U]

LU=L.U;

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i+1,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns.

sol=Solve[eqlist,varlist]

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

One's on Diagonal of Upper Triangular Matrix (Crout's Method)

Here is Crout's Method to LU decompose a Hilbert Matrix:

n=6;

A=HilbertMatrix[n];
MatrixForm[A]

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]]
MatrixForm[L]

Here we have added the additional n=6conditions that the lower diagonal matrix has 1's on the diagonal.

U =Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]]
Do[U[[i,i]]=1,{i,1,n}]
MatrixForm[U]

LU=L.U;

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i+1,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns.

sol=Solve[eqlist,varlist]

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

(Subscript[l, ii])=(Subscript[u, ii]) (Choleski's Method)

Here is Choleski's method to LU decompose a Hilbert Matrix:

n=6;

A=HilbertMatrix[n];
MatrixForm[A]

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]]
MatrixForm[L]

Here we have added the additional n=6conditions that the lower diagonal elements equal the upper diagonal.

U = Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]]
Do[U[[i,i]]=L[[i,i]],{i,1,n}]
MatrixForm[U]

LU=L.U;

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i+1,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

So we have (n^2) equations in (n^2) unknowns.

sol=Solve[eqlist,varlist];

Here we see that we have successfully factored the original matrix:

MatrixForm[L/.Flatten[sol]]
MatrixForm[U/.Flatten[sol]]
MatrixForm[L.U/.Flatten[sol]]

Building a Doolittle's Module

Here we build a module to implement Doolitte’s Method to solve a system of equations.

n=4;

A=Array[a,{n,n}];
MatrixForm[A]

L =Normal[SparseArray[{i_,j_}/;i≥j->l[i,j],{n,n}]]
Do[L[[i,i]]=1,{i,1,n}]
MatrixForm[L]

Here we have added the additional n=6 conditions that the lower diagonal matrix has 1's on the diagonal.

U =Normal[SparseArray[{i_,j_}/;i≤j->u[i,j],{n,n}]]
MatrixForm[U]

LU=L.U;

eqlist={};
varlist={};
Do[eqlist=Append[eqlist,LU[[i,j]]==A[[i,j]]],{i,1,n},{j,1,n}]
Do[varlist=Append[varlist,L[[j,i]]],{i,1,n},{j,i+1,n}]
Do[varlist=Append[varlist,U[[i,j]]],{i,1,n},{j,i,n}]

TableForm[eqlist]
varlist
Length[eqlist]
Length[varlist]

This Module makes use of the Gaussian Elimination nature of the solution of the equations. The matrices (M^((k))) perform the Gaussian elimination:

(M^((k)))=(({{1, 0, 0, …, 0", 0, 0, 0}, {0, 1, 0, …, 0, 0, 0, 0}, { ⋮,   ⋮, ⋮ , ⋱ , ⋮  , ⋮ ,  ⋮,  ⋮}, {0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, -(Subscript[m, k+1,k]), 1, 0, 0}, {⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮}, {0, 0, 0, 0, -(Subscript[m, n,k]), 0, 0, 1}}))  ,  (L^((k)))=((((M^((k)))))^(-1))=(({{1, 0, 0, …, 0", 0, 0, 0}, {0, 1, 0, …, 0, 0, 0, 0}, { ⋮,   ⋮, ⋮ , ⋱ , ⋮  , ⋮ ,  ⋮,  ⋮}, {0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, (Subscript[m, k+1,k]), 1, 0, 0}, {⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮, ⋮}, {0, 0, 0, 0, (Subscript[m, n,k]), 0, 0, 1}}))

(Subscript[m, j,i])=((Subscript[((a^((i)))), j,i])/(Subscript[((a^((i)))), i,i]))

 A  = A
 ((((M^((1)))))^(-1)) (M^((1))) A  = A
 ((((M^((1)))))^(-1)) ((((M^((2)))))^(-1)) (M^((2))) (M^((1))) A  = A
 ((((M^((1)))))^(-1)) ((((M^((2)))))^(-1))  ((((M^((3)))))^(-1)) (M^((3)))(M^((2))) (M^((1))) A  = A
⋮
L U = A

L =  ((((M^((1)))))^(-1)) ((((M^((2)))))^(-1))  ((((M^((3)))))^(-1))… ((((M^((n-1)))))^(-1))= (L^((1)))(L^((2))) (L^((3))) … (L^((n-1)))
U = (M^((n-1)))… (M^((3)))(M^((2))) (M^((1))) A

MatLUDoo[A0_]:=
Module[ {A=A0,p,k,i,m,U,L,n=Length[A0]}, 
U=A0; 
L=IdentityMatrix[n]; 
Print[MatrixForm[L],MatrixForm[U]," = ",MatrixForm[A0]];
For[ p=1, p≤n-1, p=p+1, 
For[ i=p+1, i≤n, i=i+1, 
m=((A[[i,p]])/(A[[p,p]]));
L[[i,p]]=m; 
A[[i]] = A[[i]] - m A[[p]];
U = A; ]
Print[MatrixForm[L],MatrixForm[U]," = ",MatrixForm[A0]]; ];
Print[MatrixForm[L],MatrixForm[U]," = ",MatrixForm[A0]];]

You can of course eliminate the output of intermediate steps if you like.

MatLUDoo[HilbertMatrix[6]]

Compare Gaussian Elimination and LU Decomposition for Solving Linear Systems

Linear System

Here is a simple system we want to solve:

Clear[i,j,x,n,A,b,x1,x2,x3,x4]
x={x1,x2,x3,x4};
n=Length[x];
A={{10,-1,2,-4},{-1,11,-1,3},{2,-1,10,-1},{5,3,-1,8}};
b={6,25,-11,15};

Here is the system in equation form:

TableForm[Table[Sum[A[[j,i]]x[[i]],{i,1,n}]==b[[j]],{j,1,n}]]

And here it is in matrix form:

matrix1=SparseArray[{{0}},{4,4}];
Table[matrix1[[i]]=Append[A[[i]],b[[i]]],{i,1,n}];
MatrixForm[matrix1]

Gaussian Elimination

Here is my favourite Gaussian elimination with partial pivoting module:

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]==0, {Print["Gaussian Elimination Failed"], Goto[end]}]];
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Label[end];Return[A];  ] ;

Let's use it to solve the linear system of equations:

sol=Array[f,n];
Do[sol[[i]]=GaussElim3[matrix1][[i,n+1]],{i,1,n}];
sol

Let's check our solution:

Table[Sum[A[[j,i]]sol[[i]],{i,1,n}]==b[[j]],{j,1,n}]

LU Decomposition (Doolittle's)

Here is Doolitte's LU Decomposition module:

MatLUDoo[A0_]:=
Module[ {A=A0,p,k,i,m,U,L,n=Length[A0]}, 
U=A0; 
L=IdentityMatrix[n]; 
For[ p=1, p≤n-1, p=p+1, 
For[ i=p+1, i≤n, i=i+1, 
m=((A[[i,p]])/(A[[p,p]]));
L[[i,p]]=m; 
A[[i]] = A[[i]] - m A[[p]];
U = A; ]
 ];
Return[{L,U}];]

Let's use it to solve the linear system of equations. First, solve L y = b for y via backsubstitution:

y=Array[f,n];
L=MatLUDoo[A][[1]];
Do[y[[i]]=((b[[i]]-(∑_(j=1)^(i-1))y[[j]]L[[i,j]])/(L[[i,i]])),{i,1,n}];
y

Second, solve U x = y for x:

x=Array[f,n];
U=MatLUDoo[A][[2]];
Do[x[[i]]=((y[[i]]-(∑_(j=i+1)^n)x[[j]]U[[i,j]])/(U[[i,i]])),{i,n,1,-1}];
x

Let's check our solution:

Table[Sum[A[[j,i]]x[[i]],{i,1,n}]==b[[j]],{j,1,n}]

LU Module for Solving Linear System

Here is a module that takes input as the matrix A and vector b, and returns the solution to the associated system of linear equations.

SolveSystemLUDoo[A0_,b0_]:=
Module[ {A=A0,b=b0,p,k,i,m,U,L,n=Length[A0],x,y,f}, 
U=A0; 
L=IdentityMatrix[n]; 
For[ p=1, p≤n-1, p=p+1, 
For[ i=p+1, i≤n, i=i+1, 
m=((A[[i,p]])/(A[[p,p]]));
L[[i,p]]=m; 
A[[i]] = A[[i]] - m A[[p]];
U = A; ]
 ];

y=Array[f,n];
Do[y[[i]]=((b[[i]]-(∑_(j=1)^(i-1))y[[j]]L[[i,j]])/(L[[i,i]])),{i,1,n}];
x=Array[f,n];
Do[x[[i]]=((y[[i]]-(∑_(j=i+1)^n)x[[j]]U[[i,j]])/(U[[i,i]])),{i,n,1,-1}];
Return[x];]

Let's use it to solve the linear system of equations. First, solve L y = b for y via backsubstitution:

SolveSystemLUDoo[A,b]

Check our result:

Table[Sum[A[[j,i]]SolveSystemLUDoo[A,b][[i]],{i,1,n}]==b[[j]],{j,1,n}]

This module can be used to replace Mathematica's built-in Solve command if the equations are linear.

x={x1,x2,x3,x4};
eqs=Table[Sum[A[[j,i]]x[[i]],{i,1,n}]==b[[j]],{j,1,n}]
Solve[eqs,x]

which is really just using the Mathematica's built-in LinearSolve command:

LinearSolve[A,b]

Here is a random system to check our module on

SeedRandom[4];
n=27;
A=HilbertMatrix[n];
b=RandomInteger[{-10,10},n];

SolveSystemLUDoo[A,b]

LinearSolve[A,b]

SolveSystemLUDoo[A,b]-LinearSolve[A,b]

Nice work.

Another example

SeedRandom[4];
n=12;
A=HilbertMatrix[n];
b=RandomReal[{0,1},n];

SolveSystemLUDoo[A,b]

LinearSolve[A,b]

Notice the two results don't agree, since the Gaussian elimination didn't proceed very well in either case. This is an ill-conditioned linear system.

MatDet[A]

Hmmm, very small indeed. Not suprising that there were issues.

Table[Sum[A[[j,i]]SolveSystemLUDoo[A,b][[i]],{i,1,n}]==b[[j]],{j,1,n}]