Tutorial file 5. MA530. Linear algebra

Examples in this file correspond to O'Neil "Advanced Engineering Mathematics" 6th ed.
	Compare the analytical solutions in the book  with the computer solutions in this file. 
	In addition, please read the following file in the Mathematica help system - Advanced Documentation: Linear Algebra

Vector visualization in 2D using <<Graphics`Arrow` package

To visualize vectors on a 2D plane in Mathematica using the <<Graphics`Arrow` package we need to specify the starting point ((Subscript[x, 0]), (Subscript[y, 0])) and the ending point ((Subscript[x, 1]), (Subscript[y, 1])).  A vector will be presented as an arrow from ((Subscript[x, 0]), (Subscript[y, 0])) to  ((Subscript[x, 1]), (Subscript[y, 1])).
	First we need to download the package:

<<Graphics`Arrow`

Now let us plot the vector (-1,4):

Show[Graphics[Arrow[{0, 0}, {-1, 4}]], 
        Axes -> True];

We can plot several vectors in different colors:

Show[
        Graphics[{Hue[0], Arrow[{0, 0}, {1, 1.2}], Hue[0.6], 
            Arrow[{1, 2}, { .3,  .2}], Hue[0.8], 
            Arrow[{ .4,  .2}, {0.5, 1.4}]}], Axes -> True, 
        AspectRatio -> 1, AxesOrigin -> { .5, 1}];

Please use this example to learn options - Hue, AspectRatio, Axes, AxesOrigin.
For example, here is the same plot but the AxesOrigin is at (0,0) and AspectRatio is "Automatic":

Show[%, AspectRatio -> Automatic, 
        AxesOrigin -> {0, 0}];

Vector visualization using <<Graphics`PlotField` and <<Graphics`PlotField3D` packages

Two other Mathematica packages <<Graphics`PlotField` and <<Graphics`PlotField3D`also can be used to visualize vectors in respectively 2 and 3 dimensions.
	These packages are more advanced, and we will use them to visualize vector fields. But,  the function ListPlotVectorField3D might help to draw typical vector pictures. 
	For example let us plot the vector (-1,4,1), see also figure 5.2. in the textbook:

<< Graphics`PlotField3D`

ListPlotVectorField3D[{{{0, 0, 0}, {-1, 4, 1}}}, 
        ColorFunction -> Hue, Boxed -> False, Axes -> True, 
        BoxRatios -> {1, 4, 1}, 
        AxesEdge -> {{-1, -1}, {1, -1}, {1, -1}}, 
        AxesLabel -> {"x axis", "y axis", "z axis"}, 
        AxesStyle -> Thickness[0.003], 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 500];

Please use these plot to learn more about options of three dimensional plots in Mathematica. Try to change the following options: ImageSize,  TextStyle, AxesStyle, AxesLabel, AxesEdge, BoxRatios, Axes, Boxed, 	ColorFunction. Use the Mathematica help system.

Here is the same plot with changing the option VectorHeads->True:

ListPlotVectorField3D[{{{0, 0, 0}, {-1, 4, 1}}}, 
        VectorHeads -> True, ColorFunction -> Hue, Boxed -> False, 
        Axes -> True, BoxRatios -> {1, 4, 1}, 
        AxesEdge -> {{-1, -1}, {1, -1}, {1, -1}},  
        AxesLabel -> {"x axis", "y axis", "z axis"}, 
        AxesStyle -> Thickness[0.003], 
        TextStyle -> {FontFamily -> "Times", 
            FontSize -> 12}, ImageSize -> 500];

Here is the same plot with changing the option VectorHeads->True:

Show[%, ViewPoint -> {-2, -2, 1}, ImageSize-> 800];

Show[%, ViewPoint -> {0, -2, 0},ImageSize-> 500];

These packages are to visualize vector fields, this notion is introduced in section 12.3 of the textbook. Here, as an example, I reproduce the figure 12.19. :

PlotVectorField3D[{Cos[x], Exp[-x]*Sin[y], z - y}, {x, -1, 
          1}, {y, -1, 1}, {z, -1, 1}, PlotPoints -> 6, 
        VectorHeads -> True, Axes -> True];

Algebra of vectors

In Mathematica vectors are represented as lists of objects, for example {a,b,c,d} is a vector.
	Vector addition and multiplication of a vector be a scalar as well as the other operations can be done both symbolically and numerically:

