Matrices and Vectors

Matrices can be declared in Mathematica as a list of lists of the form A={{a,b},{c,d}}, where {a, b} is the first row of the matrix and {c,d} is the second row of the matrix. The Mathematica front end also provides an Insert->Table/Matrix submenu for creating and editing arrays with any specified number of rows and columns. Once you have such an array, you can edit it to fill in whatever elements you want.

A={{1,0,0},{1,1,1}}
A=(({{1, 0, 0}, {1, 1, 1}}))
v=(({{1}, {2}, {3}}))

If you compile these examples, we can see that Mathematica treats a matrix like a list of lists. To see these matrices in standard matrix form, use the command MatrixForm[ ]. A vector in Mathematica is a list. Column vectors and row vectors are represented identically.

MatrixForm[IdentityMatrix[3]]
w={1,2,3}
MatrixForm[w]

The inverse of a matrix

Inverse[A] computes the inverse of an nxn square matrix A, if the matrix is invertible, and displays an error message otherwise. A matrix whose determinant is different from 0 is invertible. Its inverse is the unique nxn square matrix B wih the property that AB = BA = IdentityMatrix[n].

M={{1,2},{3,4}};
Inverse[M]
Print[MatrixForm[M]MatrixForm[Inverse[M]],"=",MatrixForm[Inverse[M]]MatrixForm[M],"=",MatrixForm[IdentityMatrix[2]]]

Matrix and Vector Operations

We can perform matrix addition A+B and scalar multiplication 3A in the usual way. However, to do matrix multiplication we cannot use the symbol * , instead we have to use a dot!

A=(({{1, 0, 0}, {1, 1, 1}}));
w={1,2,3};
MatrixForm[3A]
MatrixForm[A+A]
v=A.w
MatrixForm[v]

For square matrices, we can compute the power of a matrix (A^n) using the command MatrixPower[matrix, power]

B={{1,1},{1,1}}; MatrixPower[B,2]

To access the  element on row i and column j of matrix A, we use A[[i, j]] or A[[ i ]] [[ j ]].

B[[2,1]]=10;
MatrixForm[B]

Linear Algebra - Eigenvectors and Eigenvalues

The function Eigenvalues[A] returns a list with the eigenvalues of the square matrix A. The function Eigensystem[A] finds both the eigenvalues and a complete linearly independent set of eigenvectors for each eigenvalue.The output is given in the form {list of eigenvalues, list of eigenvectors}. If the geometric multiplicity of an eigenvalue is less than its algebraic multiplicity, then the list of eigenvectors is completed with 0 vectors.

Consider an nxn matrix A. The Jordan Decomposition Theorem says there exists an invertible matrix S such that (S^(-1))A S=J , where J is called the Jordan canonical form of A (J is an upper triangular matrix whose diagonal entries are the eigenvalues of A, and immediately above the diagonal some entries are equal to 1). The function JordanDecomposition[A] returns a list of the form {S, J}, where J is the Jordan canonical form of A. S is one choice of matrix satisfying (S^(-1))A S=J, although it’s not the usually the one that we would get when doing the computations by hand.

The Matrix Exponential

Given any matrix A, we can compute the matrix exponential (e^At)by using the Mathematica command MatrixExp[A*t]. If (S^(-1))A S=J, where J is the Jordan canonical form of A, then (e^At)= (Se^Jt)(S^(-1)). We can test this relation for our previous example.

M = (({{2, 0, 1}, {0, 1, -1}, {-1, 0, 0}}));
Clear[t]
Print["This is the matrix exponential of M"]
MatrixExp[M*t]

Print["For a better display, use MatrixForm"]
MatrixForm[MatrixExp[M*t]]

Print["We will check that \!\(TraditionalForm\`\*SuperscriptBox[\(e\), \(Mt\)] = \\ \*SuperscriptBox[\(Se\), \(Jt\)] \*SuperscriptBox[\(S\), \(-1\)]\)"]
{S,J}=JordanDecomposition[M];

