8.4: Bifurcations in Discrete-Time Models
- Page ID
- 7813
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)
( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\id}{\mathrm{id}}\)
\( \newcommand{\Span}{\mathrm{span}}\)
\( \newcommand{\kernel}{\mathrm{null}\,}\)
\( \newcommand{\range}{\mathrm{range}\,}\)
\( \newcommand{\RealPart}{\mathrm{Re}}\)
\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)
\( \newcommand{\Argument}{\mathrm{Arg}}\)
\( \newcommand{\norm}[1]{\| #1 \|}\)
\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)
\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)
\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)
\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)
\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vectorC}[1]{\textbf{#1}} \)
\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)
\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)
\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)
\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)
\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)
\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)The bifurcations discussed above (saddle-node, transcritical, pitchfork, Hopf) are also possible in discrete-time dynamical systems with one variable:
\[x_{t} =F(x_{t-1}) \label{(8.35)} \]
The Jacobian matrix of this system is, again, a 1×1 matrix whose eigenvalue is its content itself, which is given by \(dF/dx\). Since this is a discrete-time model, the critical condition at which a bifurcation occurs is given by
\[\left|\frac{dF}{dt}\right|_{x=x_{eq}} =1. \label{(8.36)} \]
Let’s work on the following example:
\[x_{t} =x_{t-1} +r -x^{2}_{t-1} \label{(8.37)} \]
This is a discrete-time analog of Eq. (8.1.3). Therefore, it has the same set of equilibrium points:
\[x_{eq}= \pm \sqrt{r} \label{(8.38)} \]
Next, we calculate \(dF/dx\) as follows:
\[\begin{align} \dfrac{dF}{dx} &=(r +x -x^{2})' \\[4pt] &=1-2x \label{(8.39)} \\[4pt] \left|\dfrac{dF}{dt}\right|_{x= \pm \sqrt{r}} &=|1\pm 2\sqrt{r}|\end{align} \label{(8.40)} \]
To find out which value of \(r\) makes this 1, we consider the following four scenarios:
- \(1+2\sqrt{r}=1 \Rightarrow r=0\)
- \(1-2\sqrt{r}=1 \Rightarrow r=0\)
- \(1+2\sqrt{r}=-1 \Rightarrow (no \ real \ solution)\)
- \(1-2\sqrt{r}=-1 \Rightarrow r=1\)
As you see, \(r = 0\) appears as a critical threshold again, at which a saddle-node bifurcation occurs. But now we see another critical threshold, \(r = 1\), also showing up, which was not there in its continuous-time version. This is why I said before that continuous-time models are simpler than discrete-time \ models for bifurcation analysis; discrete-time models can show a new form of bifurcation that wasn’t possible in their continuous-time counterparts. To learn more about this, we can study the stability of the two equilibrium points by checking whether \(|dF/dx|\) is greater or less than 1 at each point. The result is summarized in Table 8.4.1.
Table \(\PageIndex{1}\): Summary of bifurcation analysis of \(x _ { t } = x _ { t - 1 } + r - x _ { t - 1 } ^ { 2 }\).
Equilibrium point | \(r < 0 \) | \(0 < r < 1\) | \(1 < r \) |
---|---|---|---|
\(x_{eq} = √r \) | doesn’t exist | stable | unstable |
\(x_{eq} = -√r \) | doesn’t exist | unstable | unstable |
The result shown in the table is pretty much the same as before up to \(r = 1\), but when \(r > 1\), both equilibrium points become unstable. This must mean that the system is not allowed to converge toward either of them. Then what is going to happen to the system in this parameter range? We can let the system show its actual behavior by numerical simulations. Here is an example:
The results are shown in Fig. 8.4.1, visualizing actual trajectories of the system’s state over time for several values of \(r\). The cobweb plot version is also shown in Fig. 8.4.2 (the code is not shown here; try implementing the code yourself!).
There are some new behaviors that were not seen in continuous-time versions of the same model. For example, at \(r = 0.5\), there was a signature of “overshooting” right before the system state converged to a stable equilibrium point. But more importantly, at \(r = 1\), the system began to fail to converge to a single stable point, and for \(r > 1\), it started oscillating between two distinct states. This is called a \(period-doubling \ bifurcation\). It is a bifurcation that typically occurs in discrete-time systems (but it is also possible in continuous-time systems of higher dimensions), where the system loses stability of a period \(T\) trajectory and begins to move in another trajectory with period \(2T\). The bifurcation observed for \(r = 1\) was a transition from an equilibrium point (which is a periodic trajectory with period \(T = 1\), i.e., the system takes the same state value in each time step)
to a trajectory with period \(2T = 2\). Interestingly, the period-doubling bifurcations happen in a cascade if you keep increasing \(r\). In this particular case, another period-doubling bifurcation was observed at \(r = 1.5\), where the period-2 trajectory lost its stability and the system moved to a period-4 trajectory. This continues as you continue to increase r. The first period-doubling bifurcation from period-1 to period-2 trajectories can still be characterized as the loss of stability of an equilibrium point, but the dominant eigenvalue destabilizing the equilibrium point must be \(negative\) in order to induce the flipping behavior. This can be mathematically written as follows:
A first period-doubling bifurcation from period-1 to period-2 trajectories occurs in a discrete-time model when the eigenvalues \(λ_i\) of the Jacobian matrix at an equilibrium point satisfy the following:
\(λ_i = −1\) for some \(i\), while\(|λ_i| < 1\) for the rest.
In the example above, \(dF/dx|_{x_{eq}}=√r = −1\) when \(r = 1\), which triggers the first period doubling bifurcation.
Consider the following discrete-time dynamical system:
\[x_{t} =(1-a)x_{t-1}+ax^{3}_{t-1} \label{(8.41)} \]
This equation has \(x_{eq} = 0\) as an equilibrium point. Obtain the value of \(a\) at which this equilibrium point undergoes a first period-doubling bifurcation.
Once the system begins to show period-doubling bifurcations, its asymptotic states are no longer captured by the locations of analytically obtained equilibrium points, as drawn in bifurcation diagrams (e.g., Figs. 8.2.1, 8.2.3, etc.). However, there is still a way to visualize bifurcation diagrams numerically by simulating the behavior of the system explicitly and then collecting the actual states the system visits for a certain period of time. Then we can plot their distributions in a diagram. The data points should be collected after a sufficiently long initial transient time has passed in each simulation, so that the system’s trajectory is already showing its “final” behavior. Here is a sample code showing how to draw such a bifurcation diagram numerically:
In this code, \(r\) is gradually varied from 0 to 2 at intervals of 0.01. For each value of \(r\), the model (Equation \ref{(8.37)}) is simulated for 200 steps, and only the second half of the state values are recorded in result
. Once the simulation is finished, the states stored in result
are plotted at a particular value of r
in the plot (note that the expression [r] * 100
in Python produces a list of one hundred r
’s). The alpha option is used to make the markers transparent so that the densities of markers are also visible.
The result is shown in Fig. 8.4.3. This bifurcation diagram still shows how the system’s state depends on the parameter, but what is plotted here are no longer analytically obtained equilibrium points but numerically sampled sets of asymptotic states, i.e., where the system is likely to be after a sufficiently long period of time. If the system is converging to a stable equilibrium point, you see one curve in this diagram too. Unstable equilibrium points never show up because they are never realized in numerical simulations. Once the system undergoes period-doubling bifurcations, the states in a periodic trajectory are all sampled during the sampling period, so they all appear in this diagram. The period doubling bifurcations are visually seen as the branching points of those curves.
But what are those noisy, crowded, random-looking parts at \(r > 1.7\)? That is chaos, the hallmark of nonlinear dynamics. We will discuss what is going on there in the next chapter.
Period-doubling bifurcations and chaos are not just for abstract, contrived mathematical equations, but they can occur in various models of real-world biological, ecological, social, and engineering phenomena. The simplest possible example would be the logistic map we introduced in Section 5.5:
\[x_{t} =rx_{t-1}(1-x_{t-1}) \label{(8.42)} \]
This is a mathematical model of population dynamics, where \(x_t\) represents the population of a species that reproduce in discrete (non-overlapping) generations. This model was used by British/Australian mathematical ecologist Robert May in his influential 1976 Nature paper [32] to illustrate how a very simple mathematical model could produce astonishingly complex behaviors.
- Conduct a bifurcation analysis of this model to find the critical thresholds of \(r\) at which bifurcations occur.
- Study the stability of each equilibrium point in each parameter range and summarize the results in a table.
- Simulate the model with several selected values of \(r\) to confirm the results of analysis.
- Draw a bifurcation diagram of this model for \(0 < r < 4\).
The stability of a period-2 trajectory of a discrete-time model \(x_t = F(x_{t−1})\) can be studied by the stability analysis of another model made of a composite function of \(F\):
\[y_{r} =G(y_{\tau -1})=F(F(t_{\tau -1})) \label{(8.43)} \]
This is because the period-2 trajectory in \(F\) corresponds to one of the equilibrium points of \(G(◦) = F(F(◦))\). If such an equilibrium point of \(G\) is being destabilized so that \(dG/dx ≈−1\), it means that the period-2 trajectory in question is losing the stability, and thus another period-doubling bifurcation into a period-4 trajectory is about to occur. Using this technique, analytically obtain the critical threshold of \(r\) in Equation \ref{(8.37)} at which the second period-doubling bifurcation occurs (from period 2 to period 4). Then, draw cobweb plots of \(y_τ = G(y_{τ−1})\) for several values of r near the critical threshold to see what is happening there.