Solving  Linear Systems, Existence and Uniqueness

Solving and Uniqueness

Consider the set of equations
	 x + 2y +  z  +  t  = a
	-x + 4y - 2z         = b
	 x + 3y + 4z + 5t = c
	 x          +  z  +   t = d
We illustrate how to use a matrix representation to write these out and solve them…

Start with the matrix of coefficients of the variables, mymatrix:

mymatrix  =  {
{1, 2, 1, 
            1}, 
{-1, 4, -2, 
            0}, 
{1, 3, 4, 5}, 
{1, 0, 1, 
            1}};
mymatrix // MatrixForm

The system of equations will only have a unique solution if the determinant of mymatrix is nonzero.

Det[mymatrix]

Now define vectors for (Overscript[x, ⇀])  and (Overscript[b, ⇀]) in \(A\+_\ x\&⇀\  = \ b\&⇀\)

myx  = {x, y, z, t};

myb  =  {a, b, c, d};

The left-hand side of the first equation will be

(mymatrix . myx)[[1]]

and the left-hand side of all four equations will be

lhs  =  mymatrix . myx;
lhs // MatrixForm

Now define an indexed variable linsys with four entries, each being one of the equations in the system of interest:

linsys[i_Integer]  :=  lhs[[i]]  ==  myb[[i]]

linsys[2]

Solving the set of equations for the unknowns (Overscript[x,  ->])

linsol  = 
      Solve[{linsys[1], linsys[2],  linsys[3],  linsys[4]}, myx]

Doing the same thing a different way, using Mathematica's LinearSolve function:

?LinearSolve

LinearSolve[mymatrix, myb]

