#LyX 1.3 created this file. For more info see http://www.lyx.org/
\lyxformat 221
\textclass article
\language english
\inputencoding auto
\fontscheme default
\graphics default
\paperfontsize default
\papersize Default
\paperpackage a4
\use_geometry 0
\use_amsmath 0
\use_natbib 0
\use_numerical_citations 0
\paperorientation portrait
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\defskip medskip
\quotes_language english
\quotes_times 2
\papercolumns 1
\papersides 1
\paperpagestyle default
\layout Title
Stabilizing a Subcritical Bifurcation in a Mapping Model of CardiacMembrane
Dynamics
\layout Author
Matthew A.
Fischer
\begin_inset Foot
collapsed true
\layout Standard
Department of Mathematics, Duke University, Durham, North Carolina 27708,
USA
\end_inset
and Colin B.
Middleton
\begin_inset Formula $^{*}$
\end_inset
\layout Date
April 25th 2006
\begin_inset Foot
collapsed true
\layout Standard
This work was supported by the NSF VIGRE grant NSFDMS9983320
\end_inset
\layout Section
Introduction
\layout Subsection
Background Information
\layout Standard
Under normal circumstances, the heart responds identically to each stimulus
in a periodic train of stimuli.
By contrast, the term
\emph on
alternans
\emph default
refers to a phaselocked periodtwo response in which alternate stimuli
evoke one of two different responses.
It is believed that alternans may be a precursor to ventricular fibrillation,
an arrhythmia causing sudden cardiac death.
\layout Standard
Several authors have studied the use of feedback to suppress alternans
\begin_inset LatexCommand \cite{key2}
\end_inset
.
Indeed, this has been successfully implemented in small pieces of paced
\emph on
in vitro
\emph default
bullfrog heart
\begin_inset LatexCommand \cite{key2}
\end_inset
.
In this paper we study theoretically a variant of such control problems.
\layout Standard
A full description of the electrical response of the heart to stimulation
is difficult, involving a very stiff, manyvariable coupled system of timedepe
ndent PDE (of reactiondiffusion type) in an exceedingly complicated 3D
geometry.
If the spatial dependence in the problem is suppressed (which is approximately
valid for small pieces of
\emph on
in vitro
\emph default
heart), the description reduces to a system of ODE, often called an
\emph on
ionic model
\emph default
since what is being described is the transport of ions across the cell
membrane.
The most important variable in this system (and the most accessible experimenta
lly) is the transmembrane potential.
Typical responses of this variable to a periodic sequence of stimuli are
shown in Figure 1.
\begin_inset Float figure
wide false
collapsed false
\layout Standard
\align center
\begin_inset Graphics
filename C:/mstuff/images/fig1alt.bmp
lyxscale 60
scale 60
keepAspectRatio
\end_inset
\layout Caption
Transmembrane potential response to periodic stimulation.
In part a), the stimulation is such that 1:1 response is observed.
At more rapid stimulation as in part b), 2:2 response or alternans is observed.
\end_inset
This figure illustrates the property of cardiac tissue that is called
\emph on
excitability
\emph default
: i.e., a small stimulus leads to a raised voltage for an extended period
of time, far longer than the duration of the stimulus itself, after which
the voltage returns to its resting value.
Such a time course of the voltage is called an
\emph on
action potential
\emph default
.
Figure 1a shows a steady situation in which the same response occurs to
every stimulus in a periodic sequence, what is called
\emph on
\emph default
1:1 response.
Such 1:1 response is typical when the stimulus period is large.
By contrast, Figure 1b shows an alternans or 2:2, response.
\layout Standard
To understand arrhythmias, one needs to examine the dynamics of stimulation
and response.
Much insight into the dynamics can be gained from a simple model introduced
by Nolasco and Dahlen
\begin_inset LatexCommand \cite{key4}
\end_inset
in 1968.
In this approximation, each action potential is characterized by a single
number, its duration (acronym: APD).
More formally, the APD denotes the length of the time interval during which
the voltage exceeds a critical voltage.
Nolasco and Dahlen proposed the approximation: there is a function
\begin_inset Formula $F$
\end_inset
such that the evolution of successive APD's under repeated stimulation
is given by
\layout Standard
\align center
\begin_inset Formula \begin{equation}
APD_{n+1}=F(DI_{n})\label{eq:1}\end{equation}
\end_inset
\layout Standard
\align left
where
\begin_inset Formula $DI_{n}$
\end_inset
, an acronym for
\emph on
diastolic interval
\emph default
, equals the time elapsed between the end of the
\begin_inset Formula $n^{th}$
\end_inset
action potential and the next stimulus.
If the stimuli are applied with period
\emph on
BCL
\emph default
(acronym for basic cycle length), then
\begin_inset Formula $DI_{n}=BCLAPD_{n}$
\end_inset
, so
\begin_inset Formula $APD_{n}$
\end_inset
evolves according to the iterated mapping
\layout Standard
\align center
\begin_inset Formula \begin{equation}
APD_{n+1}=F(BCLAPD_{n})\label{eq:2}\end{equation}
\end_inset
\layout Standard
Based on experiment, previous authors proposed the function
\begin_inset LatexCommand \cite{key4}
\end_inset
\layout Standard
\align center
\begin_inset Formula \begin{equation}
F(DI)=A_{max}Ce^{DI/\tau},\label{eq:3}\end{equation}
\end_inset
\layout Standard
\align left
where
\begin_inset Formula $A_{max}$
\end_inset
,
\begin_inset Formula $C$
\end_inset
and
\begin_inset Formula $\tau$
\end_inset
are constants.
Various other forms have also been proposed.
In particular, it was shown in
\begin_inset LatexCommand \cite{key5}
\end_inset
, using asymptotics, that the behavior of the ODE in a simple ionic model
can be approximately described by (
\begin_inset LatexCommand \ref{eq:1}
\end_inset
) with
\layout Standard
\align center
\begin_inset Formula \begin{equation}
F(DI)=\tau_{close}\ln\{\frac{1(1h_{min})e^{DI/\tau_{open}}}{h_{min}}\}\label{eq:4}\end{equation}
\end_inset
\layout Standard
\align left
where
\begin_inset Formula $\tau_{close}$
\end_inset
,
\begin_inset Formula $h_{min}$
\end_inset
, and
\begin_inset Formula $\tau_{open}$
\end_inset
are parameters derived from the ionic model.
\layout Standard
Despite the oversimplification of the approximation (
\begin_inset LatexCommand \ref{eq:2}
\end_inset
), it can successfully model alternans.
In explaining this, we assume the reader is familiar with the elementary
theory of 1D iterated mappings, as discussed for example in Chapter 10
of
\begin_inset LatexCommand \cite{key6}
\end_inset
.
Suppose that the function
\begin_inset Formula $F(DI)$
\end_inset
is monotone increasing.
(Both (
\begin_inset LatexCommand \ref{eq:3}
\end_inset
) and (
\begin_inset LatexCommand \ref{eq:4}
\end_inset
) have this property.) Then, as illustrated in Figure 2, for any
\begin_inset Formula $BCL$
\end_inset
, there is a unique fixed point
\begin_inset Formula $APD_{*}$
\end_inset
such that
\begin_inset Float figure
wide false
collapsed false
\layout Standard
\align center
\begin_inset Graphics
filename C:/mstuff/images/fig2alt.bmp
scale 70
keepAspectRatio
\end_inset
\layout Caption
Graph of
\begin_inset Formula $F(BCLAPD_{n})$
\end_inset
and the line
\begin_inset Formula $APD_{n+1}=APD_{n}$
\end_inset
.
It can be seen that the sole fixed point
\begin_inset Formula $APD_{*}$
\end_inset
of
\begin_inset Formula $F(BCLAPD_{n})$
\end_inset
is at its intersection with the line
\begin_inset Formula $APD_{n+1}=APD_{n}$
\end_inset
.
\end_inset
\layout Standard
\align center
\begin_inset Formula \begin{equation}
APD_{*}=F(BCLAPD_{*})\label{eq:5}\end{equation}
\end_inset
\layout Standard
If such a fixed point is stable, then
\begin_inset Formula $APD_{n}\rightarrow APD_{*}$
\end_inset
is a possible asymptotic behavior for
\begin_inset Formula $n\rightarrow\infty$
\end_inset
of a sequence generated by (
\begin_inset LatexCommand \ref{eq:2}
\end_inset
).
If
\begin_inset Formula $BCL$
\end_inset
is large, then the derivative
\begin_inset Formula $F^{\prime}(BCLAPD_{*})$
\end_inset
at the fixed point is small, in particular less than unity, so the fixed
point is indeed stable
\begin_inset LatexCommand \cite{key6}
\end_inset
.
However, as
\begin_inset Formula $BCL$
\end_inset
decreases,
\begin_inset Formula $F^{\prime}(BCLAPD_{*})$
\end_inset
increases, and if
\begin_inset Formula $F^{\prime}(BCLAPD_{*})$
\end_inset
passes through unity, the iteration suffers a perioddoubling bifurcation.
\layout Standard
This behavior is conveniently summarized in a bifurcation diagram (see Figure
3), which for each
\begin_inset Formula $BCL$
\end_inset
, plots the limit points of sequences
\begin_inset Formula $\{ APD_{n}\}$
\end_inset
generated by (
\begin_inset LatexCommand \ref{eq:2}
\end_inset
).
Figure 3a illustrates a supercritical bifurcation, the more familiar alternativ
e and the only possibility for the mapping (
\begin_inset LatexCommand \ref{eq:3}
\end_inset
)
\begin_inset Float figure
wide false
collapsed false
\layout Standard
\align center
\begin_inset Graphics
filename C:/mstuff/images/fig3alt.bmp
lyxscale 60
scale 60
keepAspectRatio
\end_inset
\layout Caption
Bifurcation diagram showing
\begin_inset Formula $APD_{*}$
\end_inset
versus
\begin_inset Formula $BCL$
\end_inset
for a supercritical bifurcation (a) and a subcritical bifurcation (b).
Unstable critical points are shown in dashed lines.
\end_inset
.
However, it was shown in
\begin_inset LatexCommand \cite{key5}
\end_inset
that, for certain parameter ranges, the bifurcation of (
\begin_inset LatexCommand \ref{eq:4}
\end_inset
) is subcritical, as illustrated in Figure 3b.
\layout Standard
In this paper we are concerned with the case of a subcritical bifurcation.
Note that there is a range
\begin_inset Formula $BCL_{bif}1$
\end_inset
and then becomes stable via a subcritical bifurcation into points on branches
2 and 3 (periodtwo points) at
\begin_inset Formula $\mu=1$
\end_inset
.
These periodtwo points are unstable throughout their range.
\layout Itemize
As the fixed point on branch 1 (fixed point
\begin_inset Formula $x=0$
\end_inset
) loses stability at
\begin_inset Formula $\mu=1$
\end_inset
, it bifurcates into points on branches 4 and 5 (fixed points).
These fixed points begin stable and lose their stability at
\begin_inset Formula $\mu=2.$
\end_inset
\layout Itemize
At
\begin_inset Formula $\mu=2$
\end_inset
points on branches 4 and 5 (fixed points) lose stability and each bifurcates
into a pair of stable periodtwo points.
The fixed point on branch 4 bifurcates into periodtwo points on branches
6 and 7 while the fixed point on branch 5 bifurcates into periodtwo points
on branches 8 and 9.
The stable periodtwo points lose their stability at
\begin_inset Formula $\mu=\sqrt{5}$
\end_inset
.
\layout Itemize
At
\begin_inset Formula $\mu=\sqrt{5}$
\end_inset
the periodtwo points on branches 6 & 7 and 8 & 9 each lose stability via
a perioddoubling bifurcation.
Beyond
\begin_inset Formula $\mu=\sqrt{5}$
\end_inset
, perioddoubling bifurcations continue as
\begin_inset Formula $\mu$
\end_inset
becomes more negative.
\layout Standard
This map has three fixed points and six periodtwo points.
The following table is a summary of their behavior:
\layout Standard
\align center
\begin_inset Tabular
\begin_inset Text
\layout Standard
Point
\end_inset