A = {a, b, c}

B = {d, f, g}

S = A + B

T = p*A

We can calculate the norm of a vector

Nr = Norm[A]

Norm[B]

Norm[T]

and dot and cross products:

DP = A . B

CP = Cross[A, B]

Now if define numerically the  components of A and B, we immediately get all previous results numerically:

a = 2; b = 3; c = -2; d = 4; f = -3; g = 1; p = 3;

A

B

S

T

Nr

DP

CP

Clear[A, B, S, T, Nr, DP, CP]

Solution of homogeneous systems of linear equations

Example 7.24 p. 273

Systems of linear equations can be solved using the functions Solve (searches for an analytical solution) and NSolve (searches for a numerical solution), for example let us use these functions to solve the equation of Example 7.24 p. 273:

S = {x1 - 3*x2 + x3 - 7*x4 + 4*x5 == 0, 
          x1 + 2*x2 - 3*x3 == 0, 
          x2 - 4*x3 + x5 == 0};

Solve[S, {x1,x2,x3,x4,x5}]

NSolve[S, {x1,x2,x3,x4,x5}]

The operator LinearEquationsToMatrices in the package  LinearAlgebra`MatrixManipulation` allow us to present this system in a matrix form:

<< LinearAlgebra`MatrixManipulation`

M = LinearEquationsToMatrices[S, {x1, x2, x3, x4, x5}]

The list M consists of two components, the first component is the matrix of the coefficients of the system S:

MatrixForm[M[[1]]]

And the second component is a vector of the right-hand sides of the equations. In this example it is the zero vector since we consider a homogeneous system:

MatrixForm[M[[2]]]

Application of the function LinearSolve to a homogeneous system is not productive, we have got just the trivial solution:

LinearSolve[M[[1]], M[[2]]]

We can make a matrix of the system by adding the right-hand side vector to the matrix of the coefficients, and apply Gaussian elimination using function RowReduce:

SM = AppendRows[M[[1]], Transpose[{M[[2]]}]]

MatrixForm[SM]

RowReduce[SM]

Another way to solve this system is to get explicitly the null space of the matrix of the coefficients:

NS = NullSpace[M[[1]]]

Now any linear combination of the two vectors spanning the null space of the matrix of coefficients is a solution of the given equation:

Simplify[
      M[[1]] . (c1 *NS[[1]] + c2*NS[[2]])]

Please compare the three considered solutions of this system; see that the answer that we get using functions Solve, RowReduce, and NullSpace are the same.

Clear[M, SM, NS, S]

Example 7.30 p. 282

We consider the following system of linear equations, see Example 7.30 p. 282 :

S = {-4*x1 + x2 - 7*x3 == 0, 
          2*x1 + 9*x2 - 13*x3 == 0, 
          x1 + x2 + 10*x3 == 0};

Firstly let us apply functions Solve and NSolve:

Solve[S, {x1,x2,x3}]

NSolve[S, {x1,x2,x3}]

These methods evident that this system has the only solution x1=0, x2=0, x3=0

Let us also consider this system in a matrix form:

<< LinearAlgebra`MatrixManipulation`

M = LinearEquationsToMatrices[S, {x1, x2, x3}]

The list M consists of two components, the first component is the matrix of the coefficients of the system S:

MatrixForm[M[[1]]]

And the second component is a vector of the right-hand sides of the equations. In this example it is the zero vector since we consider a homogeneous system:

MatrixForm[M[[2]]]

We use Gaussian elimination to solve this system:

SM = AppendRows[M[[1]], Transpose[{M[[2]]}]]

MatrixForm[SM]

RowReduce[SM]

Therefore the matrix of coefficients reduces to the identity matrix, providing that only the trivial solution exists.

The third way to see that this system has only trivial solution is to calculate the null space or the determinant of the matrix of the coefficients:

NS = NullSpace[M[[1]]]

Det[M[[1]]]

Clear[M, NS, S]

Solution of nonhomogeneous systems of linear equations

Example 7.33 p. 286

Consider the following system of linear equations  (Example 7.33 p. 286):

S = {-3*x1 + 2*x2 + 2*x3 == 8, 
          x1 + 4*x2 - 6*x3 == 
            1, -2*x2 + 2*x3 == -2};

Solve[S, {x1,x2,x3}]

NSolve[S, {x1,x2,x3}]

