Math 331.1 •  Fall 2005

Vector fields of systems of 1st order ODEs

Copyright © 2005 Murray Eisenberg.  All rights reserved.

Preparation—load package

The function PlotVectorField used in this notebook is not one of the built-in Mathematica functions.  Rather, it is in one of the standard add-on packages—files containing sets of related, special-purpose function definitions.

If, when you opened this notebook, you already clicked the Yes button when asked if you wanted to evaluate its initialization cells, then you do not have to evaluate the next cell.  Otherwise, you must evaluate the next cell.

Needs["Graphics`PlotField`"]

If ever you want to use PlotVectorField without opening this notebook, it is essential to first evaluate the expression:
	
	Needs["Graphics`PlotField`"]

Example 1: simple harmonic oscillator

In this example you examine graphically the vector field associated with the second-order ODE for a simple harmonic oscillator:

	x ' '  =  - x.
	
The equivalent system of two first-order ODEs is:

	\(TraditionalForm\`x\ '\  = \ 
        v, 
\(\(v\ '\)\(\ \)\(=\)\(\ \)\(-x\)\(\ \)\)\)

Hence the associated vector field is:

	F(x,  v)  =  (v,  -x)

The first attempt at plotting this vector field is the following:

PlotVectorField[{v,  -x},  {x,  -4,  4}, {v,  -4,  
          4},  Axes  ->  True,  AxesLabel  ->  {x,  v},  
        ScaleFactor  ->  None];

To make that a bit simpler to understand, use a smaller number of points in the grid at which the vectors are shown:

PlotVectorField[{v,  -x},  {x,  -4,  4}, {v,  -4,  
          4},  Axes  ->  True,  AxesLabel  ->  {x,  v},  
        ScaleFactor  ->  None,  
        PlotPoints  ->  9];

The resulting plot is still a bit messy—the arrows are so long they overlap.  If you omit the ScaleFactor option (which is same as using the default value ScaleFactor -> Automatic, then the arrows are shortened by Mathematica:

vfield  =  
        PlotVectorField[{v,  -x},  {x,  -4,  4}, {v,  -4,  
            4},  Axes  ->  True,  
          AxesLabel  ->  {x,  v}];

Now the arrows are too short.

For such reasons, we usually plot the direction field of the system rather than its actual vector field.  In the direction field, the arrows all have the same length.

dirField  =  
        PlotVectorField[{v,  -x},  {x,  -4,  4}, {v,  -4,  
            4},  Axes  ->  True,  AxesLabel  ->  {x,  v},  
          ScaleFunction  ->  (0.75 &)];

In the input cell above, the expression (0.75&) stands for a function that has constant value 0.75.  That is the length each arrow has.

We are going to plot together a solution in the phase plane together with the direction field.  This will allow us to see how the direction field does, in fact, show the form of the solution curve.

For initial conditions 

	x(0)  =  0,          x ' (0)  =  
        2, 
	
we know the solution:

x[t_]  :=  2 Sin[t]
v[t_]  := 2  Sin'[t]

The function ParametricPlot displays the solution curve.  (The option AspectRatio -> Automatic forces the scales on the horizontal and vertical axes to be the same.)

curve  =  
        ParametricPlot[{x[t],  v[t]},  {t,  0,  4  π},  
          AspectRatio  ->  Automatic];

The initial point is the following graphics object (you will not see any graphics displayed!):

initial  =  
      Graphics[{PointSize[0.025],  Point[{0,  2}]}]

Here are the direction field, the solution curve in the phase plane, and the initial point together:

Show[{dirField,  curve,  initial}];

If you wish, in the cell defining dirField, try changing the ScaleFunction length from 0.75 to something else, or try changing the PlotPoints value.

To emphasize the solution curve, make it a different color or thickness:

redcurve  =  
        ParametricPlot[{x[t],  v[t]},  {t,  0,  4  π},  
          AspectRatio  ->  Automatic,  
          PlotStyle  ->  RGBColor[1,  0,  0]];

Show[{dirField,  redcurve,  initial}];

thickcurve  =  
        ParametricPlot[{x[t],  v[t]},  {t,  0,  4  π},  
          AspectRatio  ->  Automatic,  
          PlotStyle  ->  Thickness[0.01]];

Show[{dirField,  thickcurve,  initial}];

Example 2: another homogeneous linear system

The system of ODEs here is:

	\(TraditionalForm\`x\ '\  = \ x\  + \ 2\ y, 

      y\ '\  = \ 4\ x\  - \ y\)

The associated vector field is:

F[{x_,  y_}]  :=  {x  +  2 y,  4 x  -  y}

For example:

F[{-1,  (3/2)}]

Find equilibrium points:

Solve[F[{x,  y}]  ==  {0,  0},  {x,  y}]

Exercise 1:  Use Mathematica to plot the direction field of F on the region of the phase plane where -3  ≤  x  ≤  3 and -3  ≤  y  ≤  3.  I suggest a length of 1 for the vectors (instead of the length 0.75 used in Example 1, above).  Assign the result to the name field.

Exercise 2:  Print the direction field just plotted.  On it sketch by hand the solution curves through each of the points (-1,  (3/2)), (-1,  2), (-1,  2  +  (5/100)), and (-1,  3).  On the printout tell what seems to happen to (x(t),  y(t)) as t tends to infinity.  (Do this before continuing with this notebook.)

Suggestions.  Before printing the plot, enlarge it by clicking anywhere in the plot and then dragging one of its corner "handles" outward.  Include in what you print the input statement that produced the plot.  Type your name(s) in a cell just before that input cell, so you can identify your printout.

Note.  To print just a few adjacent cells, hightlight their cell brackets by dragging the mouse over them.  Then from the File menu choose Print Selection.

Now let's let Mathematica find some solutions of this system of ODEs.  You could do it directly in terms of the pair of functions x(t) and y(t), like this:

Clear[x,  y,  t]

DSolve[{x'[t]  ==  x[t]  +  2 y[t],  y'[t]  ==  
          4 x[t]  -  y[t],  x[0]  ==  -1,  
        y[0]  ==  (3/2)},  {x[t],  y[t]},  t]

Again the result is a pair of replacement rules, one for x(t) and the other for y(t).  Actually, the pair of rules is nested inside an outer pair of braces, which we remove by using First:

soln1  =  
      First[ DSolve[{x'[t]  ==  x[t]  +  2 y[t],  y'[t]  ==  
              4 x[t]  -  y[t],  x[0]  ==  -1,  
            y[0]  ==  (3/2)},  {x[t],  y[t]},  t] ]

To get actual functions of t, use the "replace" operator ( /. ) as usual:

{x[t],  y[t]}  /.  soln1

We plot the solution curve in the phase plane parametrically for, say, 0  ≤  t  ≤  1.  And we color the curve red to distinguish it from other solution curves that you will create with Mathematica below.

curve1  =  
        ParametricPlot[{x[t]  /.  soln1,  y[t]  /.  soln1},  {t,  0,  
            1},  PlotRange  ->  {{-3,  3},  {-3,  3}},  
          PlotStyle  ->  RGBColor[1, 0, 0]];

Here is that solution curve shown on the direction field:

Show[{field,  curve1}];

Exercise 3:  (a) Similarly find the solutions soln2, soln3, soln4 for each of the sets of initial conditions:
	
	(2)	x(0)  =  -1,  y(0)  =  2
	(3)	x(0)  =  -1,  y(0)  =  2 +  (5/100)
	(4)	x(0)  =  -1,  y(0)  =  3
	
(b) Display the plots curve2, curve3, curve4 in the phase plane of each of the corresponding  solution curves.  Use different colors for each so that you can distinguish them from one another and from curve1.  (You will need to adjust the final value of t in the parametric plots so as to keep the curve within the square region -3  ≤  x  ≤  3 and -3  ≤  y  ≤  3.  Some experimentation will be required to do that.)

(c) Next, Show everything together—the direction field field and all four solution curves.

(d) Finally, compare the long-term behavior (as t  tends to infinity) of these solution curves, paying special attention to the curves curve2 and curve3 whose initial points were near one another.