Print["The matrix exponential of the Jordan canonical form J is "]
MatrixForm[MatrixExp[J*t]]

Print["The result of the matrix Multiplication (Se^Jt) (S^(-1)) is"]
MatrixForm[S.MatrixExp[J*t].Inverse[S]]

We can now solve the system of differential equations x’ = Mx using matrix exponential (e^At), or using a fundamental matrix given by  (Se^Jt). The general solution can be expressed as x(t)=(e^At)c, but also as x(t)=(Se^Jt)d, where the coefficients c and d are related by the equation d=(S^(-1))c.

Clear[c,d,t];

Print["Two different ways of writing the general solution of the system x'=Mx are"]
c={c1,c2,c3};
sol=MatrixExp[M*t,c];
MatrixForm[sol]

d={d1, d2, d3};
sol=S.MatrixExp[J*t].d;
MatrixForm[sol]

Nonhomogeneous Systems of first order differential equations in Mathematica

We can use the method of Variation of Parameters to solve a nonhomogeneous system of the form x’=Ax+f(t). A particular solution is given by x(t)=(e^At)∫(e^(-At))f(t)dt=ϕ(t)∫((ϕ(t))^(-1))f(t)dt, where ϕ(t) is any fundamental matrix of the homogeneous system x’=Ax, and the integral denotes an antiderivative. This computation can be automated in Mathematica using the command Integrate[f], which returns an antiderivative of f.

Consider the system x' = A x + f (t), where A and f(t) are given below:

A = (({{2, 1, 0}, {0, 2, 0}, {0, 0, 5}}));
f[z_]:={Sin[z],3,2z}
sol=MatrixExp[A*t].(Integrate[MatrixExp[-A*t].f[t],t]);
MatrixForm[sol]

Systems of homegeneous first order differential equations in Mathematica

There are several ways to Solve linear homogeneous systems in Mathematica, and we have already seen one above, using matrix exponentials and fundamental matrices. Another way is to use the function DSolve[..]

DSolve

DSolve[eqn,y,x]
solves a differential equation for the function y, with independent variable x.
DSolve[{eqn1,eqn2,…},{y1,y2,…},…]
solves a list of differential equations for y1, y2, ... .

DSolve returns results as lists of rules. This makes it possible to return multiple solutions to an equation. For a system of equations, possibly multiple solution sets are grouped together inside curly brackets { }. You can use the rules to substitute the solutions into other calculations.

Clear[x,y,z,t]
Print["The general solution of the system x'=x+y, y'=4x+y is"]
sol=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t]},{x[t],y[t]},t]

Print["The unique solution of the initial value problem x'=x+y, y'=4x+y, x(0)=2, y(0)=0 is"]
solivp=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t],x[0]==2,y[0]==0},{x[t],y[t]},t]

The solution is given as a transformation rule x[t] -> function1, y[t] -> function2. To use x[t] and y[t] in later computations, for instance to plot the solution of the initial value problem, we use the replacement operator /.

ParametricPlot[{x[t],y[t]} /. solivp,{t,-2,1}]

Generalized Eigenvectors

We can also use Mathematica to solve the system the classical way, by Analyzing the eigenvalues of the coefficient matrix, finding a complete set of linearly independent eigenvectors and generalized eigenvectors, and writing the general solution of the system as a combination of linearly independent solutions. This method is outlined below.

Let’s consider the system x’=Mx whose coefficient matrix of the system is exactly the matrix M from the section above. The matrix has a repeated eigenvalue 1, with algebraic multiplicity 3, and geometric multiplicity 1. We have one linearly independent eigenvector v1={0,1,0}, as we discover by using the command Eigensystem[..]

M = (({{2, 0, 1}, {0, 1, -1}, {-1, 0, 0}}));
e=Eigensystem[M]
v1= e[[2]][[1]] (* the only eigenvector *)

MatrixForm[JordanDecomposition[M][[2]]]

