Hopf Bifurcation in a Predator-Prey Model

Author: Patri Forwalter-Friedman

Abstract: The behavior of the solutions of a Dynamic System is often strongly dependent upon its parameters. As one varies a parameter continuously, equilibrium points can come and go, spawning limit cycles which then may survive or fade away. An example is Hopf Bifurcation in a predator-prey model. Using animation, we examine the bifurcation as a parameter changes, first with a single trajectory and then with multiple trajectories. Finally, a two-variable animation is created which shows how another parameter in the system affects the bifurcation.

Note: High Memory / CPU requirements!

> with(DEtools):

> with(plots):

> de1 := diff(x(t),t) = (-a + b*y(t)/(c + k*y(t)))*x(t);

> de2 := diff(y(t),t) = (d - e*y(t) - f*x(t)/(c+k*y(t)))*y(t);

Warning, the name changecoords has been redefined

Using the values a=0.5, b = d = e = f = 1, c = 0.3, and letting k be our bifurcation parameter, we have the system:

> bifde1a := subs(a=0.5, b=1, d=1, e=1, f=1, c=0.3, de1);

> bifde2a := subs(a=0.5, b=1, d=1, e=1, f=1, c=0.3, de2);

> bifde1b := (parameter) -> subs(k=parameter, bifde1a);

> bifde2b := (parameter) -> subs(k=parameter, bifde2a);

Notice that our two differential equations in the system are given to Maple as functions of the parameter k . Each DE, when given a value for the parameter, becomes the DE corresponding to that parameter. We will use this feature in conjunction with Maple's powerful tools for animation in order to display our result in a way that lets the viewer see solutions of the the system change dynamically as k varies.

First we create a list of plots with the seq function. Note that each plot has the same initial point, but a different k value:

NOTE:. On most systems, it will take a little while to generate this sequence of plots.

> hopfk := seq(DEplot({bifde1b('i/15'), bifde2b('i/15')}, [x(t), y(t)], 0..140,[[x(0)=0.45,y(0)=0.225]], y=0..2.5, x=0..2.5, arrows=NONE, stepsize=1, linecolour=COLOUR(RGB, 0.5, 0, 0.5)),i=0..25):

Now that the graph has been created, we want to animate it. The use of the subs command allows us to change the line thickness of our resulting plot, as Maple's default size is too thick and makes viewing a plot difficult. The command display , with the insequence=true option, will animate this list of plots.

> hopfk := subs(THICKNESS(3)=THICKNESS(0),[hopfk]):

> display(hopfk, insequence=true);

The single trajectory plot above gives an excellent idea of the behavior of a solution of the system. In the beginning the solution tends to an equilibrium point. As k increases, the equilibrium point moves in the +x , + y direction, and suddenly explodes into a limit cycle, which grows for a while, and then shrinks back towards a point, moving the whole time. Finally, after the cycle turns back into a point, the point moves in the -x, +y direction, settling at a point on the y-axis at a time after our simulation finishes.

To see the behavior of the system more clearly, we may be tempted to use multiple trajectories. Although this lets us see even more information, it tends to be computationally intensive for animated plots, and the screen can get cluttered. If the time and computational power is available, using multiple initial points can be an invaluable tool, especially since it is usually not clear which initial condition one can choose to see the behavior of the entire system. Unlike the above example, such a condition may not even exist. If one can be found, however, it will simplify matters greatly. Below is an example, using the same system of DE's as above, with 16 initial points for each of 17 values of k .

NOTE:. On most systems, it will take a while to generate this sequence of plots.

> hopfk := seq(DEplot({bifde1b('h/10'), bifde2b('h/10')}, [x(t), y(t)], 0..50,[ '['x(0)=i/3','y(0)=j/3'] $i=0..3' $ j=0..3 ], y=0..2.5, x=0..2.5, arrows=NONE, stepsize=0.6, linecolour=COLOUR(RGB, 0.5, 0, 0)),h=0..16):

> display([op(subs(THICKNESS(3)=THICKNESS(0),[hopfk]))], insequence=true);

Now suppose we wish to see how the system's k -dependence changes when we alter one of the other parameters. To do so, we can create an animation with k changing, as above, and then change another parameter, perhaps a , create a new animation, and so on. We will end up with a sequence of animations, each one with k changing over some range, but with a different value for a in each. If we then chain those animations together and display them in sequence, we should be able to see how the parameter affects the system. This is one way of seeing how two parameters affect the behavior of the solutions of a planar system, without using three dimensional projections.

First, we need to make our DE's functions of two parameters. We define two new functions, each of which is passed a value for k and a value for a , and returns the DE corresponding to that value.

> bifde1c := (kparam, aparam) -> subs(k=kparam, a=aparam, b=1, d=1, e=1, f=1, c=0.3, de1);

> bifde2c := (kparam, aparam) -> subs(k=kparam, a=aparam, b=1, d=1, e=1, f=1, c=0.3, de2);

Next we create a function that will make a sequence of plots with k changing, given a value for a .

> hopfa := (aparam) -> seq(DEplot({bifde1c('h/8', aparam), bifde2c('h/8', aparam)}, [x(t), y(t)], 0..120, [[x(0)=0.45,y(0)=0.225]], y=0..2.5, x=0..2.5, arrows=NONE, stepsize=0.9, linecolour=COLOUR(RGB, 0, 0, 0.8)),h=0..18):

Now we create and display a sequence of animations, with a ranging from 0.3 to 0.6, in increments of 0.05.

NOTE:. On most systems, it will take a good long while to generate this sequence of plots.

> hopf_big := hopfa(0.3), hopfa(0.35), hopfa(0.4), hopfa(0.45), hopfa(0.5), hopfa(0.55), hopfa(0.6):

> display(hopf_big, insequence=true);

From this, we can see that the lower the value of a , the more k must increase in order to cause the effects we saw earlier. Also, when a is lower the collapse of the intermediate limit cycle occurs at a higher k value. At low a values, the equilibrium point quickly turns into a large limit cycle, which then shrinks slowly. At high a values, the equilibrium point quickly turns into a small limit cycle, which then shrinks quickly.

This worksheet was created by Patri Forwalter-Friedman at Harvey Mudd College. Special thanks to the Mellon grant for funding this project. The author may be contacted at pforwalt@hmc.edu.