Mathematica Solution 
of the
Local Stability Conditions

for
The Competition Equations:

dx = m x (1 - (x/h) - ((b y)/k)) dtdy = n y (1 - ((c x)/h) - (y/k)) dt

by
K. D.  Stroyan
University of Iowa

copyright 1997 by Academic Press, Inc. - All rights reserved.

Notebook Overview

This notebook is used to help with computatons needed in using the local stability theorem.  This theorem states:

Let f[x,y]and g[x,y] be smooth functions with f[xe,ye)]= g[xe,ye)]= 0.  The coefficients given by the partial derivatives evaluated at the equilibrium define the characteristic equation of the equilibrium,

(ax - r)(by - r) - aybx = r2 - (ax + by)r + (axby - aybx) = 0

Suppose that the real parts of both of the roots of this equation are negative.  Then there is a real neighborhood of (xe,ye) or a real positive E>0 such that when a solution satisfies

dx/dt = f[x,y]
		dy/dt = g[x,y]
with initial condition in the neighborhood, Abs[x(0) - xe]<E and Abs[y(0) - ye]<E, then
		(lim_(t -> ∞))x(t) = xe       and       (lim_(t -> ∞)) y(t) = ye.

The stability result means that when the real parts of the characteristic roots are negative, one can look in the microscope infinitely long and never tell the difference between the linear and nonlinear flows.

References to the Text

In Section 24.2 we study linear equilibria in two dimensions.  This NoteBook uses that theory to obtain the characteristic equation.  Mathematica finds its roots.

Step-by-step Computation

Refer to Example 24.1 in the textbook and make sure you understand what each Mathematica  command is computing.  Note that in this example, as in the textbook, we use the equilibrium point (0,5) to find the roots.  You should do the computation for the other equilibrium points.

Define the Equations

Clear[f, g, x, y]; f = y (1 - (x/7) - (y/5))g = x (1 - (x/8) - (y/6))

Find the Equilibria

Solve[{f == 0, g == 0}, {x, y}]

The Symbolic Matrix of Partials

\(aX\  = \ D[f, \ x]; \ 
aY\  = \ D[f, \ y]; \ 
bX\  = \ D[g, \ x]; \ 

    bY\  = \ D[g, \ y]; \ 
Print[TableForm[{{aX, \ aY}, \ {bX, \ bY}}]]\)

The Specific Matrix at (0,5):

x = 0; y = 5; charPoly = (aX - r) (bY - r) - aY bX

Find the Characteristic Roots for the (0,5)-equilibrium

N[Solve[charPoly == 0, r]]

Characteristic roots at (0,0):

The characteristic equation and roots at another equilibrium are:

x = 0; y = 0; charPoly = (aX - r) (bY - r) - aY bXN[Solve[charPoly == 0, r]]

Practice Exercise

Find the characteristic roots associated to each of the equilibria of 
	dx/dt = f = y(1 - x/7 - y/5)
	dy/dt = g = x(1 - x/8 - y/6)

At (8,0):

At (28,-15)

Combined Numerical Computation

The next example shows how to calculate the characteristic roots of all the equilibria of "competition" equations
	
(dx/dt) = x (1 - (x/7) - (6/9) (y/6))(dy/dt) = y (1 - (7/8) (( x)/7) - (y/6))
		
In this case the equilibrium with both coordinates is stable.  This means that the "species"  x  and  y  cooexist stably near this point.

Clear[f, g, x, y, aX, aY, bX, bY]; charPoly = (aX - r) (bY - r) - aY bXf = x (1 - (x/7) - (y/9)); \(g = y\ \((1 - x\/8 - y\/6)\); 
aX\  = \ D[f, \ x]; \ 

    aY\  = \ D[f, \ y]; \ 
bX\  = \ D[g, \ x]; \ 
bY\  = \ D[g, \ y]; \ 

    Print["General partials{{aX,aY},{bX,bY}}"]; 

    \(Print[{{aX, aY}, {bX, bY}}]; \)\)Print["Equilibria at"]; equil = Solve[{f == 0, g == 0}, {x, y}]

At (0,0):

x = 0; y = 0; charPolySolve[charPoly == 0, r]

At (28,9)/5:

x = (28/5); y = (9/5); charPolySolve[charPoly == 0, r]

At (7,0):

x = 7; y = 0; charPolySolve[charPoly == 0, r]

At (0,6):

x = 0; y = 6; charPolySolve[charPoly == 0, r]

All in One Swell Foop

Clear[x, y]charPoly == 0 /. {x, y} -> {x, y} /. equil{{x, y}, Solve[charPoly == 0, r]} /. {x, y} -> {x, y} /. equil

Practice Exercise

Find and classify as attracting, repelling or neither, all equilibria of 
	
			dx/dt = x(1 - x/5 - y/8)
			dy/dt = y(1 - x/6 - y/7)

