In this lab we will use a program called "xppaut" to examine the bifurcations of a planar system. The program "xppaut" is specialized to solve and analyze the dynamics of ordinary differential equations and maps. The program was written by Prof. B Ermentrout (U. Pitt, Math Dept) and, unlike other software, it is free. It comes with a full documentation and a web based tutorial. "xppaut" is pronounced "X-P-P-AUT" which is a shortened acronym for "X window based program for Phase Planes and AUTo". The later subprogram "AUTO" creates bifurcation diagrams using numerical continuation algorithms (psuedoarclength continuation). It would take a long time to provide the background to explain these algorithms in detail. Roughly, for fixed points, a root finder is needed to estimate the rate of change of the fixed point in the parameter(s). Recall that if the Jacobian is nonsingular, the "direction" fixed points move could be estimated. AUTO keeps track of stability information by numerical algorithms which compute eigenvalues and eigenvectors. The detection of Hopf points and the ensuing continuation of periodic orbits involves solving certain boundary value problems. This is especially involved and we will not discuss it here.
Some older PDF documentation .
To use xppaut, the user must first create a file with a *.ode suffix which defines the problem. If the code has the name "file.ode", one executes the following UNIX command in a terminal window:
xppaut file.ode &
We are going to use xppaut to examine the model described below:
The code for this problem looks like:
Lines starting with # are comment lines.
Lines starting with "p" are parameter definition lines
x' is the derivative of x(t) in t
The last line must be "d" for "done"
Execute the code in a terminal window:
xppaut Burst_FS.ode &
You should have something like:
and a bunch of other windows. Close all but the one you see above.
First let's get the window set up for phase portraits. Right now it is set up plotting solutions as functions of time t. To change it:
(1) Select "Viewaxes"
(2) Select "2D"
Change the X-axis and Y-axis in this menu to look like the following, then select OK:
From the main xppaut (black) window menu:
(3) Select "Nullcline"
(4) Select "New"
These are the x and y nullclines for the current parameter values. Their intersection locates fixed points (equilibria or steady states).
(5) Select "Initialconds"
(6) Select "(M)ouse"
(7) Click your mouse on an initial condition location
(8) Repeat (5)-(6)
You should have something like the following showing lots of trajectories.
Note there appears to be a stable spiral at the nullcline intersection. xppaut can do numerical linearization about fixed points. This is done thru the "Sing pts" menu (some people refer to fixed points as singular points...not very common terminology).
(9) Select "Sing pts"
(10) Select "(M)ouse"
(11) Point you mouse arrow close to the fixed point
(12) Answer "yes"
The resulting "Equilibria" window then summarizes stability, location and eigenvalue information. Here "c- =2" indicating that there are two complex eigenvalues with a negative real part.
You can then use the "Del" key next to the "BackSpace" key to delete old parameter values. When you're happy with the values, select "Ok".
Alternately one can set up one of the three "sliders" in white at the bottom.
(13) Click on any "Par/Var?" slider at the bottom
(14) Change the contents of the window to match the following:
Now the parameter "mu" can be changed by moving the little black slider bar between the ranges (-1,6) that were just set up.
Exercise 1: Flow
With mu=-1 create a phase portrait: (15) "Erase"
(16) "Nullcline",..."New" (note the eminant Saddle-Node!!)
(17) "Dir.field/flow",..."Flow" and make sure the "Grid" (at top left) is set to 5 (10 takes too long). Then "Return". Should get:
The last menu item is used to determine flow or superimpose vector field. I don't really recommend using it too much since it takes alot of time and is sometimes so crowded that the result is un-informative. You may want to just used "Initialconds" then "Mouse" like we did earlier.
Exercise 2: Saddle-Node
Repeat (15)-(16) above for mu values between -1 and about 1.5. You should see the x and y nullcines intersect in a Saddle-Node bifurcation near mu=1.
Exercise 3: Stable-Unstable Manifolds
Set mu=1.52 (or there abouts) using the slider and repeat (15)-(16) to get the new nullclines. There are now three fixed points all along the x axis.
(18) Using the "Sing Pts" menu and then "Menu" determine the stability of the MIDDLE fixed point. When it asks you to draw the "Invariant Sets" answer yes. There are 4. Between each you need to hit the "Esc" key and sometimes "Ok" if the trajectory goes out of bounds. The "Invariant sets" in this case are the stable and unstable manifolds associated with the saddle. They are color coded. Can you determine which color is the stable manifold?
Exercise 4: Homoclinic Bifurcations
Repeat Exercise 3 for values in the range (3,4). Somewhere near mu=3.5 the stable and unstable manifolds intersect. When this happens as a parameter is varied it is called a Homoclinic Bifurcation. Try to track the periodic orbit as you do this for different parameter values.
In order to generate a bifurcation diagram the AUTO part of xppaut must be initiated with a known fixed point from which to continue. In our case, we can verify that when mu=-15 a fixed point is at (x,y)=(3,0). In fact we just used xppaut to compute this above. To initiate AUTO do the following:
(19) Select "Initial conds"
(20) Select "New"
(21) At the top of the main (black) window enter the x coordinate of the fixed point x=3, then return. Then the y=0 value, "return".
Now we're ready.
(22) Select "File"
(23) Select "Auto"
A new bifurcation diagram window will appear. It looks like:
We will create a diagram of "x" versus "mu". We need to stipulate the range of for the window appearance and of the "run" separately.
(24) Select "Axes" under the Auto Menu. Then select "hI-low"
(25) Change the numbers so that the range of "mu" on the "X" axis is (-16,6), i.e., the Xmin=-16, Xmax=6. The range of "x" in the differential equations is actually the range of "Y" on this diagram. Set Ymin=-3,Ymax=4.
(26) We need to change some of the numerical parameters. Select "Numerics" and then change the parameters to match those below:
Here we changed the parameters:
Nmax = 500 (maximum number of steps)
Npr = 500 (print out information every 500 steps)
Dsmax=0.1 (maximum step size in continuation algorithm)
Par Min = -15 (min value of mu for run)
Par Min = 5.5 (max value of mu for run)
Note that Par Min and Par Max do not need to match the window size we changed in items (24)-(25) above.
Here goes:
(27) Select "Run"
(28) Select "Steady State"
These are the fixed points. The THICK line segments indicate STABLE fixed points. Thin line segments indicate unstable fixed points.
Note that there are two Saddle-Node bifurcation points (labelled 3 and 4 in the diagram). At the point "2" in the diagram I showed in class there is a Hopf bifurcation. xppaut knows this from the calculations it just did.
Now lets "continue" from the Hopf bifurcation. What I mean here is we are going to have the AUTO part of xppaut compute the periodic solutions eminating from this Hopf point.
(29) Select "Grab"
(30) Hit the "Tab" key once
Notice at the bottom "HB" is listed under the "Ty" column. This means that AUTO has detected a Hopf Bifurcation. Other important listings are that mu=0.999 at the HB and the period of orbits there are 2.094 (both numerical approximations).
(31) Hit the "return" key. You have now told AUTO you wish to continue from that labelled point, i.e., the one you just "grabbed".
(32) Select "Run"
(33) Select "Periodic"
You should now have a complete bifurcation diagram showing the branch of (stable) periodic orbits (dark dots) eminating from the HB and terminating on the middle branch of fixed points (which in class I showed were saddles). It should look like:
(36) "File"
(37) "Quit"