ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL
Unit 7 -- Application: The Spruce Budworm
Industrial Mathematics Institute
Columbia, SC 29208
URL: http://www.math.sc.edu/~meade/
E-mail: meade@math.sc.edu
Copyright © 2001 by Douglas B. Meade
All rights reserved
-------------------------------------------------------------------
>
Outline of Unit 7
>
Initialization
> restart;
> with( DEtools ):
> with( plots ):
Warning, the name changecoords has been redefined
> with( linalg ):
Warning, the name adjoint has been redefined
Warning, the protected names norm and trace have been redefined and unprotected
> with( PDEtools ):
>
A classical example of bifurcation in nature is the interaction between the spruce budworm and balsam fir forests in North America. The basic model for the budworm population is a logistic model with a predation term:
> f[logistic] := r*B*(1-B/K);
> f[predation] := beta*B^2/(alpha^2+B^2);
![f[logistic] := r*B*(1-B/K)](unit0724.gif){kind=small}
![f[predation] := beta*B^2/(alpha^2+B^2)](unit0725.gif)
>
where, in the absence of predation,
is the intrinsic growth rate and
is the carrying capacity. All four parameters,
,
,
, and
, are positive. The corresponding differential equation is
> ode := Diff( B, t ) = f[logistic] - f[predation];
>
Via dimensional analysis, the problem can be reduced to one involving only two parameters. Introduce new, dimensionless, dependent and independent variables,
, defined by
> new_var_eq := B(t)=alpha*y(tau),
> t=alpha/beta*tau;
>
and new (dimensionless) parameters
and
defined so that
> new_par_eq := r=beta/alpha*R,
> K=alpha*Q;
>
The corresponding differential equation in terms of the new variables and parameters is
> ode2 := dchange( {new_vars,new_pars}, eval(ode/beta,B=B(t)), [y,tau,Q,R] );
>
which can be simplified to
> ode2 := map( collect, ode2, {R,alpha} );
>
The advantage of this ODE is that it involves only two parameters.
>
To begin the bifurcation analysis, it would be nice to be able to identify all equilibrium solutions for the nondimensionalized ODE
> ode2;
>
While it is easy to see that
is an equilibrium solution, no other equilibrium solution is immediately obvious.
> equil_eq := eval( rhs(ode2), y(tau)=y ) = 0:
> factor( equil_eq );
>
What is apparent is that, in addition to the trivial equilibrium, there can be up to three additional equilibria. While Maple is capable of finding explicit formulas for the remaining equilibria, these expressions are not likely to be terribly useful.
> solve( equil_eq, {y} );
>
The bifurcation analysis presented in
Unit 6, Section B can be used to identify bifurcation points in terms of the parameters
and
. The two necessary conditions are
> bif_eq1 := equil_eq;
> bif_eq2 := diff( eval(rhs(ode2),y(tau)=y), y ) = 0;
>
The parametric solutions to these two equations are
> bif_sol := solve( {bif_eq1,bif_eq2}, {R,Q} );
>
From the fact that the original unknowns and parameters are positive,
,
and
should also be positive. Thus, it is observed that there are physically realistic equilibria only when y>1.
> p1 := plot( eval([Q,R,y=1.001..100],bif_sol), labels=['Q','R'],
> view=[0..200,0..1], numpoints=200 ):
> p2 := textplot([[15,0.1,`Region I`],[50,0.7,`Region I`],
> [100,0.25,`Region II`]]):
> display([p1,p2]);
![[Maple Plot]](unit0762.gif)
>
In Region I there is only one non-trivial equilibria while in Region II there are three non-trivial equilibria. At every point of the boundary between Regions I and II, i.e. , the red curve, there are two non-trivial equilibria. To understand this conclusion, note that the equilibria must satisfy
The function on the left-hand side is linear while the one on the right-hand side does not depend on the parameters. From a graph of these functions it is possible to visually see how the number of equilibria change with the parameter values.
> P := proc(R,Q,T)
> if nargs=2 then T := `` end if;
> plot( [R*(1-y/Q),y/(1+y^2)], y=0..10, title=T )
> end proc:
> display( array([P(0.5,5,`Region I`),P(0.5,7.3,`Boundary`),P(0.5,10,`Region II`)]) );
![[Maple Plot]](unit0764.gif)
>
To see the transition from three, to two, to one equilibrium, the following animation is helpful:
> display( seq( P(0.5, 10-q/10, sprintf("P= 0.5, Q=%4.1f",10.-q/10)), q=0..50), view=0..0.5, insequence=true );
![[Maple Plot]](unit0765.gif)
>
For additional information about this problem -- both biological and mathematical -- please consult the references.
>
1. Ludwig, Jones, and Holling, ``Qualitative Analysis of Insect Breakout Systems: The Spruce Budworm and Forest,'' J. Animal Ecology (1978), 47 , 315-332.
2. Murray, J., Mathematical Biology , Springer-Verlag, 1993.
3. Strogatz, S., Nonlinear Dynamics and Chaos , Addison-Wesley, 1994.
4. McKelvey, S., Spruce Budworm Model, URL: http://www.stolaf.edu/people/mckelvey/envision.dir/spruce.html.
>
[Back to ODE Powertool Table of Contents]
>