Using LinearEquationsToMatrices we can present this system in a matrix form:

<< LinearAlgebra`MatrixManipulation`

M = LinearEquationsToMatrices[S, {x1, x2, x3}]

The list M consists of two components, the first component is the matrix of the coefficients of the system S:

MatrixForm[M[[1]]]

And the second component is a vector of the right-hand sides of the equations:

MatrixForm[M[[2]]]

Now we can solve this systems using LinearSolve:

LinearSolve[M[[1]], M[[2]]]

Another way to solve this system is to apply the function RowReduce:

SM = AppendRows[M[[1]], Transpose[{M[[2]]}]]

MatrixForm[SM]

RowReduce[SM]

Please compare the considered different approaches to solve this system with the analytical solution in the textbook to see that the answer is the same in all cases.

Clear[M, SM, NS, S]

Example 7.36 p. 290

Consider the following system of linear equations  (Example 7.36 p. 290):

S = {x1 - x3 + 2*x4 + x5 + 6*x6 == -3, 
          x2 + x3 + 3*x4 + 2*x5 + 4*x6 == 1, 
          x1 - 4*x2 + 3*x3 + x4 + 2*x6 == 0};

Solve[S, {x1,x2,x3,x4,x5,x6}]

NSolve[S, {x1,x2,x3}]

Using LinearEquationsToMatrices we can present this system in a matrix form:

<< LinearAlgebra`MatrixManipulation`

M = LinearEquationsToMatrices[S, {x1, x2, x3}]

The list M consists of two components, the first component is the matrix of the coefficients of the system S:

MatrixForm[M[[1]]]

And the second component is a vector of the right-hand sides of the equations:

MatrixForm[M[[2]]]

Now we can solve this systems using LinearSolve:

LinearSolve[M[[1]], M[[2]]]

Another way to solve this system is to apply the function RowReduce:

SM = AppendRows[M[[1]], Transpose[{M[[2]]}]]

MatrixForm[SM]

RowReduce[SM]

Please compare the considered different approaches to solve this system with the analytical solution in the textbook to see that the answer is the same in all cases.

Clear[M, SM, NS, S]

Using the matrix inverse to solve AX=b, Example 7.40 p. 297

Consider the following system of linear equations  (Example 7.36 p. 290):

S = {2*x1 - x2 + 3*x3 == 4, x1 + 9*x2 - 2*x3 == -8, 
          4*x1 - 8*x2 + 11*x3 == 15};

Solve[S, {x1,x2,x3}]

NSolve[S, {x1,x2,x3}]

Using LinearEquationsToMatrices we can present this system in a matrix form:

<< LinearAlgebra`MatrixManipulation`

M = LinearEquationsToMatrices[S, {x1, x2, x3}]

The list M consists of two components, the first component is the matrix of the coefficients of the system S:

MatrixForm[M[[1]]]

And the second component is a vector of the right-hand sides of the equations:

MatrixForm[M[[2]]]

Now we can solve this systems using LinearSolve:

LinearSolve[M[[1]], M[[2]]]

Another way to solve this system is to apply the function RowReduce:

SM = AppendRows[M[[1]], Transpose[{M[[2]]}]]

MatrixForm[SM]

RowReduce[SM]

One more way to solve this system is to use the matrix inverse since the matrix of the coefficients is nonsingular, the solution is X=A^(-1)B.
First we find the inverse of the matrix of the coefficients:

Ainv = Inverse[M[[1]]]

Then we find the solution of the system:

Ainv . M[[2]]

Please compare the considered different approaches to solve this system with the analytical solution in the textbook to see that the answer is the same in all cases.

Clear[M, SM, NS, S]

Determinants

Definition and properties

Let us consider an arbitrary 3x3 matrix:

A = {{a11, a12, a13}, {a21, a22, a23}, {a31, a32, a33}};

MatrixForm[A]

The determinant of A is:

Det[A]

Please compare this result with the analytical formula on page 302 of the textbook. See that these formulas are identical.

Now let us consider some determinant properties:

Theorem 8.1. p. 303: If A has a zero row, then detA=0:

Det[{{0, 0, 0}, {a21, a22, a23}, {a31, a32, a33}}]

Theorem 8.2. p. 303: Let B be formed from A by multiplying row k be a scalar C. Then det B= C det A :

