Qualitative Behavior of Non-linear Systems:  A Competing Species Example

Suppose that x(t) and y(t) are the populations (counted in , say, millions) of two different species living in the same environment at time t.Assuming that the two species compete for food, water, and space,one might imagine that the populations are governed
roughly by the system of ODE's 

   x' = x(1 - x- y)
   y' = y(2-4x-y)
	
A fundamental idea for studying the qualitative behavior of solutions to nonlinear systems such as this one is to first determine the constant solutions, then consider the behavior of solutions near these constant solutions, and finally describe the entire phase portrait in terms of its behavior relative to the constant solutions.  We proceed by
first setting up the system and finding  its critical points with the help of a direction field.

\(\(Needs["Graphics`PlotField`"]
\)\)

f[x_, y_]  =  x (1 - x - y)\(g[x_, y_]\  = \ y \((2 - 4  x - y)\)	\)

\(\(PlotVectorField[{f[x, y], g[x, y]}, {x, 0,  .5}, {y,  .5, 1}, 
      Axes -> \ False, 
	ScaleFunction -> \((1&)\), 
      Ticks -> None, Frame -> True, AspectRatio -> 1, \ 
      PlotPoints\  -> \ 20]; \)\)

Why was the above choice of ranges for x and y a reasonable one?  Think in terms of the situation we are trying to model.
		
Critical points can be approximately detected with a direction field
by looking for those places in which the vector field seems to be turning quite a lot.  In this picture there appear to be three
such places--roughly speaking, the points (.4,.5), (0,1.75),
and (.9,0).  We can improve these guesses with a little help:

FindRoot[{f[x, y]  ==  0,  g[x, y]  ==  0},  {x,  .4}, {y,  .5}] 
    FindRoot[{f[x, y]  ==  0,  g[x, y]  ==  0},  {x,  .9}, {y, 0}] 
    \(FindRoot[{f[x, y]\  == \ 0, \ g[x, y]\  == \ 0}, \ {x, 0}, {y, 1.75}]\ 
    
\)

That is, (1/3, 2/3), (1,0), and (0,2) are the critical points.  Note that we missed the critical point (x,y) = (0,0).
	
Now we want to examine the behavior of solutions near each critical
point in turn.  First we pick on (1/3,2/3).  The direction field nears this point looks like:

\(\(PlotVectorField[{f[x, y], g[x, y]}, {x, 1/3 -  .01, 1/3\  +  .01}, {
        y, 2/3 -  .01, 2/3 +  .01}, 
	AxesOrigin -> {1/3, 2/3}, 
	
      Axes -> \ True, 
	ScaleFunction -> \((1&)\), Ticks -> None, 
      Frame -> True, AspectRatio -> 1, \ PlotPoints\  -> \ 20]; 
    \)\)

This looks a bit like the vector field of a linear homogeneous
constant coefficient system with one positive and one negative
eigenvalue.  Now we look at a phase portrait.  Since the system
is non-linear, we must use NDSolve to generate solution curves.
We define a function to produce these solution given (in order)

     starting time     t1
     stopping time     t2
            x(0)                  a
            y(0)                  b

Note that this will be most meaningful if t1<0<t2.  Also notice
that the semi-colon is important in the following definition.
Without it, Mathematica would try to evaluate the NDSolve command
and then crash when it didn't have definite values for t1,t2,a, and
b.

\(\(soln[t1_, t2_, a_, b_]\  := 
      \ {x[t], y[t]} /. 
        First[
			
          NDSolve[\ {\(x'\)[t] == f[x[t], y[t]], \(y'\)[t] == g[x[t], y[t]], 
              
					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x[0]\  == \ a, \ 
              y[0]\  == \ b}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {x[t], y[t]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, t1, t2}]
			]
          
\)\)

\(\(
\(Show[


      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			
            soln[\(-4\), 4, 1/3 + s, 2/3 +  .01], {s, \(- .01\),  .01,  .001}
              
		]], 
	\ \ \ {t, \(-4\), 4}
], 

      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			
            soln[\(-4\), 4, 1/3 + s, 2/3 -  .01], {s, \(- .01\),  .01,  .001}
              
		]], 
	\ \ \ {t, \(-4\), 4}
], 

	
      AspectRatio -> Automatic, 
	AxesOrigin -> {1/3, 2/3}, 
	
      PlotRange -> {{1/3 -  .01, 1/3 +  .01}, {2/3 -  .01, 2/3 +  .01}}, 
	
      Ticks -> None, 
	Frame -> True	\ \ \ 
]; \)			\)\)

Again, this looks very much like the picture of linear homogeneous
           system of saddle type .   In particular, solutions that begin close to (1/3, 2/3) tend to move 
          away from this point when t gets large .   For this 
          reason, we refer to (1/3, 2/3) as an unstable  critical point .   To see why the phase portrait near this point looks the way it does, we compute a linearization of the  system about the point
           (1/3, 2/3) .  That is, \(near\ this\ point\ we\ 
          \(have : 
		

            \((x - 1/3)\)'\  ≈ \ 
              \(f\_x\) \((x - 1/3)\)\  + \ \(f\_y\) \((y - 2/3)\)\)\)(y - 2/3)'  ≈  
          (Subscript[g, x]) (x - 1/3)  +  (Subscript[g, y]) (y - 2/3)
where all derivatives on the right side are evaluated at 
          the point (1/3, 2/3) .   In matrix form,  
        this is written