The matrix has a repeated eigenvalue 1, with algebraic multiplicity 3, and geometric multiplicity 1. We have one linearly independent eigenvector v1={0,1,0}. We need two more generalized eigenvectors. We can use the command Solve[equation, variables] to solve the two equations (M-λ(Subscript[I, 3]))(Subscript[v, 2])=(Subscript[v, 1])and  (M-λ(Subscript[I, 3]))(Subscript[v, 3])=(Subscript[v, 2]).

Clear[v2,x]
v2=Table[x[j],{j,3}];
Solve[(M-IdentityMatrix[3]).v2==v1,v2]

Therefore v2 has the form v2 ={x[1],x[2],x[3]}={1, x[2], -1} = x[2]*v1 + {1, 0, -1}. We can pick v2 to be {1, 0, -1}.

v2={1,0,-1};
Clear[v3,x]
v3=Table[x[j],{j,3}];
Solve[(M-IdentityMatrix[3]).v3==v2,v3]

Therefore v3 has the form v3 = {1, x[2], 0} = x[2]*v1 + {1, 0, 0}. We can pick v3 to be {1, 0, 0}.

A complete set of linearly independent eigenvectors and generalized eigenvectors is :

v1={0,1,0}
v2={1,0,-1}
v3={1,0,0}

StreamPlot and Phase Portraits

As before, we can visualize the phase portraits of systems of differential equations using the command StreamPlot[..]. To plot the phase portrait of the system x’=x+y, y’=y, we write:

StreamPlot[{x+y,y},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{1,1},Green},{{-3,2},Orange},Automatic}}]

Another system: x’=y, y’=-x

StreamPlot[{y,-x},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{0,1},Red},{{1,1},Green},Automatic}}]

Interactive plot of the vector field of various systems of linear differential equations, using StreamPlot[..] and Manipulate[..]

This plot shows the various types of stability of the unique equilibrium point of a 2x2 system x’=mx, defined by a 2x2 matrix m.

Manipulate[Row[{Column[{Text["MATRIX"],MatrixForm[m],Text["EIGENVALUES"], N[Eigenvalues[m]]}, Spacings->2], StreamPlot[m.{x,y},{x,-1,1},{y,-1,1},StreamScale->Large,StreamColorFunction->"Rainbow",ImageSize->Large]}
],
{{m,({{1,0},{0,2}})},{({{1,0},{0,2}})->"Nodal source",({{1,1},{0,1}})->"Degenerate source",({{0,1},{-1,1}})->"Spiral source",({{-1,0},{0,-1}})->"Nodal sink",({{-1,0},{0,-2}})->"Nodal sink",({{-1,1},{0,-1}})->"Degenerate sink",({{0,1},{-1,-1}})->"Spiral sink",({{0,1},{-1,0}})->"Center",({{1,0},{0,-2}})->"Saddle"}
}
]

Visualize Bifurcations using Manipulate[..]:

Consider a system of differential equations which depends on a real parameter α, say

({{{, ({{x'=-x+α*y, (1)}, {y'=-x-y, (2)}})}})

To see how the system changes when we change the parameter α in the interval [-4,4], we  use the command 
Manipulate[Plotting-Function-which-depends-on-parameter-α, {α,-4,4}]

Manipulate[ 
StreamPlot[{-x+α*y,-x-y},{x,-3,3},{y,-3,3},StreamScale->Large,StreamColorFunction->"Rainbow",ImageSize->Large],
{α,-4,4}]

We do the same thing, but plotting the eigenvalues of the matrix, together with the Phase Portrait.

Manipulate[Row[{Column[{Text["EIGENVALUES"], N[Eigenvalues[{{-1,α},{-1,-1}}]]}], StreamPlot[{-x+α*y,-x-y},{x,-3,3},{y,-3,3},StreamScale->Large,StreamColorFunction->"Rainbow",ImageSize->Large]}],{α,-4,4}]