And yet another way, based on \(TraditionalForm\`x\&⇀ = \(\(A\+_\^\(-1\)\) 
            A\+_\ x\&⇀\  = \ \(A\+_\^\(-1\)\) 
            b\&⇀\)\)

Inverse[mymatrix] . myb // MatrixForm

When determinants are zero

Create a matrix with one row as a linear combination of the others

myzeromatrix = 
{mymatrix[[1]], 
mymatrix[[2]], 
p*
              mymatrix[[1]]  +  q*mymatrix[[2]]  +  
            r*mymatrix[[4]], 
mymatrix[[4]]};
myzeromatrix // MatrixForm

Det[myzeromatrix]

LinearSolve[myzeromatrix, myb]

This was not expected to have a solution because one of four equations in the system was a linear combination of others in the system. Effectively, we were asking Mathematica to solve a system of three equations in four unknowns. The rank of a square matrix of coefficients is equal to the number of linearly independent equations in the system. The null space of the matrix will be empty when the equations are all linearly independent.

MatrixRank[mymatrix]
MatrixRank[myzeromatrix]

NullSpace[mymatrix]
NullSpace[myzeromatrix]

Try solving this inhomogeneous system of equations using Solve:

zerolhs =  myzeromatrix . myx

zerolinsys[i_Integer]  :=  
      zerolhs[[i]]  ==  myb[[i]]

Table[zerolinsys[i], {i, 4}] // MatrixForm

zerolinsol  = Solve[Table[zerolinsys[i], {i, 4}], myx]

No solution, as expected,  Let's see what happens if we ask Mathematica to solve the homogeneous problem:

zerolinsol  = 
      Solve[Table[
          zerolinsys[i] /. {a -> 0, b -> 0, c -> 0, 
              d -> 0}, {i, 4}], myx]

In this case, Mathematica gives a relationship between the variables, but because there are fewer equations than variables, there is still no unique solution.

Determinants

Clear[rv]

Start by building a routine to make vectors containing six random numbers on the interval {-1,1}:

rv[i_]  :=  rv[i]  = 
        Table[Random[Real, {-1, 1}], {j, 6}]

rv[1]

rv[2]

Now use rv to make a 6 x 6 matrix, then find its determinant:

randmat  =  Table[rv[i], {i, 6}]

Det[randmat]

Switching two rows changes the sign but not the magnitude of the determinant:

Det[{rv[2], rv[1], rv[3], rv[4], rv[5], rv[6]}]

Det[{rv[3], rv[1], rv[2], rv[4], rv[5], rv[6]}]

Multiply one row by a constant and calculate determinant:

Clear[a]
Det[{a*rv[2], rv[1], rv[3], rv[4], rv[5], rv[6]}]

Multiply two rows by a constant and calculate determinant:

Det[{a*rv[2], a*rv[1], rv[3], rv[4], rv[5], rv[6]}]

Multiply all rows by a constant and calculate determinant:

Det[a {rv[2], rv[1], rv[3], rv[4], rv[5], rv[6]}]

Clear[a, b, c, d, e]

lindepvec  =  
      a*rv[1]  +  b*rv[2]  +  c*rv[3]  +  d*rv[4]  +  
        e*rv[5]

Example of numerical precision: if one row of a 6 x 6 matrix is a linear combination of the other five rows, its determinant should evaluate to zero…

Det[{rv[1], rv[2], rv[3], rv[4], rv[5], lindepvec}]

Chop[Det[{rv[1], rv[2], rv[3], rv[4], rv[5], lindepvec}]]

Showing that the determinant of a product of matrices is the product of determinants

Creating a symbolic matrix

symvec  =  {a, a, a, c, c, c};

permuts  =  Permutations[symvec]
permuts // Dimensions

\(\(symbmat\  = \ {
	
          permuts[\([1]\)], 
	
          permuts[\([12]\)], 
	
          permuts[\([6]\)], 
	
          permuts[\([18]\)], 
	
          permuts[\([17]\)], 
	
          permuts[\([9]\)]
};\)\)
symbmat // MatrixForm

detsymbmat  =  Simplify[Det[symbmat]]

Creating a matrix of rational numbers

randommat  =  
        Table[Table[
            Random[Integer, {-100, 100}]/Random[Integer, {-100, 100}], {i, 6}], {j, 6}];
randommat // MatrixForm

detrandomat  =  Det[randommat]

checka  = Det[symbmat . randommat] // Simplify

detrandomat*detsymbmat   ==  checka

Does the determinant of a product depend on the order of multiplication?

checkb  = Det[randommat . symbmat] // Simplify

checka == checkb

However, the product of two matrices depends on which matrix is on the left and which is on the right

(randommat . symbmat  -  symbmat . randommat) // Simplify // 
      MatrixForm

Linear Transformations

We now demonstrate the use of matrix multiplication for manipulating an object, specifically an octohedron. The Octahedron is made up of eight polygons and the initial coordinates of the vertices were set to make a regular octahedron with its main diagonals parallel to axes x,y,z. The faces of the octahedron are colored so that rotations and other transformations can be easily tracked.

<< Graphics`Polyhedra`

Show[Polyhedron[Octahedron]]

Show[Polyhedron[Octahedron], 
      ViewPoint -> {-0.007,  -1.995,  -0.135}]

Polyhedron[Octahedron] // InputForm

The object coloct is defined below to draw an octahedron and it invokes the Polygon function to draw the triangular faces by connecting three points at specific numerical coordinates. Note that each triangles' three vertices are each specified by a vector.

Clear[coloct];
coloct = 
      Graphics3D[{
{Hue[0/8], 
            Polygon[{{0. , 0. , 1.4142135623730951}, {1.4142135623730951, 
                  0. , 0. }, {0. , 1.4142135623730951, 0. }}]}, {Hue[1/8], 
            Polygon[{{0. , 0. , 1.4142135623730951}, {0. , 
                  1.4142135623730951, 0. }, {-1.4142135623730951, 0. , 
                  0. }}]}, {Hue[2/8], 
            Polygon[{{0. , 0. , 
                  1.4142135623730951}, {-1.4142135623730951, 0. , 
                  0. }, {0. , -1.4142135623730951, 0. }}]}, {Hue[3/8], 
            Polygon[{{0. , 0. , 
                  1.4142135623730951}, {0. , -1.4142135623730951, 
                  0. }, {1.4142135623730951, 0. , 0. }}]}, {Hue[4/8], 
            Polygon[{{1.4142135623730951, 0. , 
                  0. }, {0. , -1.4142135623730951, 0. }, {0. , 
                  0. , -1.4142135623730951}}]}, {Hue[5/8], 
            Polygon[{{1.4142135623730951, 0. , 0. }, {0. , 
                  0. , -1.4142135623730951}, {0. , 1.4142135623730951, 
                  0. }}]}, {Hue[6/8], 
            Polygon[{{0. , 
                  0. , -1.4142135623730951}, {0. , -1.4142135623730951, 0. }, {-1.4142135623730951, 0. , 0. }}]}, {Hue[
              7/8], Polygon[{{0. , 1.4142135623730951, 0. }, {0. , 
                  0. , -1.4142135623730951}, {-1.4142135623730951, 
                  0. , 0. }}]}}, {AmbientLight -> GrayLevel[0], 
          AspectRatio -> 1, Axes -> False, 
          AxesEdge -> Automatic, AxesLabel -> None, 
          AxesStyle -> Automatic, Background -> Automatic, 
          Boxed -> True, BoxRatios -> Automatic, 
          BoxStyle -> Automatic, ColorOutput -> Automatic, 
          DefaultColor -> Automatic, Epilog -> {}, 
          FaceGrids -> None, ImageSize -> Automatic, 
          Lighting -> True, 
          LightSources -> {{{1. , 0. , 1. }, 
                RGBColor[1, 0, 0]}, {{1. , 1. , 1. }, 
                RGBColor[0, 1, 0]}, {{0. , 1. , 1. }, RGBColor[0, 0, 1]}}, 
          Plot3Matrix -> Automatic, PlotLabel -> None, 
          PlotRange -> All, PlotRegion -> Automatic, 
          PolygonIntersections -> True, Prolog -> {}, 
          RenderAll -> True, Shading -> True, 
          SphericalRegion -> True, Ticks -> Automatic, 
          ViewCenter -> Automatic, 
          ViewPoint -> {1.3, -2.4, 2. }, 
          ViewVertical -> {0. , 0. , 1. }, 
          DefaultFont :> $DefaultFont, 
          DisplayFunction :> $DisplayFunction, 
          FormatType :> $FormatType, 
          TextStyle :> $TextStyle}]

Show[coloct, Lighting -> False]

Now a matrix tmat is defined that will operate on the coordinate vectors fed to the Polygon function in coloct. In this first example, each z coordinate will be changed to -z:

Clear[tmat];
tmat =  {{1, 0, 0}, {0, 1, 0},  {0, 
            0, -1}};
tmat // MatrixForm

The Replace is used to modify the coordinates fed to the Polygon function in coloct by matrix multiplication of the vertex vectors:

Show[
      coloct /. {Polygon[{a_List ,  b_List ,  c_List}] -> 
            Polygon[{tmat . a, tmat . b, tmat . c}]}, 
      Lighting -> False]

The program seetrans does the same thing by feeding a matrix tranmat to operate on the vertex vectors before Polygon is executed:

seetrans[tranmat_] :=  
      Show[coloct /. {Polygon[{a_List ,  b_List ,  c_List}] -> 
              Polygon[{tranmat . a, tranmat . b, tranmat . c}]}, 
        Lighting -> False]

When tranmat is the identity matrix, the octagon is rendered in its initial orientation.

seetrans[{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}]

The next command rotates the octagon by 45° about the vertical (z) axis…

seetrans[{{Cos[Pi/4], Sin[Pi/4], 0}, {Sin[-Pi/4], Cos[Pi/4], 
          0}, {0, 0, 1}}]

And the next command both rotates it 180° about x and y and stretches it by a factor of 5 along z.

seetrans[{{-1, 0, 0}, {0, -1, 0}, {0, 0, 5}}]