Solving  Linear Systems, Existence and Uniqueness

Solving Systems of Equations

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 (OverVector[x])  and (OverVector[b]) in (UnderBar[A]) (OverVector[x]) = (OverVector[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]

Solving Matrix-Vector Equations

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

?LinearSolve

LinearSolve[mymatrix,myb]

And yet another way, based on (OverVector[x])=((UnderBar[A])^(-1))(UnderBar[A]) (OverVector[x]) = ((UnderBar[A])^(-1))(OverVector[b])

Inverse[mymatrix].myb//MatrixForm

And yet even another way, a very efficient LinearSolveFunction can be produced by LinearSolve. This function will operate on any rhs vector of the appropriate length. This would be an efficient way to find the numerical solution to a known matrix, but for many different rhs (OverVector[b]).

mymatrixsol= LinearSolve[mymatrix];

The result can be applied as a function calling a vector :

mymatrixsol[myb]
Simplify[mymatrixsol[myb]]

Singular Matrices: Zero Determinants

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

zerolinsolhet =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:

zerolinsolhom =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 and Numerical Approximations to Zero

Clear[rv]

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

rv[i_]:=rv[i]=Table[RandomReal[{-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}]

However, numerical precision does

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

Determinants and the Order of Matrix Multiplication

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 random  rational numbers

RandomMat=Table[Table[((RandomInteger[{-100,100}])/(RandomInteger[{-100,100}])),{i,6}],{j,6}];
RandomMat

DetRandomMat = Det[RandomMat]

CheckA =Det[SymbMat.RandomMat]//Simplify

DetRandomMat*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

Visualization Example: Polyhedra

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.

Needs["PolyhedronOperations`"];
Show[PolyhedronData["Octahedron"]]

Above, the color of the three dimensional object derives from the colors in the light sources. For example, note that there appears to be a blue light pointing down from the left. The lights stay fixed as we rotate the object. If Lighting -> None, then the polyhedron's faces will appear to be black.

Show[PolyhedronData["Octahedron"],Lighting->None]

We can extract data from the Octahedron, and build our own with individually colored faces. We will need the individual colors to identify what happens to the polyhedron under linear transformaions.

PolyhedronData["Octahedron","Faces"]

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 that we obtain from the Octahedron data. Because we will turn off lighting, we will ask that each of the faces glow, using the Glow graphics directive

octa = {p[1],p[2],p[3],p[4], p[5],p[6]} = PolyhedronData["Octahedron","Faces"][[1]];
colface[i_] :=Glow[Hue[i/8]] ;
ColOct=
{{colface[0],Polygon[{p[4],p[5],p[6]}]},{colface[1],Polygon[{p[4],p[6],p[2]}]},{colface[2],Polygon[{p[4],p[2],p[1]}]},{colface[3],Polygon[{p[4],p[1],p[5]}]},{colface[4],Polygon[{p[5],p[1],p[3]}]},{colface[5],Polygon[{p[5],p[3],p[6]}]},{colface[6],Polygon[{p[3],p[1],p[2]}]},{colface[7],Polygon[{p[6],p[3],p[2]}]}};

Show[Graphics3D[ColOct],Lighting->None]

Linear Transformations: Matrix Operations on Polyhedra

Here we write a function that will take a matrix and operate on all the points (i.e., position vectors) of all the triangles (The function uses a pattern for the arguments to Polygon, and then matrix multiplies the matches). The function will also take some text to help us visually identify which linear operation has been performed.

transoct[tmat_,description_String] :={ColOct/.{Polygon[{a_List,b_List,c_List}]->Polygon[{tmat.a,tmat.b,tmat.c}]},Text[Style[MatrixForm[tmat]],{0,0,-.25}],Text[Style[description,Darker[Red]],{0,0,.25},Background-> White]}

Show[Graphics3D[transoct[{{1,0,0},{0,1,0},{0,0,-1}},"mirror-[001]"]],Lighting->None]

We build a set of octahedra with different types of linear operations. We will visualize them in a few steps.

identity=(({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}));  rot90[001]=(({{Cos[π/2], Sin[-π/2], 0}, {Sin[π/2], Cos[π/2], 0}, {0, 0, 1}})); ref[010]=(({{1, 0, 0}, {0, -1, 0}, {0, 0, 1}}))
o[1,1] = transoct[identity,"original"];
o[1,2] = transoct[rot90[001],"90-[001]"];
o[1,3] = transoct[ref[010],"m-[010]"];
o[2,1] = transoct[ref[010].rot90[001],"90-[100] then m-[010]"];
o[2,2] = transoct[rot90[001].ref[010],"m-[010] then 90-[100]"];

Mathematica has many common geometric transformations available as functions: Here, we extract information about counter-clockwise rotation about [110].

RotationTransform[Pi,{1,1,0}]

o[2,3] = transoct[(({{0, 1, 0}, {1, 0, 0}, {0, 0, -1}})),"180-[110]"];

We use Manipulate and GraphicsGrid to visualize the effect of all the transformations.  We use Manipulate's sliders to set the view point for all images simultaneously. To do this, we right a function that will place the viewpoint someplace on a sphere of radius 3.

sc[θ_,ϕ_] := 3{Cos[θ] Sin[ϕ], Sin[θ] Sin[ϕ], Cos[ϕ]}
Manipulate[GraphicsGrid[Table[Show[Graphics3D[o[i,j]],Lighting->None,ViewPoint->sc[θ,ϕ],ImageSize->{200,200},PlotRange->{{-1,1},{-1,1},{-1,1}}],{j,3},{i,2}]],{{θ,2.1},0,2 π},{{ϕ,-1.4},-π/2,π/2}]

Linear Transformations: Symmetry Operations on Crystal Structures

Group theory arises from the mathematics of relating the symmetry of objects.  In crystallography, it is applied to the classification of different types of arrangments of points (i.e. crystal lattices) and the symmetry operations include multiplication of matrices representing common symmetry operations. In this example, we will construct an FCC unit cell and then show that the cell remains invariant under different symmetry operations.  The set of all operations that leave a lattice invariant uniquely defines the space group of the lattice.  We build an FCC from graphics primiitves, and then translate it so the center is at the origin.

corners = Flatten[Table[{i,j,k},{i,0,1},{j,0,1},{k,0,1}],2]
faces= Join[Permutations[{0.5,0.5,0}],Permutations[{0.5,0.5,1}]]
fccsites = Join[faces,corners]
srad = (Sqrt[2])/4;
FCC = Table[Sphere[fccsites[[i]],srad],{i,1,14}]
axes ={ {RGBColor[1,0,0,.5],Cylinder[{{0,0,0},{2,0,0}},.05]},{RGBColor[0,1,0,.5],Cylinder[{{0,0,0},{0,2,0}},.05]},{RGBColor[0,0,1,.5],Cylinder[{{0,0,0},{0,0,2}},.05]}};
fccmodel = Translate[Join[FCC,axes],{-.5,-.5,-.5}]
Graphics3D[fccmodel]

We visualize the 3-fold rotation about [111], and the roto-inversion (OverBar[3 ])(an inversion of all axis through the center, followed by a 3-fold rotation)

bbox = 1.25{{-1,1},{-1,1},{-1,1}};

gPlain = Graphics3D[fccmodel, PlotRange->bbox,SphericalRegion-> True]

fccRot= Rotate[fccmodel,2π/3,{1,1,1}];

gRot = Graphics3D[fccRot, PlotRange->bbox,SphericalRegion-> True]

fccInvRot = Rotate[GeometricTransformation[fccmodel,(({{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}}))],2π/3,{1,1,1}];

gRotInv= Graphics3D[fccInvRot, PlotRange->bbox,SphericalRegion-> True]

vp=ViewPoint/.Options[Graphics3D]
vv=ViewVertical/.Options[Graphics3D]

GraphicsRow[{Show[gRot,ViewPoint->Dynamic[vp],ViewVertical->Dynamic[vv]],Show[gRot,ViewPoint->Dynamic[vp],ViewVertical->Dynamic[vv]],Show[gRotInv,ViewPoint->Dynamic[vp],ViewVertical->Dynamic[vv]]},ImageSize-> Full]