Continuing a Saddle-Node Bifurcation of Periodics with AUTO

Joe McKenna

Summary

An example system that undergoes a saddle-node bifurcation of periodics is analyzed. Relations between parameters and variables for which the system undergoes the bifurcation are derived and accurately predict AUTO calculations. An explanation of using AUTO to continue the saddlenode bifurcation of periodics in two parameters is given. The general approach is to locate a Hopf bifurcation on a stationary branch, locate a saddle-node bifurcation of periodics on the periodic branch emanating from the Hopf bifurcation, and continue the saddle-node bifurcation of periodics in two parameters. If the periodic branch or saddle-node bifurcation of periodics is known a priori, the first or first two steps are not necessary to complete the third, resp.

Bifurcation Analysis of Example System

Consider the ODE system in polar coordinates

\{ \begin{matrix}ṙ & = r(r(b − r) − c) for b > 0 \\ θ̇ & = r. \end{matrix}  (1)

It can be seen from the graph of ṙ vs. r that (1) undergoes a saddle-node bifurcation of periodics. An explanation is provided below.

Figure 1: ṙ vs. r for (1).

The stationary radii of (1), given by ṙ = 0, are r_{0} = 0 and the roots that result from translating

r(b − r) down by c, that is, r_{±} = \frac{b}{2} ± \sqrt{\frac{b^{2}}{4} − c}. Note also, since θ̇ ≥ 0, the flow everywhere in the xy-plane except at the origin is counter-clockwise.

For 0 < c < \frac{b^{2}}{4} , there are three distinct stationary radii 0 = r_{0} < r_{−} < r_{+} for which ṙ(r_{0}), ṙ(r_{+}) < 0 and ṙ(r_{−}) > 0 (Fig 1, blue), so in the xy-plane, (1) has a stable equilibrium point (0, 0) surrounded by unstable and stable limit cycles with radii r_{−} < r_{+}, resp.

As c → \frac{b^{2}}{4}^{−}, r_{+} and r_{−} collide at r_{snp} = \frac{b}{2}, for which ṙ(r_{snp}) = 0 and r̈(r_{snp}) < 0 (Fig 1, purple), so in the xy-plane, (1) has a stable equilibrium point (0, 0) surrounded by a semi-stable (attracts for initial conditions with r ≥ r_{snp}, repels otherwise) limit cycle with radius r_{snp}.

For c > \frac{b^{2}}{4} , the only real solution to ṙ = 0 is r_{0} = 0, for which ṙ(r_{0}) < 0 (Fig 1, cyan), so in the xy-plane, (1) has as its only equilibrium the stable equilibrium point (0, 0).

We have shown that a saddle-node bifurcation of periodics occurs at (r_{snp}, c_{snp}) = (\frac{b}{2}, \frac{b^{2}}{4} ). It can also be seen from Figure 1 that (0, 0) becomes stable as it gives rise to an unstable limit cycle when c increases past zero, that is, the system undergoes a subcritical Hopf bifurcation at (r, c) = (0, 0). Below are the system’s phase portraits corresponding to the c values in Figure 1.

Figure 2: Phase portraits corresponding to Figure 1. The x and y range of each graph is 3b and all flows are counter-clockwise.

Now that we have an analytic relation in (b, c) and (c, r = max y) for which the system undergoes a saddle-node bifurcation of periodics, we could graph the bifurcation locus in the bc-plane and the cy-plane by simply plotting c = \frac{b^{2}}{4} and y = \sqrt{c}, resp., but we will use AUTO to generate these graphs.

Figure 3 is the bifurcation diagram of (1) with c as a bifurcation parameter in a neighborhood of zero. We will explain how the data plotted in the diagram was generated using AUTO and how it was used to continue the locus of saddle-node bifurcations of periodics.

Figure 3: Bifurcation diagram.

Using AUTO

For each run, AUTO expects two files, one to specify the system equations and initial conditions and one to specify the constants that dictate the program’s behavior.

The Equations File

For this example, the system equations and initial conditions appear in snp.f90, in the two functions FUNC and STPNT, resp. The values of x, y, b, and c are referred to as U(1), U(2), PAR(1), and PAR(2), resp. The equations (2) appearing in FUNC are derived from (1) as follows. Transforming (1) into rectangular coordinates, we write

( \begin{matrix}x & y \\ −y & x \end{matrix} ) (ẋ ; ẏ) = r^{2} (\begin{matrix}p − c \\ r \end{matrix})

where p(r) = r(b−r), since ṙ = \frac{1}{r} (x y) (ẋ ẏ)^{T} and θ̇ = \frac{1}{r^{2}} (−y x) (ẋ ẏ)^{T} for r = \sqrt{x^{2} + y^{2}}. Multiplying by the inverse of the lefthand side matrix and rearranging, we get

(ẋ ; ẏ) = ( \begin{matrix}x & −y \\ y & x \end{matrix} ) (\begin{matrix}p − c \\ r \end{matrix})
(ẋ ; ẏ) = ( \begin{matrix}p − c & −r \\ r & p − c \end{matrix} ) (\begin{matrix}x \\ y \end{matrix}) .  (2)

The variables and parameters initialized in STPNT are chosen to first locate a Hopf bifurcation on a branch of stationary solutions. Knowing that for b > 0, a Hopf bifurcation occurs as c increases past zero at the equilibrium point (x, y) = (0, 0), we initialize b = 1, c = −0.1, x = 0, and y = 0, and intend to continue the solution (0, 0) with c increasing to 0.5. After locating the Hopf bifurcation, we continue in c the periodic solutions that arise from it to locate a saddlenode bifurcation of periodics. Finally, we continue in b and c the saddle-node bifurcation of periodics.

The Constants Files

For this example, the constants that dictate the program’s behavior appear in three different files c.snp.ss, c.snp.ps, and c.snp.snp (AUTO expects each constants file to begin with the c. prefix) for computing stationary solutions, periodic solutions, and saddle-nodes of periodics, resp. The full meaning of each constant is explained in Chapter 10 of the AUTO documentation (auto.pdf) but important values for the current example are explained here.

In c.snp.ss, the value IPS=1 (problem specification index) indicates that a stationary solution is to be continued and Hopf bifurcations are to be detected. The value IRS=0 (restart index) indicates that the continuation is to start from the initial conditions specified in snp.f90. The value ICP=['c'] (continuation parameter index) indicates that c is the continuation parameter. The value UZSTOP={'c':.5} indicates that continuation should stop when c reaches the value 0.5. Note there are other ways to stop continuation such as specifying a maximum number of steps with NMX, maximum number of bifurcations with MXBF, maximum number of bifurcations of a specific type with SP (for example, SP=['HB1'] would stop the continuation at the first Hopf bifurcation), etc.

In c.snp.ps, the value IPS=2 indicates that a periodic solution is to be continued. The value ILP=1 (limit point index) turns on detection of limit points, which allows us to locate the saddlenode bifurcation of periodics. The value ICP=['c','PERIOD'] indicates that c and period are continuation parameters. By default, AUTO records maximum values of the system variables when continuing periodic solutions. Setting the value IPLT=-2 indicates that the minimum values of y = U(2) are also to be recorded. The value THL={'PERIOD':0.0} neglects period in determining the continuation stepsize to avoid problems near homoclinic orbits, a standard precaution taken when continuing periodic solutions.

In c.snp.snp, the value ICP=['c','b','PERIOD'] indicates that c, b, and period are continuation parameters. The value ISW=2 (branch switching index) allows AUTO to compute a locus of saddle-node bifurcations of periodics.

The AUTO Script

The commands that generate the bifurcation diagram data are summarized in the script snp.auto. They can be run all-at-once from a UNIX command line with auto snp.auto or step-by-step from the AUTO command line with demofile('snp.auto'). Here, we explain the effect of each command in snp.auto.

First, the stationary branch is stored in the bifurcation diagram object ss by running ss=run(e ='snp',c='snp.ss'). The arguments e='snp' and c='snp.ss' specify the names of the equations and constants file, resp. In the terminal output from this command, you should see that AUTO identified a Hopf bifurcation and labelled it 2. Next, the periodic branch is stored in ps by running ps=run(ss('HB1'),c='snp.ps'). The argument ss('HB1') overrides the value of IRS in the constants file c.snp.ps and instructs AUTO to continue from the first Hopf bifurcation located on the stationary branch ss. In the terminal output from this command, you should see that AUTO identified a saddle-node bifurcation of periodics and labelled it 4. Next, starting data for continuing the saddle-node bifurcation of periodics is stored in snpstart by running snpstart=run(ps('LP1'),c='snp.snp'). Last, the locus of saddle-node bifurcations

of periodics is stored in snp by running snp=merge(run(snpstart)+run(snpstart,DS='-')).

This command merges continuations of the saddle-node bifurcation of periodics for both increasing and decreasing values of continuation parameters. The continuation for solely increasing or decreasing values of continuation parameters could be stored in snp by running

snp=run(snpstart,DS='-') or snp=run(snpstart), resp.

The data from each branch are combined and saved in the files b.snp, s.snp, and d.snp in a format explained in the AUTO documentation by running save(ss+ps+snp,'snp'). The temporary files created by AUTO are deleted by running clean.

Plotting the Data

The locus of saddle-node bifurcations of periodics is plotted in Figure 4. The data appear in both b.snp and the last data block of tex/gnuplot/unstable.dat. While using AUTO, any stored bifurcation diagram object, say x, can be graphed with PyPlot by running plot(x) from the AUTO command line. Also, the AUTO documentation describes a few methods for exporting data.

I’ve written a small program cleanData that reformats the data in b.snp into a number of

data blocks in two files tex/gnuplot/stable.dat and tex/gnuplot/unstable.dat. Each data

block corresponds to a connected portion of a branch with the same stability. The source cleanData.cpp may be modified for a different project by recompiling it after changing the name of the input and output data files specfied in the source file.

The formats of stable.dat and unstable.dat are conducive to graphing with gnuplot. For example, the scripts that generated Figures 1 - 4, appearing in plot1.p - plot4.p in tex/gnuplot, use the two data files. These can be run from a UNIX command line with gnuplot plot.p or from the gnuplot command line with load 'plot.p'. Also, the data plotted in Figure 2 were generated by integrating (2) with CVODE by running make (requires CVODE) then snp in

tex/gnuplot/fig2data.

Figure 4: The locus of saddle-node bifurcations of periodics.