DetermB = 
      Det[{{C*a11, C*a12, C*a13}, {a21, a22, a23}, {a31, a32, 
            a33}}]

Simplify[DetermB/Det[A]]

Theorem 8.3. p. 304: Let B be formed from A by interchanging two rows. Then det B= - det A :

DetermB = 
      Det[{{a21, a22, a23}, {a11, a12, a13}, {a31, a32, a33}}]

Simplify[DetermB/Det[A]]

Corollary 8.1 p. 304: If two rows of A are the same, then detA=0 :

DetermB = 
      Det[{{a11, a12, a13}, {a11, a12, a13}, {a31, a32, a33}}]

Corollary 8.2 p. 304: If for some scalar C, row k of A is C times row i, then detA=0 :

DetermB = 
      Det[{{a11, a12, a13}, C*{a11, a12, a13}, {a31, a32, a33}}]

Theorem 8.4. p. 305: Let A and B be nxn matrices. Then, det AB=  det A det B :

B = {{b11, b12, a13}, {b21, b22, b23}, {b31, b32, b33}};

Simplify[Det[A . B] - Det[A] Det[B]]

Theorem 8.6. p. 306: det A=  det (A^t) :

Det[A] == Det[Transpose[A]]

Clear[A, B, DetermB]

Example 8.1. p. 309

Evaluate det(A) for the following matrix, please compare this with an analytical solution in the book.

A = {{-6, 0, 1, 3, 2}, {-1, 5, 0, 1, 7}, {8, 3, 2, 1, 7}, {0, 
          1, 5, -3, 2}, {1, 15, -3, 9, 4}}

MatrixForm[A]

Det[A]

Clear[A]

Examples 8.2. and 8.3 p. 312

Evaluate det(A) for the following matrix. Please compare this with analytical solutions in the book.

A = {{-6, 3, 7}, {12, -5, -9}, {2, 4, -6}}

MatrixForm[A]

Det[A]

Clear[A]

Example 8.4. p. 315

Evaluate det(A) for the following matrix. Please compare this with an analytical solution in the book.

A = {{15, -7, 4, 7, 3}, {0, 12, -6, 3, 9}, {0, 0, Sqrt[2], 
          15, -4}, {0, 0, 0, Pi, 22}, {0, 0, 0, 0, e}}

MatrixForm[A]

Det[A]

Clear[A]

Eigenvalues and eigenvectors

Example 9.3. p. 326

Find eigenvalues and eigenvectors of A, please compare this with an analytical solution in the book.

A = {{1, -1, 0}, {0, 1, 1}, {0, 0, -1}}

MatrixForm[A]

One way to find eigenvalues is to calculate the characteristic polynomial first.

P[x_] = CharacteristicPolynomial[A, x]

Simplify[P[x]]

and then to solve an equation P(x)=0:

Solve[P[x] == 0]

However, the function Eigenvalues gives the result in one step:

Eigenvalues[A]

The function Eigenvectors gives the eigenvectors, please note that the eigenvectors for the given A are linearly dependent, and compare this result with the textbook to see how this result can emerge analytically:

Eigenvectors[A]

The function Eigensystem gives both the eigenvalues and the eigenvectors:

Eigensystem[A]

ClearAll[A, P]

Example 9.4. p. 327

Find eigenvalues and eigenvectors of A, please compare this with the analytical solution in the book.

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

MatrixForm[A]

Eigenvalues[A]

P = Eigenvectors[A]

Note that in this case we have two linearly independent eigenvectors:

Det[Eigenvectors[A]]

Eigensystem[A]

ClearAll[A, P]

Diagonalization

Here I list basic facts concerning diagonalization of matrices, and please also read carefully paragraph 9.2 in the textbook.

Definition 9.3. A diagonal matrix is a square matrix having all off-diagonal elements equal to zero.

Definition 9.4. An nxn matrix A is diagonalizable if there exist an nxn matrix P such that (P^(-1)) AP is a diagonal matrix.  When such P exists, we say that P diagonalized A.

Theorem 9.4. Let A be an nxn matrix. Then A is diagonalizable if it has n linearly independent eigenvectors. Further, if P is the nxn matrix having these eigenvectors as columns, then (P^(-1)) AP is the diagonal matrix having the corresponding eigenvalues down its main diagonal.

Example 9.6. p. 332

Let us diagonalize the following matrix:

A = {{-1, 4}, {0, 3}}

