Application Center - Maplesoft

App Preview:

Lesson 7: Application: The Spruce Budworm

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Lesson07.mw

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 7 -- Application: The Spruce Budworm

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

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 Lesson 7

7.A Biological Background

7.B Bifurcation Analysis

7.C References

>

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( PDEtools ):

Warning, the name changecoords has been redefined

>

7.A Biological Background

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)

f[predation] := beta*B^2/(alpha^2+B^2)

>

where, in the absence of predation, r is the intrinsic growth rate and K is the carrying capacity. All four parameters, r , K , alpha , and beta , are positive. The corresponding differential equation is

> ode := Diff( B, t ) = f[logistic] - f[predation];

ode := Diff(B, t) = r*B*(1-B/K)-beta*B^2/(alpha^2+B^2)

>

Via dimensional analysis, the problem can be reduced to one involving only two parameters. Introduce new, dimensionless, dependent and independent variables, y = y(tau) , defined by

> new_var_eq := B(t)=alpha*y(tau),

>              t=alpha/beta*tau;

new_var_eq := B(t) = alpha*y(tau), t = alpha*tau/beta

>

and new (dimensionless) parameters R and Q defined so that

> new_par_eq := r=beta/alpha*R,

>              K=alpha*Q;

new_par_eq := r = beta*R/alpha, K = alpha*Q

>

The corresponding differential equation in terms of the new variables and parameters is

> ode2 := dchange( {new_var_eq, new_par_eq}, subs(B=B(t), ode/beta), [y,tau,Q,R] );

ode2 := diff(y(tau), tau) = (beta*R*y(tau)*(1-y(tau)/Q)-beta*alpha^2*y(tau)^2/(alpha^2+alpha^2*y(tau)^2))/beta

>

which can be simplified to

> ode2 := map( collect, ode2, {R,alpha} );

ode2 := diff(y(tau), tau) = y(tau)*(1-y(tau)/Q)*R-y(tau)^2/(1+y(tau)^2)

>

The advantage of this ODE is that it involves only two parameters.

>

7.B Bifurcation Analysis

To begin the bifurcation analysis, it would be nice to be able to identify all equilibrium solutions for the dimensionless ODE

> ode2;

diff(y(tau), tau) = y(tau)*(1-y(tau)/Q)*R-y(tau)^2/(1+y(tau)^2)

>

While it is easy to see that y = 0 is an equilibrium solution, no other equilibrium solution is immediately obvious.

> equil_eq := eval( rhs(ode2), y(tau)=y ) = 0:

> factor( equil_eq );

-y*(-R*Q-R*y^2*Q+R*y+R*y^3+y*Q)/(Q*(1+y^2)) = 0

>

It is apparent from the numerator 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 so cumbersome that they are not likely to be very useful.

> solve( equil_eq, {y} );

{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...{y = 0}, {y = 1/6*((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2)*((4*R^3+12*R^2*Q+8*Q^2*R^3+12*Q^2*R-20*Q^3*R^2+4*Q^4*R^3+4*Q^3-Q^4*R)/R)^(1/2))*R^2)^(1/3)/R+2/3*(-3*R-3*Q+Q^2*R)/((72*R*Q-36*Q^2+8*Q^3*R+12*3^(1/2...

>

However, they do serve to show that the equation

> equil_eq;

y*(1-y/Q)*R-y^2/(1+y^2) = 0

>

defines y = y(Q, R) either implicitly, or explicitly.  All equilibrium values y are a function of the two parameters Q and R , and therefore lie on a surface above the QR -plane.  This surface is seen in Figure 7.1.

> implicitplot3d(equil_eq, Q=0..100, R=0..1, y=0..10, axes=box, grid=[20,20,20], title="Figure 7.1");

[Plot]

>

Above a point in the QR -plane where this surface folds over itself, there are multiple equilibria.  Careful inspection of this surface shows there are regions in the QR -plane where only one equlibrium solution exists, and regions where there are clearly three equilibrium solutions.  Exactly along the leading edge of a fold in the surface, there will be just two equilibrium solutions.  It is along such folds that the bifurcation points are found.

The bifurcation analysis presented in Lesson 6, Section B can be used to identify bifurcation points in terms of the parameters Q and R . The two necessary conditions are

> bif_eq1 := equil_eq;

> bif_eq2 := diff( eval(rhs(ode2),y(tau)=y), y ) = 0;

bif_eq1 := y*(1-y/Q)*R-y^2/(1+y^2) = 0

bif_eq2 := (1-y/Q)*R-y*R/Q-2*y/(1+y^2)+2*y^3/(1+y^2)^2 = 0

>

The parametric solutions to these two equations are

> bif_sol := solve( {bif_eq1,bif_eq2}, {R,Q} );

bif_sol := {Q = 2*y^3/(y^2-1), R = 2*y^3/(1+2*y^2+y^4)}

>

Since the original unknowns and parameters are positive, y , Q and R should also be positive. Thus, there are physically realistic equilibria only when y*`>` .

Figure 7.2 shows the regions I and II  in the first quadrant of the QR -plane where there will be either one or three nontrivial equilibrium solutions, respectively.

> 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], title="Figure 7.2");

[Plot]

>

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

R*(1-y/Q) = y/(1+y^2)

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 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`)]) );

[Plot]

>

To see the transition from three, to two, to one equilibrium, the animation in Figure 7.3 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, title="Figure 7.3" );

[Plot]

>

For additional information about this problem -- both biological and mathematical -- please consult the references.

>

7.C 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 ]

>