ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 27 -- Application: Competition Models

Outline of Unit 27

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( student ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

The classical competition model describes an ecosystem with two species that would grow exponentially in isolation from the other species. But, when the two species interact, both populations suffer. If and denote the size of the two species at time , then the model takes the form

> U := [ X(t), Y(t) ]:

> f := (X,Y) -> alpha * X - beta * X*Y:

> g := (X,Y) -> delta * Y - epsilon * X*Y:

> sys1 := op(equate( diff( U, t ), [ f( op(U) ), g( op(U) ) ] ));

> ic1 := X(0) = X[0], Y(0) = Y[0];

>

The unbounded growth of the surviving species in the competition model with exponential growth raises the same concerns as the Malthusian model for a single population. One way to attempt to address this issue is to use logistic models for the competition-free growth terms.

> f := (X,Y) -> alpha * ( 1 - X / M ) * X - beta * X*Y:

> g := (X,Y) -> delta * ( 1 - Y / N ) * Y - epsilon * X*Y:

> sys2 := op(equate( diff( U, t ), [ f( op(U) ), g( op(U) ) ] ));

> ic2 := X(0) = X[0], Y(0) = Y[0];

>

The equilibrium solutions to this model are

> equil := solve( {f(x0,y0)=0,g(x0,y0)=0}, {x0,y0} );

>

There are now four equilibria to consider. Linearization about each equilibrium solution proceeds as usual:

> for X0 in [equil] do

> b := eval( [ f(x0,y0), g(x0,y0) ], X0 );

> A := simplify( eval( evalm(jacobian( [f(x0,y0),g(x0,y0)], [x0,y0] )), X0 ) );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [x0,y0]`=eval([x0,y0],X0) );

> print( A = evalm( A ) );

> print( lambda = [eval(eigenvals( A ))] );

> print( `-------------------------------------------------------------` );

> end do:

>

The equilibrium at the origin is still an unstable node (source). The two new equilibria consist of a finite population of only one species. The additional assumption < ensures that the equilibrium ( 0, ) is a stable node (sink); likewise, ( , 0 ) is a stable node (sink) when < . These two assumptions ensure that < 0 and hence that the fourth equilibrium point is physically possibly ( i.e. , in the first quadrant). The eigenvalues for this case are rather complicated. Determination of the sign of the eigenvalues is more easily accomplished by looking at the trace and determinant of the coefficient matrix for the linearization about this point:

> A = evalm( A );

> trA := simplify( trace( A ) );

> detA := det( A );

> assume(alpha-beta*N<0, delta-epsilon*M<0, alpha>0, delta>0 );

> is( trA < 0 );

> is( detA < 0 );

> unassign( 'alpha', 'beta', 'M', 'delta', 'epsilon', 'N' ):

>

Because the determinant is negative, the eigenvalues must have different signs and the "coexistence" equilibrium is a saddle point for the linear, and nonlinear, systems.

To perform a graphical investigation of the competition model with logistic growth, select numerical values for the parameters

> param2 := { alpha=2, beta=0.5, M=5,

> delta=1, epsilon=0.5,

> N=10, X[0]=5, Y[0]=6 };

>

and obtain a numerical solution to the IVP

> sol2 := dsolve( eval( { sys2, ic2 }, param2 ), U, type=numeric ):

>

A plot of the two components of this system is

> odeplot( sol2, [ [t,X(t)], [t,Y(t)] ], 0..20,

> legend=["Species X", "Species Y"] );

>

This solution exhibits competitive exclusion with the additional property that the surviving species remains bounded.

The direction field shows the general behavior of solutions. The specific initial conditions have been chosen to illustrate how minor changes in the initial conditions can lead to significant differences in the qualitative behavior of solutions.

> ic3 := [ [ X(0)=0.1, Y(0)=0.7 ], [ X(0)=0.1, Y(0)=0.4 ],

> [ X(0)=5, Y(0)=6.3 ], [ X(0)=5, Y(0)=6.6 ] ]:

> DEplot( {eval( sys2, param2 )}, U, t=0..20, X=0..6, Y=0..12,

> ic3, stepsize=0.1 );

>

Note that there is a curve along which the solution approaches the coexistence equilibrium. In practice, however, it is not possible to locate a point on this curve or to expect a numerical solver to be able to track these solutions.

>