MATH 302-010
Ordinary Differential Equations
Prof. D. A. Edwards
Copyright 2024

Linear Algebra and ODE Systems

We begin with our standard command:

Quit[]

Since we are working with matrices and vectors, we also insert the following command, which makes them display as they do in class:

$PrePrint=If[MatrixQ[#] || VectorQ[#],MatrixForm[#],#]&;

Basic Operations

A vector in Mathematica can be entered as a list:

xvec={2,1,3}
yvec={4,3,2}
zvec={2,5}

To perform vector addition, subtraction, or scalar multiplication, we treat the vectors just as if they were regular variables.  Of course, we have to make sure the dimensions match.

3*xvec
tvec=t*xvec
wvec=xvec-2*yvec
yvec+zvec

Differentiation of a vector works just like for a scalar:

D[tvec,t]

In Mathematica, a matrix is defined as a list of lists.  Each internal list corresponds to a row:

a1={{2,3},{1,-5}}
a2={{1,0},{-1,3}}

Matrix addition is straightforward.  You add matrices just as you would numbers, as long as their dimensions match:

a1+a2

For scalar-matrix multiplication, you again use the * command:

3*a1-6*a2

To perform matrix-vector multiplication, we use the . command instead of the * command, which gives us a nonsensical result:

a1.zvec
a1*zvec

We may also need to calculate the inverse of a matrix, which is done in Mathematica using the Inverse command:

Inverse[a1]

Exercises

1.  Calculate the sum of r=(({{1
-2}, {3}})) and s=(({{0
2}, {-4}})).

Solution

r = {1,-2,3}
s={0,2,-4}
rps = r+s

2.  Calculate the inverse of (({{3, -1}, {4, 4}})).

Solution

mat2={{3,-1},{4,4}}
Inverse[mat2]

Solving Linear Systems

There are two basic ways to solve systems of linear equations in Mathematica.  First, we can input the equations one by one, being careful to use two equals signs in each equation:

eq1=2*x+3*y==0
eq2=x-5*y==4

Once we have entered the equations, we may use the Solve command as in MATH 242.  The difference is that we are solving a set of equations for a set of variables, so that is the syntax we must use:

Solve[{eq1,eq2},{x,y}]

However, we note that this system of equations is the same as (Subscript[A, 1])x=b, where (Subscript[A, 1]) is given above.  Therefore, we can also use the inverse, once we define b:

b1={0,4}
Inverse[a1].b1

If the system has an infinite number of solutions (which is the case with solving for eigenvectors), Solve will show that.  Consider the following system:

eq4=x-2*y==3
eq5=-2*x+4*y==-6
Solve[{eq4,eq5},{x,y}]

Here we see that there is a free variable.  There is a command (LinearSolve) that works directly with the matrix-vector form of the system, but it is not good at solving equations with free variables, so we don’t use it.

Exercise

1.  Solve the following system:
                         3 x + 2 y = 7
4 x - 5 y = -11

Solution

prob1a=3*x+2*y==7
prob1b=4*x-5*y==11
Solve[{prob1a,prob1b},{x,y}]

Determinants, Real Eigenvalues, and Eigenvectors

To calculate the determinant of a matrix, we use the Det[Matrix] command:

a={{2,1},{2,2}}
Det[a]

To calculate just the eigenvalues of a matrix, we use the Eigenvalues command:

evec=Eigenvalues[a]

The eigenvalues are listed in vector format. To get a particular eigenvalue, we just use the Vector[[i]]notation to pick out a value:

evec[[2]]

If you want to know the eigenvectors as well, use the Eigensystem command, which returns both the eigenvalues and the eigenvectors:

esys = Eigensystem[a]

The eigenvalues are stored in the first row, and the corresponding eigenvectors are stored in the second row.  The eigenvectors may be listed in any order, and any multiple of an eigenvector is also an eigenvector, so Mathematica’s results may not be the same as what you do by hand.  Separating out the component parts, we have

lvec=esys[[1]]
zmat=esys[[2]]

To select out an eigenvector, we must pick a row of zmat, which again we do using the [[i]] notation:

z2=zmat[[2]]

Note that the result is back in the proper vector format.  We can check that this is an eigenvector by multiplying directly:

l2=lvec[[2]]
a.z2
Simplify[l2*z2]

Exercise

1. Find the determinant, eigenvalues, and eigenvectors of (({{0, 4}, {13, 9}})).

Solution

mat1={{0,4},{13,9}}
Det[mat1]
Eigensystem[mat1]

Solutions of Homogeneous Equations

Mathematica solves systems of ODEs by a combination of the techniques used to solve systems of algebraic equations and single ODEs.  In particular, if given a set of equations:

sys1a=x'[t]==3*x[t]+2*y[t]
sys1b = y'[t]==-2*x[t]-2*y[t]

use DSolve, listing both variables to be found as a set, as shown below:

sol1=DSolve[{sys1a,sys1b},{x[t],y[t]},t]

Note that two arbitrary constants appear, but they appear in BOTH functions, since they are related to one another by the system.  If you have particular initial conditions to solve, you can add them---as many or as few as you have:

ic1a=x[0]==3
DSolve[{sys1a,sys1b,ic1a},{x[t],y[t]},t]

ic1b=y[0]==5
DSolve[{sys1a,sys1b,ic1a,ic1b},{x[t],y[t]},t]

We can check if our solution makes sense by looking at the eigenvalues of the underlying matrix:

mat1 = {{3,2},{-2,-2}}
Eigenvalues[mat1]

Since the eigenvalues are 2 and -1, the exponents should be as well.  This can be verified by expanding the solution:

Expand[sol1]

Exercise

1.  Find the general solution of the following system:

x'=x-y
y'=-2x+7y

Then find the solution subject to x(0)=3, y(0)=5.

Solution

p1eq1 = x'[t]==x[t]-y[t]
p1eq2 = y'[t]== -2*x[t]+7*y[t]
DSolve[{p1eq1,p1eq2},{x[t],y[t]},t]

p1ic1 = x[0]==3
p1ic2 = y[0]==5
DSolve[{p1eq1,p1eq2,p1ic1,p1ic2},{x[t],y[t]},t]

Phase Plane (Real Eigenvalues)

To plot the phase plane, we normally would like to see trajectories and direction fields.  This can most easily be handled by the StreamPlot command.

IMPORTANT: To use StreamPlot, we create a vector of the right-hand side of each equation, and eliminate any t-dependence.  For our previous system, the right-hand sides are

rhs={3*x+2*y,-2*x-2*y}

Once we have the right-hand side, we use StreamPlot[equations,xrange,yrange]:

StreamPlot[rhs,{x,-1,1},{y,-1,1}]

This gives us a direction field, and since the arrows line up tip to tail, we see that we have trajectories as well.

Exercise

1.  Graph the phase plane of the following system:

x'=x-y
y'=-2x+7y

Solution

We write the system in the form that StreamPlot needs:

p1r1 = x-y
p1r2 = -2*x+7*y
StreamPlot[{p1r1,p1r2},{x,-2,2},{y,-2,2}]