Demo # 11
Systems of ODEs in Mathematica

Mathematica can be used to solve systems of equations.  consider the system
	(dx/dt) = -(1/4) x + 2  y
	(dy/dt) = -2 x - (1/4) y
with the initial conditions
	x(0) = -1
	y(0) = 0
We can solve these equations within Mathematica using DSolve:

Clear[x, y];
deq1 = x'[t] == -(1/8) x[t] + 2 y[t];
deq2 = y'[
            t] == -2 x[t]  - (1/4) y[t];

ic1 = x[0] == -1;
ic2 = y[0] == 0;

soln = DSolve[{deq1, deq2, ic1, ic2}, {x[t], y[t]}, t]

These solutions look pretty complicated.  Let's plot them to see what is going on

Plot[Evaluate[x[t] /. soln], {t, 0, 10}]

Plot[Evaluate[y[t] /. soln], {t, 0, 10}]

So the solutions are a pair of damped oscillating solutions that appear to be 90 degrees out of phase with each other.  We can make a phase portrait by making a parametric plot of the y vs x solutions.

plt1 = 
      ParametricPlot[Evaluate[{x[t], y[t]} /. soln], {t, 0, 10}, 
        AspectRatio -> 1, Frame -> True, Axes -> True, 
        AxesLabel -> {"x", "y"}]

<< Graphics`PlotField`

We can gain further insight into the solutions by plotting the direction field of the equations.  The direction field is obtained by eliminating t between the two ODEs
	\(TraditionalForm\`\(dy/dt\)\/\(dx/dt\) = \(dy\/dx = \(\(-2\) x - y/4\)\/\(\(-x\)/4 + 2  y\)\)\)
and can be plotted in Mathematica using PlotVectorField[{dx,dy}, ...]

plt2 = 
      PlotVectorField[{(-x/4 + 2 y), (-2 x  - 
              y/4)}, {x, -1, 1}, {y, -1, 1}, Frame -> True, 
        Axes -> True, AxesLabel -> {"x", "y"}]

Let's see how this lines up with particular solution we found by plotting the trajectory for the particular solution on top of the dirction field.

Show[plt1, plt2]

This seems to be a pretty powerful method for visualizing the solutions.  Let's look at what happens for some other choices of coefficient matrices.

$Post := If[MatrixQ[#1],  MatrixForm[#1],  #1] &;

M = {{-(1/8), 2}, {-2 ,  -(1/4)}}

Let's find the eigenvalues of this matrix

Eigenvalues[M]

... and their numerical values

Eigenvalues[M] // N

The eigenvalues are complex with negative real parts.  The imaginary parts indicate an oscillatory soluton while the negative real parts indicate an exponentially decaying solution.  So we expect damped oscillations moving towards zero.  This is consistent with what we found in our plots above.

Eigenvectors[M]

Let's try another coefficient matrix:

a = {{0, 1}, {2 ,  3}}
Eigenvalues[a]
Eigenvalues[a] // N
Eigenvectors[a]
plt3 = 
      PlotVectorField[{ y, 2  x  + 3  y}, {x, -1, 1}, {y, -1, 1}, 
        Frame -> True, Axes -> True, 
        AxesLabel -> {"x", "y"}]

This time we have purely real eigenvalues and eigenvectors.  One eigenvalue is positive and the other is negative.  This indicates a saddle point, as is verified in the plot above.

The eigenvectors give the directions that show the flow directly to and from the fixed point at the origin.  Let's plot them up to see how this works graphically.

<< Graphics`Arrow`

evec1a = 
        Arrow[{(1/4) (-3 - (Sqrt[17])), 1}, {0, 
            0}];
evec1b = 
        Arrow[{-(1/4) (-3 - (Sqrt[17])), -1}, {0, 
            0}];

plt4 = Show[Graphics[{Hue[0], evec1a, evec1b}]]

evec2a = 
        Arrow[{0, 0}, {(1/4) (-3 + (Sqrt[17])), 
            1}];
evec2b = 
        Arrow[{0, 
            0}, {-(1/4) (-3 + (Sqrt[17])), -1}];

plt5 = Show[Graphics[{Hue[0.7], evec2a, evec2b}]]

Show[plt3, plt4, plt5]

As expected, the eigenvectors identify the flows to an from the fixed point, which is a saddle point in this example.

Let's do another example.

b = {{0, -1}, {2 ,  0}}
Eigenvalues[b]
Eigenvalues[b] // N
Eigenvectors[b]
pltA = 
      PlotVectorField[{- y, 2  x }, {x, -1, 1}, {y, -1, 1}, 
        Frame -> True, Axes -> True, 
        AxesLabel -> {"x", "y"}]

Purely imaginary eigenvalues indicate undamped oscillating solutions.  The circulating phase plot is consistent with this.  We can check it by finding a particular solution.  From the phase plot above, it looks like an convenient initial condition would be x(0) = 0 and y(0) = 1.

solnA = 
      DSolve[{x'[t] == -y[t], y'[t] == 2 x[t], 
          x[0] == 0, y[0] == 1}, {x[t], y[t]}, t]

pltB = 
      ParametricPlot[Evaluate[{x[t], y[t]} /. solnA], {t, 0, 10}, 
        AspectRatio -> 1, Frame -> True, Axes -> True, 
        AxesLabel -> {"x", "y"}]

Show[pltA, pltB]

As expected, the orbit is closed.