\begin_inset Text
\layout Standard
Type
\end_inset

\begin_inset Text
\layout Standard
Exists for
\end_inset

\begin_inset Text
\layout Standard
Stable for
\end_inset

\begin_inset Text
\layout Standard
Maps to
\end_inset

\begin_inset Text
\layout Standard
1
\end_inset

\begin_inset Text
\layout Standard
Fixed point
\end_inset

\begin_inset Text
\layout Standard
all
\begin_inset Formula $\mu$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
1
\end_inset

\begin_inset Text
\layout Standard
2
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
Never
\end_inset

\begin_inset Text
\layout Standard
3
\end_inset

\begin_inset Text
\layout Standard
3
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
Never
\end_inset

\begin_inset Text
\layout Standard
2
\end_inset

\begin_inset Text
\layout Standard
4
\end_inset

\begin_inset Text
\layout Standard
Fixed point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $2<\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
4
\end_inset

\begin_inset Text
\layout Standard
5
\end_inset

\begin_inset Text
\layout Standard
Fixed point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $2<\mu<1$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
5
\end_inset

\begin_inset Text
\layout Standard
6
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\sqrt{5}<\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
7
\end_inset

\begin_inset Text
\layout Standard
7
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\sqrt{5}<\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
6
\end_inset

\begin_inset Text
\layout Standard
8
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\sqrt{5}<\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
9
\end_inset

