Linear Algebra with Mathematica

This notebook demonstrates how carry out various linear algebra computations using Mathematica.  To evaluate a cell of Mathematica commands, select it and then Shift-Return or Enter.  Open or close sections of the notebook by double-clicking.

Solving systems of linear equations

To use Mathematica effectively you need to work with lists.   Lists are entered using curly brackets with entries separted by commas.  You can extract the elements in the list using a double square bracket notation [[ ]].  For example here is a list of three numbers and the notation for extracting the third number.

l = {4,-1,2}
l[[3]]

A system of equations is just a list of equations.  Equations are written using a double equal sign == .
Here is a systems of two linear equations in two unknowns x, y.  You can give the system a name so you don’t have to type it again later.

eqs = {x + 2*y ==3 , 2*x - y == 0}

There is a command for solving systems of equations.  Here we solve the system above for the two variables x, y.  The answer is in the form of a list of substitution rules.  There is only one solution so there is only one solution rule which can be extracted as follows.

solutionrules = Solve[eqs,{x,y}]
solutionrule = solutionrules[[1]]

Finally you can get the actual list of numbers by applying the substitution rule to the vector of variables {x,y}.  You apply a subtitution rule using the /. notation

solution = {x,y}/.solutionrule

A more efficient way work with linear equations is to enter a coefficient matrix A and the vector of constants b.  A vector is just a list of numbers and a matrix is a list of the vectors which form the rows of the matrix.  The system above would be represented by

A = {{1,2},{2,-1}}
b= {3,0}

Matrices can be printed out nicely using the command MatrixForm.  Here are two ways to do it.  The // notation applies a function to whatever comes before it.

MatrixForm[A]
A//MatrixForm

Now the system of equations can be obtained as A.{x,y} ==b (using multiplication of A times the vector of variables {x,y}).

eqs = A.{x,y}==b
Solve[eqs,{x,y}]

Another alternative is to use the LinearSolve command which takes a matrix A and a vector b as arguments.

LinearSolve[A,b]

Finally, you could work with the augmented matrix [ A | b ].  The first command defines the matrix without printing it (because of the semicolon).  The second command prints it nicely.

augmented = {{1,2,3},{2,-1,0}};
augmented//MatrixForm

Now the command RowReduce gives the (unique) reduced row echelon form.  Using this, the solution can just be read off.

rref = RowReduce[augmented]
rref//MatrixForm

Example --Approximate solution of a differential equations using linear algebra

Suppose we want to find a function y(t) solving the differential equation y’’ + y = t for 0 <= t<= 1 and also satisfying the boundary conditions  y(0)= y(1) = 0.   This is an easy second order linear differential equation whose exact solution is as follows:

y[t_] := t  - Sin[t]/Sin[1]
exactplot =Plot[y[t],{t,0,1}]

We will use it as an example of how to find approximate solutions with linear algebra.  The method is also applicable to harder problems where there is no nice formula for the exact solution.  Divide the interval [0,1] at the points {t0, t1, t2, t3, t4} = {0, 1/4, 1/2, 3/4, 1}.  The values of y at these points will be called {y0, y1, y2, y3, y4}.   We have the boundary values y0 = 1 and y4 = 0 and we will find the approximate values of y1, y2, y3 by solving 3 linear equations.  Here are the equations:

{t1,t2,t3}= {1/4,1/2,3/4};
eqs = { (y2+y0-2*y1)/h^2 + y1 == t1,(y3+y1-2*y2)/h^2 + y2 ==t2,(y4+y2-2*y3)/h^2 + y3 == t3}

Insert the boundary values and the value of h and then solve using Mathematica.

eqs1 = eqs/.{h->1/4,y0->0,y4->0}
Solve[eqs1,{y1,y2,y3}]
N[%]

The exact values are

{y[.25],y[.5],y[.75]}

Example -- equilibrium of a Leontieff input-output model from economics

Suppose we have three economic sectors each producing various commodities.  Each sector need the other commodities to produces its own output (internal demands).  In addition there are demands for the commodites from consumers.  A matrix called A encodes the internal demands per unit of production.  Here is an example with three commodities.

