MATLAB TUTORIAL, part 2.3: Planar Systems

Planar Systems
In first part of the course, we discussed the direction field for first order equations. The construction of a direction field is equally useful in the study of autonomous systems of two first-order equations. It provides an overall view of where the solution curves go, and the arrows show which way the system moves as time increases. We define here matlab code for plotting such direction fields.

We consider systems of first order ordinary differential equations in normal form on the plane
\[ \begin{cases} \dot{x} &= f (x , y) , \\ \dot{y} &= g (x , y) , \end{cases} \]
where we use notation x and y for dependent variables. Of course, the above system can be written in compact vector form
\[ \frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x} ) , \]
\[ {\bf x} (t) = \begin{bmatrix} x (t) \\ y (t) \end{bmatrix} , \qquad {\bf f} (x , y) = \begin{bmatrix} f (x , y) \\ g (x , y) \end{bmatrix} \]
are 2-column vectors. In engineering and physics, it is a custom to denote a derivative with respect to time variable t by dot: \( \dot{\bf x} = {\text d}{\bf x} / {\text d} t. \)

If an initial position of the vector \( {\bf x} (t) \) is known, we get an initial value problem:

\[ \frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}({\bf x} ) , \qquad {\bf x} (t_0 ) = {\bf x}_0 , \]
where \( {\bf x}_0 \) is a given column vector. Recall that an autonomous differential equation is a system of ordinary differential equations which does not depend on the independent variable (t in our case). It is of the form

\[ \dot{\bf x} = {\bf f} ( {\bf x} ) . \]
Here dot stands for the derivative with respect to time variable t and \( {\bf f} ({\bf x}) \) takes values in n-dimensional Euclidean space and t is usually time.


2.4.1. Planar Case


Consider the following non-linear system of ordinary differential equations

\[ \dot{x} = \left( x-2y \right) \left( x+2y-4 \right) , \qquad \dot{y} = x \left( 1+y \right) . \]
The critical points are found by solving the algebraic equations
\[ \left( x-2y \right) \left( x+2y-4 \right) =0 , \qquad x \left( 1+y \right) =0 . \]
One way to satisfy the second equation is by choosing x = 0. Then the first equation becomes \( 2y \left( 2y-4 \right) =0 , \) so y = 0 or y = 2. More solutions can be found by choosing y = -1 in the second equation. Then the first equation becomes \( \left( x+2 \right) \left( x-6 \right) =0 , \) so x = -2 or x = 6. Thus, we have obtained the four critical points: (0,0), (0,2), (-2,-1), and (6,-1).

Figure shows a direction field ■

For a two-dimensional autonomous system with at least one asymptotically stable critical point, it is often of interest to determine where in the phase plane the trajectories lie that ultimately approach a given critical point. Let P be a point in the xy-plane with the property that a trajectory passing through P ultimately approaches the critical point as \( t \mapsto \infty . \) Then this trajectory is said to be attracted by the critical point. Further, the set of all such points is called the basin of attraction or the region of asymptotic stability of the critical point. A trajectory that bounds a basin of attraction is called a separatrix because it separates trajectories that approach a particular critical point from other trajectories that do not do so. Determination of basin of attraction is important in understanding the large-scale behavior of the solutions of an autonomus system.

Example. Consider the system

\[ \dot{x} = \left( 3+x \right) \left( y-2x \right) , \qquad \dot{y} = y \left( 2+x-x^2 \right) . \]
Equating slope functions to zero and solving corresponding algebraic equations, we find critical points to be (0,0), (-1,-2), (3,6), and (-3,0). ■

The trajectoreis of a two-dimensional autonomous system

\[ \dot{x} = f\left( x,y \right) , \qquad \dot{y} = f \left( x,y \right) \]
can sometimes be found by solving a related first order differential equation
\[ \frac{{\text d}y}{{\text d}x} = \frac{g(x,y)}{f(x,y)} , \]
which is a first order equation in the variables x and y. Observe that such a reduction is not usually possible for nonautonomous system when f and g depend also on time varibale t. If the above equation \( {\text d}y / {\text d}x = g(x,y) / f(x,y) \) can be solved by any of the methods of Part I, and if we write solutions (implicitly) in the form
\[ H(x,y) = c, \]
with some arbitrary constant c, then it defines trajectories of the given autonomous system. In other words, the trajectories lie on the level curves of H(x,y). Keep in mind that there is no general way of solving first order equation \( y' = g(x,y)/f(x,y) \) to obtain the function H(x,y), so this approach is applicable only in special cases.

Example. Find trajectories of the system

\[ \dot{x} = 2-y , \qquad \dot{y} = 9 - x^2 . \]
From the equations
\[ 2-y =0 , \qquad 9 - x^2 =0, \]
we find that the critical points of the given system are the points (3,2) and (-3,2). To determine the trajectories, note that the system is equivalent
\[ \frac{{text d}y}{{text d}x} = \frac{9 - x^2}{2-y} . \]
Separating the variables and integrating, we find that solutions satisfy
\[ H(x,y) = 2y - \frac{y^2}{2} - 9x + \frac{x^3}{3} = c, \]
where c is an arbitrary constant. A computer plotting routine is helpful in displaying the level curves of H(x,y). The point (-3,2) is center anad (3,2) is a saddle point.


Direction fields

Since Matlab has no friendly subroutine to plot direction fileds for ODEs, we present several codes that allow to plot these fields directly. Many others can found on the Internet.