\begin_inset Text
\layout Standard
9
\end_inset

\begin_inset Text
\layout Standard
Periodtwo point
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
\begin_inset Formula $\sqrt{5}<\mu<2$
\end_inset
\end_inset

\begin_inset Text
\layout Standard
8
\end_inset

\end_inset
\layout Subsection
Feedback
\layout Standard
Recall that in Hall and Gauthier
\begin_inset LatexCommand \cite{key2}
\end_inset
feedback is used to suppress alternans.
In this paper we analyze an analogous problem: stabilizing a subcritical
bifurcation.
We would like to be able to devise a method using feedback to observe alternans
when theoretically possible but not observable under normal circumstances.
Such a method could be used to test predictions of models, some of which
predict subcritical bifurcations
\begin_inset LatexCommand \cite{key5}
\end_inset
.
\layout Standard
To stabilize the periodtwo points, we wish to perturb iterates of (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
) by applying feedback, where an example of applying feedback is below in
(
\begin_inset LatexCommand \ref{eq:7}
\end_inset
).
In Hall and Gauthier
\begin_inset LatexCommand \cite{key2}
\end_inset
, feedback is an adjustment of the BCL, which is an argument of (
\begin_inset LatexCommand \ref{eq:4}
\end_inset
).
This is an indirect method of adjusting the
\begin_inset Formula $n+1st$
\end_inset
iterate.
In our analysis, we will adjust the
\begin_inset Formula $n+1st$
\end_inset
iterate directly in the following manner:
\layout Standard
\begin_inset Formula \begin{equation}
x_{n+1}=f(x_{n})+\gamma g(x_{n},x_{n1},...,x_{nk})\label{eq:7}\end{equation}
\end_inset
\layout Standard
\align left
Where
\begin_inset Formula $g(x_{n},x_{n1},...,x_{nk})$
\end_inset
is a corrective function that adjusts up or down the
\begin_inset Formula $n+1st$
\end_inset
iterate directly using information about previous iterates and
\begin_inset Formula $\gamma$
\end_inset
, a constant gain coefficient.
\layout Standard
In our quest to devise a corrective scheme we make the following requirements:
\layout Enumerate
\begin_inset Formula $g(x_{n},x_{n1},...,x_{nk})$
\end_inset
must stabilize the previously unstable periodtwo points and destabilize
the previously stable fixed point.
In particular,
\begin_inset Formula $g(x_{n},x_{n1},...,x_{nk})$
\end_inset
should vanish at the periodtwo points to ensure that correction stops
after the system becomes steadystate at these points.
\layout Enumerate
\begin_inset Formula $g(x_{n},x_{n1},...,x_{nk})$
\end_inset
must use as little information about the function (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
) as possible.
This is to allow for our corrective scheme to be less sensitive to experimental
error and to discrepancies between models and reality.
\layout Enumerate
Correction should be used as often as possible.
This requirement is necessary to combat the error and noise inherent in
experimental measurements and therefore the information used in correction.
\layout Standard
A corrective function used by previous authors in a similar context to stabilize
an unstable fixed point involves the difference between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n1st$
\end_inset
iterates:
\layout Standard
\begin_inset Formula \begin{equation}
x_{n+1}=f(x_{n})+\gamma(x_{n}x_{n1})\label{eq:8}\end{equation}
\end_inset
\layout Standard
\align left
One can quickly see, however, that this scheme is illsuited because the
corrective function cannot vanish at the periodtwo points
\begin_inset Foot
collapsed false
\layout Standard
\align left
To fix this problem for the map (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
), one could modify the relation between the
\begin_inset Formula $nth$
\end_inset
and
\begin_inset Formula $n1st$
\end_inset
iterates:
\layout Standard
\align center
\begin_inset Formula $x_{n+1}=f(x_{n})+\gamma(x_{n}+x_{n1})$
\end_inset
\layout Standard
\align left
Such a function, however, requires that the periodtwo points be symmetric
about zero.
If the periodtwo points weren't symmetric about zero but about some nonzero
fixed point, introducing a constant term into the corrective function could
force the function to vanish at the periodtwo points but the fixed point
would have to be known.
This would require too much information about (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
) to be known, thus reducing the generality of the corrective scheme.
\end_inset
.
\layout Standard
In devising a new corrective function, it is helpful to first think about
what information is known about (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
).
Here, we will assume that all that is known about (
\begin_inset LatexCommand \ref{eq:6}
\end_inset
) is that it has one stable fixed point and two unstable periodtwo points.
By definition, the periodtwo points satisfy
\begin_inset Formula $f(f(x))=x$
\end_inset
.
If the system were to reach equilibrium at these points, the difference
between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n2nd$
\end_inset
iterates would be zero.
Therefore, an intuitive suggestion for the corrective function
\begin_inset Formula $g(x_{n},x_{n1},...,x_{nk})$
\end_inset
would be the difference between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n2nd$
\end_inset
iterates:
\layout Standard
\align center
\begin_inset Formula \begin{equation}
x_{n+1}=f(x_{n})+\gamma(x_{n}x_{n2})\label{eq:9}\end{equation}
\end_inset
\layout Standard
Although this seems like a viable option for a corrective scheme, it ultimately
fails.
A proof that this scheme fails will be provided in section 2.2.
In this paper, we prove that the most natural candidates for a corrective
scheme scheme do not work and devise a new scheme that is capable of stabilizin
g unstable periodtwo points.
\layout Section
Corrective Schemes
\layout Subsection
Stability Analysis
\layout Standard
Before analyzing particular corrective schemes, we need to understand how
to diagnose their stability as well as the stability of the map without
feedback.
The fixed points of an iterated map will be stable iff the eigenvalues
\begin_inset Formula $\lambda_{i}$
\end_inset
of the Jacobian satisfy the following requirement
\begin_inset LatexCommand \cite{key7}
\end_inset
:
\layout Standard
\align center
\begin_inset Formula \begin{equation}
\lambda_{i}<1,\forall i\label{eq:10}\end{equation}
\end_inset
\layout Standard
\align left
Note that the periodtwo points of
\begin_inset Formula $f(x_{n})$
\end_inset
are fixed points of the map composed with itself
\layout Standard
\begin_inset Formula \begin{equation}
h(x)=f(f(x))\label{eq:11}\end{equation}
\end_inset
\layout Standard
\align left
For
\begin_inset Formula $1<\mu<1$
\end_inset
the location of the periodtwo points can be solved analytically:
\layout Standard
\begin_inset Formula \begin{equation}
x^{(1)}=\sqrt{1\mu},\; x^{(2)}=\sqrt{1\mu}\label{eq:12}\end{equation}
\end_inset
\layout Standard
\align left
Linearizing (
\begin_inset LatexCommand \ref{eq:11}
\end_inset
) and evaluating at the alternans points yields
\layout Standard
\begin_inset Formula \begin{equation}
h^{\prime}(\pm\sqrt{1\mu})=(912\mu+4\mu^{2})^{2}=(f^{\prime}(\sqrt{1\mu})^{2}=f^{\prime}(\sqrt{1\mu})^{2}\label{eq:13}\end{equation}
\end_inset
\layout Standard
\align left
Note that
\layout Standard
\begin_inset Formula \begin{equation}
f^{\prime}(\pm\sqrt{1\mu})1=4(\mu1)(\mu2)\label{eq:14}\end{equation}
\end_inset
\layout Standard
\align left
It can be seen from (
\begin_inset LatexCommand \ref{eq:14}
\end_inset
) that
\begin_inset Formula $h'(\pm\sqrt{1\mu})>1$
\end_inset
for
\begin_inset Formula $1<\mu<1$
\end_inset
because
\begin_inset Formula $f'(\pm\sqrt{1\mu})>1$
\end_inset
for
\begin_inset Formula $1<\mu<1$
\end_inset
.
\layout Standard
In this paper, we will propose various correction schemes, which we will
later classify as onestep, twostep and threestep correction schemes.
The onestep scheme will simply be (
\begin_inset LatexCommand \ref{eq:9}
\end_inset
), where the difference between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n2nd$
\end_inset
iterates is used as a corrective function multiplied by a constant gain
parameter
\begin_inset Formula $\gamma$
\end_inset
.
\layout Standard
The twostep scheme will use a corrective function of the same form on each
iterate, the difference between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n2nd$
\end_inset
iterates, but will switch between two constant gain parameters
\begin_inset Formula $\gamma_{1}$
\end_inset
and
\begin_inset Formula $\gamma_{2}$
\end_inset
.
The twostep scheme can be stated as follows:
\layout Standard
\begin_inset Formula \[
x_{n+1}=f(x_{n})+\gamma_{1}(x_{n}x_{n2})\]
\end_inset
\begin_inset Formula \[
x_{n+2}=f(x_{n+1})+\gamma_{2}(x_{n+1}x_{n1})\]
\end_inset
where
\begin_inset Formula $n$
\end_inset
is even.
Note that the onestep scheme is a special case of the twostep scheme
with
\begin_inset Formula $\gamma_{1}=\gamma_{2}$
\end_inset
.
The twostep scheme can be rewritten in vector notation:
\layout Standard
\begin_inset Formula \begin{equation}
\left\begin{array}{c}
x_{n+2}\\
x_{n+1}\\
x_{n}\end{array}\right=m(x_{n},\: x_{n1},\: x_{n2})=\left\begin{array}{c}
\mu Y_{n}Y_{n}^{3}+\gamma_{2}(Y_{n}x_{n1})\\
Y_{n}\\
x_{n}\end{array}\right\label{eq:15}\end{equation}
\end_inset
\layout Standard
\noindent
where
\begin_inset Formula $Y_{n}=\mu x_{n}x_{n}^{3}+\gamma_{1}(x_{n}x_{n2})$
\end_inset
\layout Standard
Next, we setup a method for diagnosing the stability of the periodtwo points
under the twostep scheme.
Note that the twostep scheme is a function of
\begin_inset Formula $x_{n}$
\end_inset
,
\begin_inset Formula $x_{n1}$
\end_inset
, and
\begin_inset Formula $x_{n2}$
\end_inset
that generates the next two iterates in the sequence.
Rather than analyze this as a multistep scheme, we may look at the twostep
scheme as a function that takes an ordered triple and advances the indices
by two (
\begin_inset LatexCommand \ref{eq:15}
\end_inset
).
The linearization of the twostep scheme evaluated at the periodtwo points
(
\begin_inset LatexCommand \ref{eq:12}
\end_inset
) is then shown below:
\layout Standard
\begin_inset Formula \begin{equation}
\left\begin{array}{c}
x_{n+2}\\
x_{n+1}\\
x_{n}\end{array}\right=\left\begin{array}{ccc}
\xi & \gamma_{2} & 3\gamma_{1}\gamma_{1}\gamma_{2}2\gamma_{1}\mu\\
3+\gamma_{1}+2\mu & 0 & \gamma_{1}\\
1 & 0 & 0\end{array}\right\left\begin{array}{c}
x_{n}\\
x_{n1}\\
x_{n2}\end{array}\right\label{eq:16}\end{equation}
\end_inset
\layout Standard
\align left
where the 1,1entry is given by
\layout Standard
\align center
\begin_inset Formula $\xi=9233\gamma_{1}3\gamma_{2}+\gamma_{1}\gamma_{2}12\mu+2\gamma_{1}\mu+2\gamma_{2}\mu+4\mu^{2}$
\end_inset
\layout Standard
\align left
The sequence of iterates produced by the twostep scheme will converge to
the periodtwo points if the eigenvalues of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) satisfy (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
).
Note that periodtwo points of
\begin_inset Formula $f(x_{n})$
\end_inset
are fixed points of
\begin_inset Formula $m(x_{n},\: x_{n1},\: x_{n2})$
\end_inset
.
\layout Standard
Next, we will determine a set of requirements on the entries of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) to allow us to determine if this matrix satisfies (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
).
Note that the characteristic polynomial of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
), a 3x3 matrix, will be cubic:
\layout Standard
\begin_inset Formula \begin{equation}
\lambda^{3}a\lambda^{2}b\lambdac\label{eq:17}\end{equation}
\end_inset
\layout Standard
\align left
We can place requirements on the real coefficients of (
\begin_inset LatexCommand \ref{eq:17}
\end_inset
) so that the eigenvalues of the Jacobian (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) satisfy (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
).
The requirements are outlined in the table below:
\layout Standard
\align center
\begin_inset Tabular
\begin_inset Text
\layout Standard
\bar under
Requirement
\end_inset

