«In this lecture, we will learn how to do a bifurcation analysis with the computer program AUTO. AUTO is built into xppaut, which is where the aut ...»
Bifurcation Analysis with AUTO
Marc R. Roussel
May 20, 2004
In this lecture, we will learn how to do a bifurcation analysis with the computer program
AUTO. AUTO is built into xppaut, which is where the aut part of the name comes from, so we
can just continue to use the tool we started to learn in the last lecture.
We will proceed by example, working with the autocatalator model we started to study last
time. Because AUTO is a tricky piece of software, the instructions given in this set of notes will be much more detailed regarding the operation of the software than would otherwise be the case.
Recall that the xpp input ﬁle is # Autocatalator.ode a’ = mu*(kappa+c) - a*bˆ2 - a b’ = (a*bˆ2 + a - b)/sigma c’ = (b-c)/delta param mu=0.1, kappa=65, delta=2e-2, sigma=5e-3 # The variables range over several orders of magnitude, so # it’s convenient to plot log(a), log(b) and log(c) aux la = log(a) aux lb = log(b) aux lc = log(c) # In order to avoid problems with the logs, # start from a point other than (0,0,0).
a(0) = 1 b(0) = 1 c(0) = 1 # This system is stiff, so we need # an appropriate integrator.
@ METHOD=stiff # The time scale of the oscillations is really fast, # and the spikes are really sharp and high, so we need # to adjust both the integration step size and # the maximum variable value allowed.
@ DT=1e-4, BOUNDS=1e4 @ MAXSTOR=1000000 done In order to use AUTO, we have to start with a known solution, preferably an equilibrium point.
We know from our previous study that there is a stable focus at µ 0¡015. We set µ to this value, ¢ and then get a trajectory. We then hit Initialconds Last to make sure that we have gotten rid ¢ of the transient. If you do a Window/zoom Fit operation, you should see a perfectly ﬂat line.
Better yet, if you look in the data viewer, the values should be perfectly constant. Now click on ¢ File Auto. This will open up the AUTO window.
AUTO needs to be set up. We need to tell it what parameter(s) we will vary, what we will ¢ plot, and so on. Start by clicking on Axes Hi-lo. This will set up AUTO to plot local minima and maxima of one of the components of the computed solutions as a function of a parameter.
A simple limit cycle will, for instance, appear as a pair of points on our diagram. We could in principle choose any of the variables of the model to plot. I choose a.1 Similarly, we could obtain bifurcation diagrams as a function of any of the parameters, but it makes sense to build on our earlier work and to vary µ. We need to set the scale of the bifurcation diagram. We found interesting behavior when µ was varied from 0.015 to 0.154. We didn’t really look at what happens after that. It would therefore make sense to vary µ from our starting point at least up to 0.2. If we look at the graphs we generated last time, we note that log a was never larger than 0. Thus, the maximum in a itself would be less than 1. In the kind bifurcation diagram we’re going to draw now, the variable goes on the Y axis and the parameter on the X axis. We therefore set up the axes
If all went well, the axis labels and limits in the AUTO window should reﬂect the values you just typed in.
We now need to set up the bifurcation diagram computation (as opposed to just the plotting
axes). To do this, click on Numerics and set the following parameters:
1 Note that it’s tempting to use the auxiliary variables representing log a, log b or log c since, as we noted last time, the variables range over several orders of magnitude. Unfortunately, AUTO produces erratic results with auxiliary variables, so it’s best to stick with the model’s basic variables.
2I should call this an Andronov-Hopf bifurcation point, but then AUTO’s abbreviation doesn’t make sense, so in this document I will revert to the older usage.
1. Click on Grab. You can move between the points using your keyboard’s arrow keys. Select the ﬁrst EP and hit the return key. This will reset AUTO’s internal state to what it was when you started. If you get into grab mode and decide you don’t want to grab any points, hitting the escape key will get you out without changing anything.
2. Click on File Reset diagram. This deletes some ﬁles which AUTO uses to keep track of what it’s doing.
3. Click on Clear to erase the contents of the AUTO window.
The ﬁrst two steps are particularly important and must be performed in that order. If you mess up at this point, it will be difﬁcult to recover, and you will have to quit AUTO and start over.
Now click on the Numerics button and increase Nmax to a large value, say 5000. (Note that the value of Nmax when you ﬁrst entered this window was 200, which was the point number of the second EP.) If you click on Run, AUTO will now compute the complete bifurcation diagram. The
result will look at least roughly like this:
Note the nearly vertical section of the branch of periodic solutions, and the corresponding large number of points calculated near µ 1¡539 626 10¡ 2 (points 6 to 208). Near this value of µ, the size of the limit cycle is increasing very rapidly and AUTO must do a lot of work to maintain accuracy. AUTO will sometimes get stuck in regions like this one. The best thing to do then is try to step over the problem by using larger values of Ds and Dsmin. Note also that the line of maxima which runs along the top of the page has counterpart minima along the bottom. They are so close to zero that they disappear in this plot. However, if you manually set a small negative value for ¢ Ymin in Axes Hi-lo, these minima will reappear.
If we try to grab the ﬁrst branch point, we immediately get the MX error status. This indicates that there probably isn’t really a branch point here. AUTO was wrong. Following the second branch point just regenerates points we already have in our diagram. We therefore don’t have to worry about this point. The same is true of the third branch point.
Since there wasn’t any excitement at small µ after all, we will continue to higher values of the parameter. Set Par max back to 0.2 in the Numerics dialog. Grab the rightmost EP on the periodic ¢ branch we have computed. Click on Run Extend. By a similar method, extend the steady-state
branch. Our diagram now has the following appearance:
We can clearly see the stable period-2 cycle emerging from the point where the period-1 cycle becomes unstable.
On the practical side, we’re clearly numbering far too many points, so we’ll increase Npr for the rest of this session. We’ll also increase Nmax, and then extend this branch of solutions. You can click on Abort to stop AUTO if it seems to have ﬁlled in a part of the diagram. We again discover a period-doubling bifurcation at µ 1¡526 725 10¡ 1 which we can continue. We can repeat this procedure until we don’t seem to be getting any more detail. Here’s a piece of our diagram showing
the period-2, period-4 and period-8 orbits:
The period-8 orbit is so close to the period-4 that we can’t actually see it in this picture.
This is an essentially complete diagram since the higher period solutions will be invisible on the scale of this graph.
There are many other types of bifurcations than the simple Andronov-Hopf and period-doubling types seen here. AUTO can also detect torus bifurcations (TR). A torus bifurcation is in some ways analogous to an Andronov-Hopf bifurcation, except that it’s a limit cycle that loses stability, this time to a solution which oscillates around the unstable limit cycle. The result is a trajectory which can be thought of as moving on the surface of a torus (a doughnut).