MatrixForm[A]

We calculate the eigenvectors

P = Eigenvectors[A]

P // MatrixForm

and check that the eigenvectors are linearly independent:

Det[P]

Now let us calculate (P^(-1)) AP:

MatrixForm[Transpose[Inverse[P]] . A . Transpose[P]]

Actually using the more general function which calculate the Jordan decomposition of a square matrix, we can get the diagonal matrix similar to the given matrix if it exists:

{m, p} = JordanDecomposition[A]

MatrixForm[m]

MatrixForm[p]

ClearAll[A, P, m, p]

Example 9.7. p. 333

Let us diagonalize the following matrix:

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

MatrixForm[A]

We calculate the eigenvectors

P = Eigenvectors[A]

P // MatrixForm

and check that the eigenvectors are linearly independent:

Det[P]

Now let us calculate (P^(-1)) AP:

MatrixForm[
      Simplify[Transpose[Inverse[P]] . A . Transpose[P]]]

Actually using another, more general function which calculates the Jordan decomposition of a square matrix, we can get the diagonal matrix similar to the given matrix if it exists:

{m, p} = JordanDecomposition[A]

MatrixForm[m]

MatrixForm[p]

ClearAll[A, P, m, p]

Example 9.8. p. 334

Let us diagonalize the following matrix:

A = {{-1, -4}, {3, -2}}

MatrixForm[A]

We calculate the eigenvectors

P = Eigenvectors[A]

Simplify[Transpose[P]] // MatrixForm

and check that the eigenvectors are linearly independent:

Det[Transpose[P]]

Now let us calculate (P^(-1)) AP:

MatrixForm[
      Simplify[Transpose[Inverse[P]] . A . Transpose[P]]]

The same result using the Jordan decomposition:

{m, p} = JordanDecomposition[A]

MatrixForm[m]

MatrixForm[p]

ClearAll[A, P, m, p]

Example 9.9. p. 334

Let us diagonalize the following matrix:

A = {{5, -4, 4}, {12, -11, 12}, {4, -4, 5}}

MatrixForm[A]

We calculate the eigenvectors

P = Eigenvectors[A]

P // MatrixForm

and check that the eigenvectors are linearly independent:

Det[P]

Now let us calculate (P^(-1)) AP:

MatrixForm[
      Simplify[Transpose[Inverse[P]] . A . Transpose[P]]]

Similar result we can get using the Jordan decomposition:

{m, p} = JordanDecomposition[A]

MatrixForm[m]

MatrixForm[p]

ClearAll[A, P, m, p]

Orthogonal and symmetric matrices

In addition to these examples please also read paragraph 9.3 in the textbook.

Definition 9.3. A square matrix A is orthogonal if and only if ((AA)^t)=(A^t) A=(Subscript[I, n])  

Therefore an orthogonal matrix is nonsingular and we can get its inverse by taking its transpose

Example 9.12. p. 340

Let us consider the following matrix:

A = {{0, 1/Sqrt[5], 2/Sqrt[5]}, {1, 0, 0}, {0, 
          2/Sqrt[5], -1/Sqrt[5]}}

MatrixForm[A]

Transpose[A] // MatrixForm

A . Transpose[A] // MatrixForm

Therefore A is an orthogonal matrix.

ClearAll[A]

Definition 9.6. A square matrix A is symmetric if A=(A^t)  

Theorem 9.10 The eigenvalues of a real, symmetric matrix are real numbers

Theorem 9.11 Let A be a real symmetric matrix. Then eigenvectors associated with distinct eigenvalues are orthogonal.

Example 9.13. p. 345

Let us consider the following matrix:

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

MatrixForm[A]

Eigenvalues[A]

P = Eigenvectors[A]

We can check that the eigenvectors are orthogonal:

P[[1]] . P[[2]]

P[[2]] . P[[3]]

P[[1]] . P[[3]]

Now let us normalize the eigenvectors, using the package <<LinearAlgebra`Orthogonalization` :

<<LinearAlgebra`Orthogonalization`

PN = Map[Normalize, P, {1}]

MatrixForm[PN]

The normalized matrix is an orthogonal matrix

Transpose[PN] . PN // MatrixForm

which diagonalizes the given matrix A:

Transpose[Inverse[PN]] . A . Transpose[PN] // MatrixForm

Clear[A, P, PN]