(({{x' (t) - 1/3}, {y' (t) - 2/3}}))  = (({{(Subscript[f, x]), (Subscript[f, y])}, {(Subscript[g, x]), (Subscript[g, y])}}))(({{x (t) - 1/3}, {y (t) - 2/3}})).

The matrix on the right side is the Jacobian matrix of the vector function  f(x, y) .  We compute this and examine the corresponding linear system .

\(Jac[x_, y_]\  = 
      \ {{D[f[x, y], x], D[f[x, y], y]}, 

        		{D[g[x, y], x], D[g[x, y], y]}}; 
MatrixForm[Jac[x, y]]\)

A  =  Jac[1/3, 2/3]; MatrixForm[A]

\(\(PlotVectorField[
      A . {x, y}, {x, \(- .01\),  .01}, {y, \(- .01\),  .01}, 
	
      Axes -> \ True, 
	ScaleFunction -> \((1&)\), Ticks -> None, 
      Frame -> True, AspectRatio -> 1, \ PlotPoints\  -> \ 20]; 
    \)\)

This looks exactly like the direction field picture for the orginal system near (1/3,2/3).

rhs =  A . {x[t], y[t]};  MatrixForm[rhs]

\(system\  = 
      \ {\(x'\)[t]\  == \ Part[rhs, 1], \ 
        \(y'\)[t]\  == 
			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Part[rhs, 2], 
        \ x[0] == a, y[0]\  == \ b}\)

\(linsoln[a_, b_]\  = \ 
      Simplify[{x[t], y[t]} /. 
          First[
		\ DSolve[system, {x[t], y[t]}, t]
		]]\)

\(\(
\(Show[


      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			linsoln[s,  .01], {s, \(- .01\),  .01,  .001}

            		]], 
	\ \ \ {t, \(-4\), 4}
], 

      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			linsoln[s, \(- .01\)], {s, \(- .01\),  .01,  .001}

            		]], 
	\ \ \ {t, \(-4\), 4}
], 
		
      AspectRatio -> Automatic, 
	AxesOrigin -> {0, 0}, 
	
      PlotRange -> {{\(- .01\),  .01}, {\(- .01\),  .01}}, 
	
      Ticks -> None, 
	Frame -> True	\ \ \ 
]; \)			\)\)

The phase portrait also matches with that of the original system.

Now we turn to the other critical points of our example, examining these in somewhat less detail than the first one.

(x,y) = (0,0).

The Jacobian matrix of (f,g) at this point is

A  =  Jac[0, 0];  MatrixForm[A]

If we look at the eigenvalues of A

Eigenvalues[A]

we see that the linear system associated to the point (0,0) has a source at the origin.  In  particular, all solutions that begin near the origin move away as t increases and (0,0) is an unstable critical point of the orginal system.  This can be seen by looking at a direction field near (0,0):

\(\(PlotVectorField[{f[x, y], g[x, y]}, {x, \(- .01\),  .01}, {y, 
        \(- .01\),  .01}, 
	AxesOrigin -> {0, 0}, 
	Axes -> \ True, 
	
      ScaleFunction -> \((1&)\), Ticks -> None, Frame -> True, 
      AspectRatio -> 1, \ PlotPoints\  -> \ 20]; \)\)

(x,y) = (1,0)

A  =  Jac[1, 0]; MatrixForm[A]

Eigensystem[A]

This time we have a sink in the associated linear system.  Thus (1,0) is a stable critical point.

\(\(PlotVectorField[{f[x, y], g[x, y]}, {x, 1 -  .01, 1\  +  .01}, {y, 
        \(- .01\),  .01}, 
	AxesOrigin -> {1, 0}, 
	Axes -> \ True, 
	
      ScaleFunction -> \((1&)\), Ticks -> None, Frame -> True, 
      AspectRatio -> 1, \ PlotPoints\  -> \ 20]; \)\)

(x,y) =(0,2)

A  =  Jac[0, 2];  MatrixForm[A]

Eigensystem[A]

Again, a sink.  The point (0,2) is also a stable critical point.

Now we know more or less how solutions of our system behave near each critical point.  It remains to understand the global behavior of the
system.  To this end we plot some trajectories on a larger scale.  We will take our initial conditions in the first quadrant along the line 2x + y = 1.

\(\(Show[
	

      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			soln[\(-2\), 3, s, 1 - 2  s], {s, 0,  .5,  .025}

            		]], 
	\ \ \ {t, \(-2\), 3}
], 

      ParametricPlot[
	\ \ 
        Evaluate[
          Table[
			soln[\(-2\), 3, s, 4 - 2  s], {s, 0, 2,  .05}

            		]], 
	\ \ \ {t, \(-2\), 3}
], 
	AspectRatio -> 1, 
	
      PlotRange -> {{0, 1}, {0, 2}}, 
	Ticks -> None, 
	
      Frame -> True	\ \ \ 
]; \)\)

Notice that there are two distinguished  solution curves implied in this picture--both passing through the saddle point at (1/3,2/3).  One of the curves joins the two sinks (1,0) and(0,2).  The other joins the source (0,0) to the saddle (1/3,2/3) and seems to separate solutions that converge to (1,0) as t goes to infinity from those that converge to(0,2).  These curves are called separatrices for the system.  There is a general theorem due to Poincare
and Bendixson that describes the possible limiting behaviors of a 2 x 2 system in terms of separatrices and critical points.  This example illustrates some of the flavor of that theorem.