Toronto Math Forum

PREFACE

In this tutorial, you will learn how to create a graph like the following one, which is the phase portrait for the system
\begin{equation*}
\left\{ \begin{array}{l} x'(t) = x\ln{(y^2)}, \\
y'(t) = y\ln{(x^2)}.
\end{array} \right.
\end{equation*}

Note that the right-hand expressions are not continuously differentiable at $(0,0)$ and there is no linearization at this point, which explains that while stable, the point $(0,0)$ does not look (after close examination) like a standard node.

In order to create such plots, you need computational software and the ability to write expressions that the software can interpret. There are three popular programs for this purpose: Waterloo Maple, MATLAB, and Wolfram Mathematica. However, I should caution you that these programs are expensive, even if you are buying a student license. There are a lot of open-source alternatives: Mathics, Sage, GNU Octave (similar to MATLAB), and Scilab, to name a few. Each software has its advantages and disadvantages. For example, it is commonly agreed that  Mathematica is good for symbolic manipulation of expressions (e.g. integration, solving differential equations, etc.) whereas MATLAB is more appropriate for working with large sets of data. This is a decision that I will leave to you. This guide is concerned with Mathematica.

1 DIRECTION FIELDS

A direction field for a first-order linear ordinary differential equation can be created using the VectorPlot function. For our purposes, the basic form is

VectorPlot[{1, <function>}, {x, <xmin>, <xmax>}, {y, <ymin>, <ymax>}]
Here, <function> is the differential equation in the form $y'(x) = f(x,y)$. By default, the direction field that you get will be a coloured vector field. If you want to get a classic direction field, you'll need to add the option VectorStyle -> Arrowheads[0]. Often, the colour unnecessarily increases image size, and if you intend the plot your solution on the direction field, then it is better to keep the direction field in black and the solution in colour. To do this, you can add the VectorStyle -> Black option. To see how these options are used, see the following example for the differential equation $y' = xy$.

VectorPlot[{1, x*y}, {x, -1, 1}, {y, -1, 1}, VectorStyle -> {Arrowheads[0], Black}]

2 PHASE PORTRAITS

A phase portrait for a two-dimensional autonomous system of ordinary differential equations can be created using the StreamPlot function. The generic form of what you will type will look like this:

StreamPlot[{<function1>, <function2>}, {x, <xmin>, <xmax>}, {y, <ymin>, <ymax>}]
Here, <function1> and <function2> correspond to $x'(t)$ and $y'(t)$, respectively. As with VectorPlot, you can add some options after these arguments. I often use only StreamStyle -> Black since the default colour is unnecessary. To exemplify this function, I will show you how to plot the phase portrait for case $(d)$ of the system in the third semester-end challenge.

StreamPlot[{-y + x*(-(Sin[Pi*Sqrt[x^2 + y^2]])^2), x + y*(-(Sin[Pi*Sqrt[x^2 + y^2]])^2)}, {x, -5, 5}, {y, -5, 5}, StreamStyle -> Black]

You can see that the integral curves on this phase portrait are not very long. To increase the length, you can make use of the StreamScale -> <size> option, where <size> is a real number in $[0,1]$. The length that you specify depends on the phase plane and on your intentions. Applying the StreamScale option not only scales the length of the curves but also the arrowheads. This is not what we want, and we can add the Arrowheads[<size>] option to avoid this. It is usually sufficient to use a size of Tiny or Small. In the following example, you can see a phase portrait of the very first system in this tutorial. The main difference is that the integral curves are much longer.

StreamPlot[{x*Log[y^2], y*Log[x^2]}, {x, -1, 1}, {y, -1, 1}, StreamScale -> 1, StreamStyle -> {Arrowheads[Small], Black}]

3 CONTOUR MAPS

A contour map for a system is made with the ContourPlot function. Its form is not very different from the above functions:

ContourPlot[{<implicitfunction>}, {x, <xmin>, <xmax>}, {y, <ymin>, <ymax>}]
Note that <implicitfunction> is simply a function in the form $H(x,y)=C$. Thus, ContourPlot is useful only if we can find $H$. For example, it was possible in the Easter challenge to determine $H$ for system $(a)$ to be $H(x,y) = \cos{x}+\cos{y}$. We can plot this with the expression

ContourPlot[{Cos[x] + Cos[y]}, {x, -2Pi, 2Pi}, {y, -2Pi, 2Pi}, ColorFunction -> "BlueGreenYellow"]

Note that, out of personal preference, I have added the option ColorFunction -> "BlueGreenYellow".

FURTHER READING

The most important advice that I can give anyone who is willing to pursue this kind of work is to read the program documentation. The documentation for Mathematica is located below, but you can also access it from within the program. You can learn a lot by simply looking at the examples that are given.

http://reference.wolfram.com/mathematica/guide/Mathematica.html

Sometimes, reading documentation is not good enough because you may not directly make a connection between a problem and a function. Then, you should search for tutorials. Many nice tutorials can be found on university websites. Here is an example.

http://www.math.mtu.edu/~msgocken/pdebook2/mathtut2.pdf

Last, but not least, is a technical comparison of Maple, MATHLAB, and Mathematica.

http://amath.colorado.edu/computing/mmm/

LINKS

Maplesoft's Maple
http://www.maplesoft.com/products/maple/

MathWorks' MATLAB
http://www.mathworks.com/products/matlab/

Wolfram Mathematica
http://www.wolfram.com/mathematica/

Mathics
http://www.mathics.org/

Sage
http://www.sagemath.org/

Scilab
http://www.scilab.org/

GNU Octave
http://www.gnu.org/software/octave/

Gnuplot
http://www.gnuplot.info

LOG

06.04.13: First version published.
09.04.13: Updated Section 2 with StreamScale option and new example.

I have plotted both spirals. It is now easier to see that in the case of $(e)$, the system has more circular character. The trajectories must conform to a missing spiral, which I think is the spiral for which $f(r)$ is discontinuous.

Attached are the stream plots of the three systems. One can refer to Easter challenge topic to see how system $(a)$ can be characterized. When we look at system $(b)$ for both cases $\alpha = \pm 1$, we find that it has the same critical points as system $(a)$.

If $\alpha = 1$, then the stable centers in $(a)$ are spiral points in $(b)$. If in $(a)$ the closed trajectories are directed counter-clockwise, then the corresponding spiral points in $(b)$ is also directed counter-clockwise and are stable. This is observed at $(0,0)$. If the closed trajectories in $(a)$ are directed clockwise, then the corresponding spiral points in $(b)$ are also directed clockwise but are unstable. This is observed at $(\pi,\pi)$.

If $\alpha = -1$, then the same analysis applies. However, the stabilities are inverted because we have changed the sign of the coefficient.

In both systems, saddle points remain saddle points, although they are rotated.

In each case $(a)$ to $(e)$, there is one critical point $(0,0)$ that is a spiral point. Let us first make some generalizations. In cases $(a)$ to $(d)$, there is a set of limit cycles. It would seem that when $f(r) = 0$, a limit cycle occurs. Thus, the limit cycles are given by the equation $r = \ell$, where $\ell = 1,2,...,n$. In case $(e)$, there are no limit cycles.

The stream plot for $(a)$ tells us that the spiral point is unstable. We also see that the limit cycles with $r = 2k+1$, where $k = 0,1,...,n$, are asymptotically stable both internally and externally. Limit cycles with $r = 2k$, where $k = 1,2,...,n$, are unstable both internally and externally. In case $(b)$, the spiral point is stable, and the limit cycles behave in an opposite manner. This is because of the inclusion of a minus sign, which changes the direction of the trajectories.

In case $(c)$, the spiral point is unstable. In fact, all limit cycles are internally stable and externally unstable. The trajectories are only able to spiral outward from $(0,0)$--this corresponds to the fact that $\sin^2{x}$ is always positive. In case $(d)$, the spiral point is stable and the limit cycles behave in an opposite manner. This is because of presence of the minus sign.

Case $(e)$ is an outlier. No sinusoidal function appears, so no periodic behaviour is observed. As an aside, there is no trajectory that goes through the critical point, since $f(r)$ is discontinuous there.

The sine function has the following consequences in this situation:

• it generates a periodic limit cycle whenever $r = \ell$,
• the sign of its output determines the orbital stability of the limit cycle.

I attached stream plots and corresponding contour maps of $f(r)$. Light areas are positive and dark areas are negative. What I would like to continue working on is how the output of the sine function relates to the direction of the trajectories in the regions that are bounded by the limit cycles.

The critical points can be determined by solving the equations
 \begin{equation*}
 0 = y + x(1 - x^2 - y^2), \qquad 0 = -x + y(1 - x^2 - y^2).
 \end{equation*}
It turns out that the only critical point is $(0,0)$. Now, let
 \begin{equation*}
 F(x,y) = y + x(1 - x^2 - y^2), \qquad G(x,y) = -x + y(1 - x^2 - y^2).
 \end{equation*}
The Jacobian matrix of these functions is
 \begin{equation*}
 J = \left( \begin{array}{cc} 1 - 3x^2 - y^2 & 1 - 2xy \\ -1 - 2xy & 1 - x^2 - 3y^2 \end{array} \right).
 \end{equation*}
Evaluating it at $(0,0)$, we find that
 \begin{equation*}
 J = \left( \begin{array}{rr} 1 & 1 \\ -1 & 1 \end{array} \right).
 \end{equation*}
The linear system at $(0,0)$ is then
 \begin{equation*}
 \frac{d}{dt} \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{rr} 1 & 1 \\ -1 & 1 \end{array} \right)
 \left( \begin{array}{c} x \\ y \end{array} \right).
 \end{equation*}