A= {{0.2,0.3,0.1},{0.1,0.3,0.2},{0.4,0.2,0.1}};
A//MatrixForm

The first row means that to produce one unit of the first commodity requires 0.2 units of itself, 0.3 units of commodity 2 and 0.1 units of commodity 3.  If the outputs of the sectors are {x1,x2,x3} then the total internal demand is A.{x1,x2,x3}

internaldemand = A.{x1,x2,x3}

Suppose the external demands are for 2 units of commodity 1, 1 unit of commodity 2 and 3 units of commodity 3.  Then there is an external demand vector {2,1,3}.  The total demand is the sum of internal and external demands:

externaldemand = {2,1,3}
totaldemand = internaldemand + externaldemand

Next we look for outputs levels which exactly equal the total demands for each commodity.  In other words
{x1,x2,x3} = totaldemand.  This is a system of linear equations with a unique solution.

eqs = {x1,x2,x3}==totaldemand
Solve[eqs,{x1,x2,x3}]

Using matrix times vector multiplication, the Leontieff equilibrium equations can be written X = A X +d or
(I - A)X = d where X={x1,x2,x3} and I is the 3x3 identity matrix.  Setting  B = I - A gives a linear system A X = d,
where d is the external demand.

A= IdentityMatrix[3]-B;
A//MatrixForm
d=externaldemand
LinearSolve[A,d]

Matrix Equivalence, Elementary Matrices and the LU decomposition.

Here is an example to show how to do row and column operations using elementary matrices.  Here is the initial matrix A.

A={{1,2,0,1,1},{3,6,1,4,0},{2,4,1,3,-1}};
A//MatrixForm

Row reduction via elementary matrices

We use elementary matrices to clear out the first row and column.  These are obtained from the 3 x 3 identity matrix by doing the corresponding operations.  The product T2.T1 could also be called an elementary matrix whose effect will be to clear out the first column.

T1={{1,0,0},{-3,1,0},{0,0,1}};
T2={{1,0,0},{0,1,0},{-2,0,1}};
{T1//MatrixForm,T2//MatrixForm,T2.T1//MatrixForm}
T2.T1.A//MatrixForm

Next we clear out the third column, using another elementary matrix to obtain EA, the reduced row echelon form. The product P = T3.T2.T1 is a nonsingular matrix which gives the whole reduction.

T3={{1,0,0},{0,1,0},{0,-1,1}};
P=T3.T2.T1;
EA=T3.T2.T1.A;
{T3//MatrixForm,P//MatrixForm,EA//MatrixForm}

LU decomposition

For later use we compute the matrix L = Inverse[P] = Inverse[T1]Inverse[T2]Inverse[T3] (opposite order).  Note that it just records minus the multipliers used in the Gauss elimination.  L is lower triangular and U= EA is upper triangular. The original matrix is  A = L.U

L = Inverse[T1].Inverse[T2].Inverse[T3];
U= EA;
{L//MatrixForm,U//MatrixForm,L.U//MatrixForm}

Column operations and rank normal form

We can also do column operations to reduce the matrix to rank normal form.  Starting from EA we multiply on the left by 5 x 5 elementary matrices to clear out the pivot rows and then permute the columns to get the normal form.

T4={{1,-2,0,-1,-1},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}};
T5={{1,0,0,0,0},{0,1,0,0,0},{0,0,1,-1,3},{0,0,0,1,0},{0,0,0,0,1}};
T6={{1,-0,0,0,0},{0,0,1,0,0},{0,1,0,0,0},{0,0,0,1,0},{0,0,0,0,1}};
{T4//MatrixForm,T5//MatrixForm,T6//MatrixForm}
Nr=EA.T4.T5.T6;
Nr//MatrixForm

The matrices P and Q constructed above can be used to reduce A to rank normal form Nr -- a matrix with an r x r identity block and the rest zeros.  r is the rank of A.  The normal form is Nr = PAQ.

Q= T4.T5.T6;
{P//MatrixForm,Q//MatrixForm,P.A.Q//MatrixForm}