Clear[f,g,x,y,aX,aY,bX,bY];
charPoly = (aX - r)(bY - r) - aY bX;
f = x(1 - x/5 - y/8);
g = y(1 - x/6 - y/7);
aX = D[f,x];
aY = D[f,y];
bX = D[g,x];
bY = D[g,y];
equil = Solve[{f==0, g==0},{x,y}];
(charPoly==0)/.{x,y}->{x,y}/.equil
{{x,y},Solve[charPoly==0,r]}/.{x,y}->{x,y}/.equil

A Symbolic Computation: Competing Species

The computations throughout this section are done symbolically.  In the Ecology project on competing species you can show geometrically that stability of a positive equilibrium for both of the two "species"  x  and  y  of the equations below depends on the competition terms  b  and  c  being less than one.  (The parameters  h  and  k  are the "carrying capacities" of the environment for the species when they are alone in the environment.)
 	
	dx/dt = m x(1 - x/h - b y/k);
	dy/dt = n y (1 - c x/h - y/k);
			
	The parameter  b  gives a competition multiplier for  y  relative to  y's  carrying capacity in competition with  x  in the  x-equation.  Likewise,  c  gives the relative competitiveness of  x  in the  y-equation.
	
		The stability for this problem can be done as follows, but it is quite a lot easier to solve geometrically, as you can see in the competing species project.

The Positive Equilibrium Point

This computation will find the equilibrium points to be used throughout the computations.

Clear[x, y, h, k, b, c]; equns = {1 == (x/h) + ((b y)/k), 1 == (y/k) + ((c x)/h)}

equil = Solve[equns, {x, y}]

equil = Simplify[equil]

Question

When are both of these points positive? (The parameters are all positive.)

The Matrix of Partial Derivatives

This step will give the matrix of the partial derivatives symbolically.

f = m x(1 - x/h - b y/k);
g = n y(1 - c x/h - y/k);
aX = D[f,x];
aY = D[f,y];
bX = D[g,x];
bY = D[g,y];
matrixOfPartials = {{aX,aY},{bX,bY}};
Print[TableForm[matrixOfPartials]]

The Matrix at the Positive Equilibrium

This step will give the matrix of the partial derivatives at the interior equilibrium point:

x = (h - b h)/(1 - b c);
y = (k - c k)/(1 - b c);
Print[TableForm[Simplify[{{aX,aY},{bX,bY}}]]]

Notice the common factor in the denominator of all these terms

Multiply this matrix by the common factor

TableForm[Simplify[(???){{aX,aY},{bX,bY}}]]

The Characteristic Polynomial at Equilibrium

The next step gives the characteristic equation for the equilibrium.  We can simplify it somewhat by multiplying by the common term in all the denominators.

charPoly = ((aX - r)(bY - r) - aY bX)
simplePoly = Collect[Collect[Together[(b c - 1)charPoly],r]/(b c -1),r]

The General Characteristic Roots

This step will give the characteristic roots for the given equilibrium point which can be used to decide what type of equilibrium point it is.

Solve[charPoly == 0, r]
Solve[simplePoly == 0, r]

When are the characteristic roots opposite signs?

Roots of a quadratic equation

 r^2 + B r + C == 0
 
have the same sign if  4 C  is positive and the discriminant is positive, then Sqrt[B^2 - 4 C] < B.  Our equation in terms of the partial derivative coefficients is

(aX - r)(bY - r) - bX aY = r^2 - (aX + bY) + (aX bY - bX aY) == 0

so in terms of the collected parameters above,

	B = -(aX + bY)
	C = (aX bY - bX aY)

C > 0?

Solve[aX bY == aY bX,{b,c}]

Simplify[(aX bY - aY bX)(1 - b c)^2/(m n)]

When is the product (1 - b)(1 - c)(1 - b c) > 0?
Why is C > 0  or  (aX bY - aY bX) > 0  under the same circumstances?

Solve[(aX + bY)^2 == 4(aX bY - bX aY),{b,c}]
Solve[(m+n)^2 == 4 m n(1-b c),{b,c}]

Notice that the simplified solution gives
	
		b c = -(m - n)^2/(4 m n)
		
and recall that  b  and  c  are positive fractional competition rates.

Exercises

Exercise 24.3.2 a

Find and classify as attracting, repelling or neither, all equilibria of 
	
			dx/dt = x(1 - x/10 - y/5)
			dy/dt = y(1 - x/7 - y/7)

Exercise 24.3.2 b

Find and classify as attracting, repelling or neither, all equilibria of 
	
			dx/dt = x(1 - x/7 - y/7)
			dy/dt = y(1 - x/10 - y/5)

Exercise 24.3.4

Find and classify as attracting, repelling or neither, all equilibria of 
	
			dx/dt = y(1 - x/10 - y/5)
			dy/dt = x(1 - x/7 - y/7)