The eigenvalues of the coefficient matrix are
 \begin{equation*}
 \text{det} \left( \begin{array}{cc} 1 - \lambda & 1 \\ -1 & 1 - \lambda \end{array} \right) = \lambda^2 - 2\lambda + 2 = 0
 \Longleftrightarrow \lambda = 1 \pm i.
 \end{equation*}
The eigenvalues are complex with positive real part. Therefore, $(0,0)$ is an unstable spiral point. We can conclude that the non-linear system behaves similarly. To determine the orientation of the spiral, we use the vector $(x,y)^T = (0,1)$ in the equation of the system to find that
 \begin{equation*}
 \frac{d}{dt} \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{rr} 1 & 1 \\ -1 & 1 \end{array} \right)
 \left( \begin{array}{c} 0 \\ 1 \end{array} \right) = \left( \begin{array}{c} 1 \\ 1 \end{array} \right).
 \end{equation*}
This tells us that the spiral is oriented clockwise. A stream plot of the system is attached. Note that there is a closed trajectory $x^2 + y^2 = 1$ that exhibits orbital stability. It is the limit cycle of the system and is approached by all trajectories as $ t\rightarrow \infty$. This can be shown rigorously by applying the same procedure that is given in Example 1 of Section 9.7--in fact, the polar system of equations that must be solved is the exact same.

I guess that for certain values of $\ell$ and $n$, we have saddle points, such as $(-\pi,0)$.

This is a nice problem. I can give a little bit of input: in system (b), the rate of change $y'(x)$ is twice that of the same rate in system (a). Perhaps we can treat (b) as a vertical expansion of (a), which is what the contour maps suggest. We can also note that each system has the same critical points $(2 \pi n,2 \pi n)$, where $n$ is an integer.

By inspection, we see that the separatrices in system (a) are lines with slopes $\pi$ and $-\pi$. In fact, by looking more closely at the contour maps, the separatrices appear to have the equations $$ y = ± \pi x + n \pi, $$ where $n$ is an integer. In system (b), the separatrices are sinusoidal functions that oscillate about an equilibrium line that is parallel to the $y$-axis. Finally, because there are infinitely many critical points in each system, there are also infinitely many separatrices.

Oh, okay. I was confused because the origin is a saddle point (as seen from the stream plot and the eigenvalues $\lambda = ±4$) and there is no apparent basin of attraction. I should note that the textbook (9th ed.) discusses separatrices only in the context of basins of attraction. Anyway, I updated the solution again and plotted more level curves for completion.