\begin_inset Text
\layout Standard
\bar under
Description
\end_inset

\begin_inset Text
\layout Standard
1.
a + b + c > 1
\end_inset

\begin_inset Text
\layout Standard
Real eigenvalues will not become larger than 1
\end_inset

\begin_inset Text
\layout Standard
2.
a  b + c < 1
\end_inset

\begin_inset Text
\layout Standard
Real eigenvalues will not become less than 1
\end_inset

\begin_inset Text
\layout Standard
3.
c(a  c)  b > 1
\end_inset

\begin_inset Text
\layout Standard
Complex eigenvalues will not have magnitude greater than 1
\end_inset

\end_inset
\layout Standard
A proof of these requirements is included as an appendix to this paper.
They are sufficient for the stability of any 3x3 linear mapping.
Thus if (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) meets requirements 1  3 outlined in the above table then the twostep
scheme is capable of stabilizing the periodtwo points.
\layout Standard
We will now address the threestep scheme and show that the same requirements
can be used to analyze its Jacobian.
The threestep scheme will also use a corrective function of the same form
on each iterate, the difference between the
\begin_inset Formula $nth$
\end_inset
and the
\begin_inset Formula $n2nd$
\end_inset
iterates, but will alternate between three constant gain parameters
\begin_inset Formula $\gamma_{1}$
\end_inset
,
\begin_inset Formula $\gamma_{2}$
\end_inset
and
\begin_inset Formula $\gamma_{3}$
\end_inset
.
The threestep scheme can be stated as follows:
\layout Standard
\begin_inset Formula \[
x_{n+1}=f(x_{n})+\gamma_{1}(x_{n}x_{n2})\]
\end_inset
\begin_inset Formula \begin{equation}
x_{n+2}=f(x_{n+1})+\gamma_{2}(x_{n+1}x_{n1})\label{eq:18}\end{equation}
\end_inset
\begin_inset Formula \[
x_{n+3}=f(x_{n+2})+\gamma_{3}(x_{n+2}x_{n})\]
\end_inset
where
\begin_inset Formula $n$
\end_inset
is a multiple of three.
The threestep scheme can also be rewritten in vector form:
\layout Standard
\begin_inset Formula \begin{equation}
(x_{n+3},\: x_{n+2},\: x_{n+1})=L(x_{n},\: x_{n1},\: x_{n2})\label{eq:19}\end{equation}
\end_inset
\layout Standard
Note that because the threestep scheme is a function of
\begin_inset Formula $x_{n}$
\end_inset
,
\begin_inset Formula $x_{n1}$
\end_inset
, and
\begin_inset Formula $x_{n2}$
\end_inset
that generates the next three iterates in the sequence, we may look at
the corrective scheme as a function that takes an ordered triple and advances
the indices by three (
\begin_inset LatexCommand \ref{eq:19}
\end_inset
).
The linearization of the threestep scheme is then shown below (actual
entries of the Jacobian are omitted because they are too complicated):
\layout Standard
\begin_inset Formula \begin{equation}
\left\begin{array}{c}
x_{n+3}\\
x_{n+2}\\
x_{n+1}\end{array}\right=\left\begin{array}{ccc}
\frac{\partial x_{n+3}}{\partial x_{n}} & \frac{\partial x_{n+3}}{\partial x_{n1}} & \frac{\partial x_{n+3}}{\partial x_{n2}}\\
\frac{\partial x_{n+2}}{\partial x_{n}} & \frac{\partial x_{n+2}}{\partial x_{n1}} & \frac{\partial x_{n+2}}{\partial x_{n2}}\\
\frac{\partial x_{n+1}}{\partial x_{n}} & \frac{\partial x_{n+1}}{\partial x_{n1}} & \frac{\partial x_{n+1}}{\partial x_{n2}}\end{array}\right\left\begin{array}{c}
x_{n}\\
x_{n1}\\
x_{n2}\end{array}\right\label{eq:20}\end{equation}
\end_inset
\layout Standard
\align left
Because
\begin_inset Formula $f(x^{*})=x^{*}$
\end_inset
for periodtwo points
\begin_inset Formula $x^{*}$
\end_inset
of this particular mapping, the periodtwo points of
\begin_inset Formula $f(x_{n})$
\end_inset
will be periodtwo points of the threestep scheme and fixed points of
the threestep scheme composed with itself.
Thus, the eigenvalues of the corrective scheme composed with itself must
satisfy (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
) for the periodtwo points to be stable.
This is only true, however, when the eigenvalues of (
\begin_inset LatexCommand \ref{eq:20}
\end_inset
) satisfy (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
).
Thus, if it is shown that if the eigenvalues of (
\begin_inset LatexCommand \ref{eq:20}
\end_inset
) satisfy (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
), the periodtwo points can be stabilized using the threestep scheme.
In addition, because this matrix is 3x3, we can use the stability requirements
outlined in the above table to diagnose the stability of the periodtwo
points under this scheme.
\layout Standard
Also note that the onestep scheme is a special case of the threestep scheme
but the twostep scheme is not a special case of the threestep scheme.
This is because there is no way to alternate between two different gain
parameters in the threestep scheme.
\layout Subsection
The Onestep and Twostep Schemes Fail
\layout Standard
Because the onestep scheme is a special case of the twostep scheme, we
show that both schemes fail by proving that the twostep scheme fails.
In this effort, we will show that the eigenvalues of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) are unstable as defined in (
\begin_inset LatexCommand \ref{eq:10}
\end_inset
) for
\begin_inset Formula $1<\mu<1$
\end_inset
.
The Jacobian of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) has the following characteristic polynomial:
\layout Standard
\align center
\begin_inset Formula $\lambda^{3}(9+3\gamma_{1}+3\gamma_{2}\gamma_{1}\gamma_{2}+12\mu2\gamma_{1}\mu2\gamma_{2}\mu4\mu^{2})\lambda^{2}$
\end_inset
\layout Standard
\begin_inset Formula \begin{equation}
(3\gamma_{1}3\gamma_{2}+2\gamma_{1}\gamma_{2}+2\gamma_{1}\mu+2\gamma_{2}\mu)\lambda+\gamma_{1}\gamma_{2}\label{eq:21}\end{equation}
\end_inset
\layout Standard
\align left
To analyze the twostep scheme without correction, we set
\begin_inset Formula $\gamma_{1}=\gamma_{2}=0$
\end_inset
.
The eigenvalues of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) with these requirements can be read off the diagonal:
\begin_inset Formula $\lambda_{1}=912\mu+4\mu^{2}$
\end_inset
,
\begin_inset Formula $\lambda_{2}=0$
\end_inset
, and
\begin_inset Formula $\lambda_{3}=0$
\end_inset
.
As previously mentioned,
\begin_inset Formula $\lambda_{1}>1$
\end_inset
for
\begin_inset Formula $1<\mu<1$
\end_inset
, the range of the subcritical bifurcation.
Thus, before implementing correction, one eigenvalue of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) is larger than 1.
\layout Standard
We will now show that even with nonzero gain parameters (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) will have at least one eigenvalue larger than 1 for
\begin_inset Formula $1<\mu<1$
\end_inset
.
We prove this assertion by showing that the coefficients of (
\begin_inset LatexCommand \ref{eq:21}
\end_inset
) violate stability requirement 1.
Requirement 1 from our stability requirement table is the following:
\layout Standard
\begin_inset Formula \begin{equation}
a+b+c>1\label{eq:22}\end{equation}
\end_inset
\layout Standard
\noindent
where
\begin_inset Formula $a$
\end_inset
,
\begin_inset Formula $b$
\end_inset
, and
\begin_inset Formula $c$
\end_inset
are the coefficients of the characteristic polynomial of a 3x3 matrix as
designated in (
\begin_inset LatexCommand \ref{eq:17}
\end_inset
).
Violation of this requirement means that at least one of the eigenvalues
of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) is greater than positive one.
To determine whether or not requirement (
\begin_inset LatexCommand \ref{eq:22}
\end_inset
) is met, we plug the coefficients of (
\begin_inset LatexCommand \ref{eq:21}
\end_inset
) into their corresponding places in (
\begin_inset LatexCommand \ref{eq:22}
\end_inset
).
The resulting inequality is
\begin_inset Formula \begin{equation}
9+12\mu4\mu^{2}>1\label{eq:23}\end{equation}
\end_inset
\layout Standard
\noindent
Thus, if the above inequality is true for a particularly value of
\begin_inset Formula $\mu$
\end_inset
, none of the eigenvalues of (
\begin_inset LatexCommand \ref{eq:16}
\end_inset
) would be larger than one for that value of
\begin_inset Formula $\mu$
\end_inset
.
We see that this inequality is not satisfied for
\begin_inset Formula $1<\mu<1$
\end_inset
, the range of our model of the subcritical bifurcation, meaning the subcritical
bifurcation cannot be stabilized by the twostep scheme.
Note that (
\begin_inset LatexCommand \ref{eq:23}
\end_inset
), a stability requirement, does not depend on
\begin_inset Formula $\gamma_{1}$
\end_inset
or
\begin_inset Formula $\gamma_{2}$
\end_inset
and therefore cannot be manipulated by correction.
In other words, perturbations that depend on
\begin_inset Formula $\gamma_{1}$
\end_inset
and
\begin_inset Formula $\gamma_{2}$
\end_inset
have no effect on keeping all eigenvalues within the positive one boundary.
Therefore, the twostep scheme fails because it is unable to keep all eigenvalu
es of the mapping less than 1.
Because the onestep scheme is a special case of the twostep scheme, the
failure of the twostep scheme necessarily means that the onestep scheme
also fails.
\layout Subsection
The Threestep Scheme is Successful Under Certain Parameter Settings
\layout Standard
We now analyze the stability of the threestep scheme (
\begin_inset LatexCommand \ref{eq:19}
\end_inset
) for
\begin_inset Formula $1<\mu<1$
\end_inset
.
We start by finding the eigenvalues of the Jacobian (
\begin_inset LatexCommand \ref{eq:20}
\end_inset
) of the threestep scheme without correction.
Setting
\begin_inset Formula $\gamma_{1}=\gamma_{2}=\gamma_{3}=0$
\end_inset
, we see that the eigenvalues before correction are:
\begin_inset Formula $\lambda_{1}=8\mu^{3}36\mu^{2}+54\mu27$
\end_inset
,
\begin_inset Formula $\lambda_{2}=0$
\end_inset
, and
\begin_inset Formula $\lambda_{3}=0$
\end_inset
.
We note that
\begin_inset Formula $\lambda_{1}<1$
\end_inset
for
\begin_inset Formula $1<\mu<1$
\end_inset
, the range of our subcritical bifurcation.
\layout Standard
Next we address whether or not nonzero gain parameters can be used in the
threestep scheme to stabilize the periodtwo points.
Due to the fact that the
\begin_inset Formula $x_{n+3}$
\end_inset
term in (
\begin_inset LatexCommand \ref{eq:19}
\end_inset
) is an order 27 polynomial in
\begin_inset Formula $x_{n}$
\end_inset
, the entries of the Jacobian and the coefficients of the characteristic
polynomial are quite complicated and have been omitted.
Nevertheless, it is possible to analyze the general threestep scheme using
computations similar to those performed in section 2.2.
The inequalities that define the set of all possible
\begin_inset Formula $(\gamma_{1},\gamma_{2},\gamma_{3})$
\end_inset
combinations are available online in a Mathematica notebook file at:
\layout Standard
\align center
www.duke.edu/~maf15/Research/ThreeStepScheme/
\layout Standard
Here we analyze only a minimal case of the threestep scheme that applies
feedback every third iterate, i.e.
\begin_inset Formula $\gamma_{1}=\gamma_{2}=0$
\end_inset
in (
\begin_inset LatexCommand \ref{eq:19}
\end_inset
).
This scheme may be rewritten as:
\layout Standard
\begin_inset Formula \[
y=f(x_{n})\]
\end_inset
\layout Standard
\begin_inset Formula \[
z=f(y)\]
\end_inset
\layout Standard
\begin_inset Formula \begin{equation}
x_{n+1}=f(z)+\gamma_{3}(zx_{n})\label{eq:24}\end{equation}
\end_inset
\layout Standard
\align left
where (
\begin_inset LatexCommand \ref{eq:24}
\end_inset
) is a subsequence of the iterates of a sequence constructed using the threeste
p scheme (
\begin_inset LatexCommand \ref{eq:19}
\end_inset
) with
\begin_inset Formula $\gamma_{1}=\gamma_{2}=0$
\end_inset
.
This subsequence corresponds to elements whose subscript is divisible by
3.
Note that in (
\begin_inset LatexCommand \ref{eq:24}
\end_inset
),
\begin_inset Formula $x_{n+1}$
\end_inset
depends only on
\begin_inset Formula $x_{n}$
\end_inset
.
Because
\begin_inset Formula $x_{n+1}$
\end_inset
is determined by a single variable map, the Jacobian of such a mapping
is a 1x1 matrix whose lone entry is its eigenvalue.
Thus,
\begin_inset Formula $\gamma_{3}$
\end_inset
should be chosen such that the following holds:
\layout Standard
\begin_inset Formula \begin{equation}
1<\frac{\partial x_{n+1}}{\partial x_{n}}_{x_{n}=\pm\sqrt{1\mu}}<1\label{eq:25}\end{equation}
\end_inset
\layout Standard
\noindent
These two inequalities determine the range of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilize the periodtwo points in terms of
\begin_inset Formula $\mu$
\end_inset
.
This range of
\begin_inset Formula $\gamma_{3}$
\end_inset
is shown in dark gray in Figure 5.
\begin_inset Formula $\gamma_{3}$
\end_inset
values above the dark gray area result in eigenvalues larger than 1 whereas
\begin_inset Formula $\gamma_{3}$
\end_inset
values below the dark gray area result in eigenvalues less than 1.
Note that as
\begin_inset Formula $\mu$
\end_inset
decreases, the range of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilizes the periodtwo points decreases.
Thus, the farther away from the bifurcation point, the smaller the range
of effective gain parameters.
This analysis suggests that stabilizing a subcritical bifurcation in experiment
with the threestep scheme will be easier closer to the bifurcation point.
\layout Standard
\begin_inset Float figure
wide false
collapsed false
\layout Standard
\align center
\begin_inset Graphics
filename C:/mstuff/images/mincase3.bmp
lyxscale 60
scale 45
keepAspectRatio
\end_inset
\layout Caption
Range of
\begin_inset Formula $\gamma_{3}$
\end_inset
for
\begin_inset Formula $1<\mu<1$
\end_inset
.
The dark gray region corresponds to values of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilize the periodtwo points.
Above this region, the eigenvalue at the periodtwo points is greater than
1 and below this region the eigenvalue is less than 1.
The light gray region corresponds to values of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilize the fixed point.
Above this region, the eigenvalue at the fixed point is greater than 1
and below this region the eigenvalue is less than 1.
\end_inset
\layout Standard
Another point of interest is the stability of the fixed point
\begin_inset Formula $x=0$
\end_inset
.
This solution will be stable as long as the following holds:
\layout Standard
\begin_inset Formula \begin{equation}
1<\frac{\partial x_{n+1}}{\partial x_{n}}_{x_{n}=0}<1\label{eq:26}\end{equation}
\end_inset
\layout Standard
\noindent
The fixed point will be unstable otherwise.
These two inequalities determine the range of
\begin_inset Formula $\gamma_{3}$
\end_inset
that preserves the stability of the fixed point for particular values of
\begin_inset Formula $\mu$
\end_inset
.
This range of possible
\begin_inset Formula $\gamma_{3}$
\end_inset
is shown in light gray in Figure 5.
\begin_inset Formula $\gamma_{3}$
\end_inset
values above the light gray area result in eigenvalues larger than 1 whereas
\begin_inset Formula $\gamma_{3}$
\end_inset
values below the light gray area result in eigenvalues less than 1.
In determining a scheme to stabilize a periodtwo point, we want the fixed
point to be unstable to ensure that the system is driven to equilibrium
at the periodtwo points rather than at the fixed point.
Therefore, we will want to choose
\begin_inset Formula $\gamma_{3}$
\end_inset
outside of the light gray region.
\layout Standard
Fortunately, for most
\begin_inset Formula $\mu$
\end_inset
the range of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilizes the periodtwo points and the range of
\begin_inset Formula $\gamma_{3}$
\end_inset
that preserves the stability of the fixed point do not overlap.
There exists a range of
\begin_inset Formula $\mu$
\end_inset
for which values of
\begin_inset Formula $\gamma_{3}$
\end_inset
that stabilize the periodtwo points also preserve the stability of the
fixed point.
This range of
\begin_inset Formula $\mu$
\end_inset
is approximately
\begin_inset Formula $1<\mu<.826155$
\end_inset
and the corresponding region is shaded black.
In this range of
\begin_inset Formula $\mu$
\end_inset
, a set of periodsix points appears, with two of the periodsix points
lying between the periodtwo points and the fixed point.
These periodsix points are not critical points of the original map and
are created by implementing the threestep scheme.
These new critical points are also unstable throughout the range
\begin_inset Formula $1<\mu<.826155$
\end_inset
.
\layout Section
Conclusions
\layout Standard
In this paper we devised a method of stabilizing a subcritical bifurcation
using feedback.
We determined that a minimum of a threestep scheme is required when using
the difference between the
\begin_inset Formula $nth$
\end_inset
and
\begin_inset Formula $n2nd$
\end_inset
iterates as the corrective function.
Although a minimal case of the threestep scheme is analyzed in this paper,
feedback may be applied on each iterate, which is in accordance with the
guideline set in section 1.3.
In addition to stabilizing the periodtwo points, the threestep scheme
destabilizes the fixed point, uses little information about the system
and has a corrective function that vanishes at the periodtwo points.
Thus, it meets all the guidelines we established before formulating a feedback
scheme.
\layout Standard
Our analysis suggests that this scheme may be most effectively implemented
closer to the bifurcation point.
The range of successful gain parameters far from the bifurcation point
is very small.
Also, closer to the bifurcation point convergence to the periodtwo points
is less sensitive to the initial value
\begin_inset Formula $x_{0}$
\end_inset
of the sequence.
Further studies should seek to increase the robustness of the corrective
scheme by decreasing its sensitivity to initial
\begin_inset Formula $x_{0}$
\end_inset
of the sequence.
\layout Section*
Acknowledgements
\layout Standard
We would like to thank Dr.
David Schaeffer, who served as sponsor for this paper, for his assistance
throughout the process of researching this topic and writing the paper.
We would also like to thank Dr.
Kraines and the PRUV program for setting up the research program and funding
through which this research took place.
\layout Bibliography
\bibitem [1]{key7}
Guckenheimer, J.
and P.
Holmes.
Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields.
(1983) SpringerVerlag, New York
\layout Bibliography
\bibitem [2]{key2}
Hall, G.
M., and D.
J.
Gauthier (2002).
Experimental Control of Cardiac Muscle Alternans.
\shape italic
Phys.
Rev.
Lett
\shape default
.
88, 198102.
\layout Bibliography
\bibitem [3]{key5}
Mitchell, C.
C., and D.
G.
Schaeffer (2003).
A TwoCurrent Model for the Dynamics of Cardiac Membrane.
\shape italic
Bulletin Math Bio
\shape default
.
65, 767793.
\layout Bibliography
\bibitem [4]{key4}
Nolasco, J.
and R.
Dahlen (1968).
A Graphic Method for the Study of Alternation in Cardiac Action Potentials.
\shape italic
J.
Appl.
Physiol
\shape default
.
25, 191196.
\layout Bibliography
\bibitem [5]{key6}
Strogatz, S.
H.
\shape italic
Nonlinear Dynamics and Chaos
\shape default
.
(1994) Westview Press, Massachusetts
\layout Standard
\begin_inset ERT
status Collapsed
\layout Standard
\backslash
newpage
\end_inset
\layout Section*
Appendix: Proof of Stability Requirements
\layout Standard
\series bold
Theorem
\series default
: The eigenvalues of a real 3x3 matrix
\begin_inset Formula $A$
\end_inset
are contained in the open unit disk in the complex plane if the coefficients
of its characteristic polynomial
\begin_inset Formula $p(\lambda)$
\end_inset
below satisfy
\layout List
\labelwidthstring 00.00.0000
1.
\begin_inset Formula $a+b+c>1$
\end_inset
\layout List
\labelwidthstring 00.00.0000
2.
\begin_inset Formula $ab+c<1$
\end_inset
\layout List
\labelwidthstring 00.00.0000
3.
\begin_inset Formula $c(ac)b>1$
\end_inset
\layout Standard
where
\begin_inset Formula $det(A\lambda I)=p(\lambda)=\lambda^{3}a\lambda^{2}b\lambdac$
\end_inset
.
\begin_inset ERT
status Collapsed
\layout Standard
\backslash
bigskip
\end_inset
\layout Standard
\noindent
\shape italic
Proof
\shape default
.
If
\begin_inset Formula $a=b=c=0$
\end_inset
, each
\begin_inset Formula $\lambda=0$
\end_inset
.
Note that inequalities 13 are satisfied when all three coefficients vanish.
As coefficients of the characteristic polynomial are varied, the roots
move continuously.
The roots will not exit the unit disk unless
\layout List
\labelwidthstring 00.00.0000
i.
A real root passes through
\begin_inset Formula $\lambda=1$
\end_inset
.
\layout List
\labelwidthstring 00.00.0000
ii.
A real root passes through
\begin_inset Formula $\lambda=1$
\end_inset
.
\layout List
\labelwidthstring 00.00.0000
iii.
A pair of complex conjugate roots crosses the unit circle off the real
axis.
\layout Standard
\noindent
Conditions 1, 2 and 3 derive from excluding possibilities i, ii, and iii,
respectively.
\layout Standard
If
\begin_inset Formula $\lambda=1$
\end_inset
is a root, then
\begin_inset Formula $a+b+c=1$
\end_inset
.
Thus if
\begin_inset Formula $a$
\end_inset
,
\begin_inset Formula $b$
\end_inset
and
\begin_inset Formula $c$
\end_inset
are varied away from zero but
\begin_inset Formula $a+b+c$
\end_inset
remains greater than
\begin_inset Formula $1$
\end_inset
, no real root can pass through
\begin_inset Formula $\lambda=1$
\end_inset
.
\layout Standard
Similarly if
\begin_inset Formula $\lambda=1$
\end_inset
is a root, then
\begin_inset Formula $ab+c=1$
\end_inset
.
Thus inequality 2 prevents a real root from passing through
\begin_inset Formula $\lambda=1$
\end_inset
.
\layout Standard
If
\begin_inset Formula $p(\lambda)$
\end_inset
has complex conjugate roots on the unit circle, then
\begin_inset Formula $p(\lambda)$
\end_inset
may be factored as:
\layout Standard
\align center
\begin_inset Formula $(\lambdae^{i\theta})(\lambdae^{i\theta})(\lambdar)=\lambda^{3}(2cos\theta+r)\lambda^{2}+(2rcos\theta+1)\lambdar$
\end_inset
\layout Standard
\noindent \align left
where r
\begin_inset Formula $\in\mathbb{R}$
\end_inset
.
Matching coefficients in
\layout Standard
\align center
\begin_inset Formula $\lambda^{3}(2cos\theta+r)\lambda^{2}+(2rcos\theta+1)\lambdar$
\end_inset
\layout Standard
\align left
and
\layout Standard
\align center
\begin_inset Formula $\lambda^{3}a\lambda^{2}b\lambdac$
\end_inset
,
\layout Standard
\noindent
we obtain the equations
\layout Standard
\align center
\begin_inset Formula $a=2cos\thetar$
\end_inset
\layout Standard
\align center
\begin_inset Formula $b=2rcos\theta+1$
\end_inset
\layout Standard
\align center
\begin_inset Formula $c=r$
\end_inset
\layout Standard
\noindent \align left
Eliminating
\begin_inset Formula $r$
\end_inset
and
\begin_inset Formula $\theta$
\end_inset
from these equations we deduce that
\begin_inset Formula $c(ac)b=1$
\end_inset
.
Thus inequality 3 prevents complex conjugate roots from crossing the unit
circle.
\the_end