\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2009(2009), No. 53, pp. 1--21.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}
\begin{document}
\title[\hfilneg EJDE-2009/53\hfil Chaos in an extended van der Pol's equation]
{Bifurcation routes to chaos in an extended van der Pol's equation
applied to economic models}
\author[L. P\v ribylov\'a\hfil EJDE-2009/53\hfilneg]
{Lenka P\v ribylov\'a}
\address{Lenka P\v ribylov\'a \newline
Dept. of Applied Mathematics\\ Masaryk University\\
Jan\'a\v ckovo n\'am. 2a\\
602 00 Brno, Czech Republic}
\email{pribylova@math.muni.cz}
\thanks{Submitted April 21, 2008. Published April 17, 2009.}
\subjclass[2000]{70K50, 37D45}
\keywords{Hopf bifurcation; period doubling; chaos}
\begin{abstract}
In this paper a 3-dimensional system of autonomous differential equations
is studied. It can be interpreted as an idealized macroeconomic model
with foreign capital investment introduced in \cite{vosvrda} or
an idealized model of the firm profit introduced in \cite{bouali}.
The system has three endogenous variables with only one non-linear term
and can be also interpreted as an extended van der Pol's equation.
It's shown that this simple system covers several types of bifurcations:
both supercritical and subcritical Hopf bifurcation and generalized
Hopf bifurcation as well, the limit cycle exhibits period-doubling
bifurcation as a route to chaos. Some results are analytical and
those connected with chaotic motion are computed numerically
with continuation programs Content, Xppaut and Maple. We present conditions
for stability of the cycles, hysteresis, explore period doubling and using Poincar\'e mapping
show a three period cycle that implies chaos.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\section{Introduction}
In this paper we consider the autonomous system of three
differential equations
\begin{equation}
\label{eq:system1}
\begin{gathered}
\dot x=ay+px(\kappa-y^2),\\
\dot y=v(x+z),\\
\dot z=mx-ry,
\end{gathered}
\end{equation}
where $x$, $y$, $z$ are real endogenous variables and $a$, $p$,
$v$, $m$ and $r$ real exogenous parameters. This system may be
interpreted as an extension of the van der Pol's equation
\begin{equation}
\label{eq:vanderpol}
\begin{gathered}
\dot x=ay+px(\kappa-y^2),\\
\dot y=vx.
\end{gathered}
\end{equation}
\subsection{Two economic models}
Vosvrda \cite{vosvrda} introduced an idealized macroeconomic
model with foreign capital investment in the form
\begin{equation}
\begin{gathered}
\dot S=aY+pS(\kappa-Y^2),\\
\dot Y=v(S+F),\\
\dot F=mS-rY,
\end{gathered}
\end{equation}
where $S(t)$ are savings of households, $Y(t)$ is the Gross
Domestic Product, $F(t)$ is the foreign capital inflow, $t$ is the
time and dot denotes the derivative with respect to $t$. Positive
parameters represent corresponding ratios: $a$ is the variation of
the marginal propensity to savings, $p$ is the ratio of the
capitalized profit, $v$ is the output/capital ratio, $\kappa$ is
the potential GDP (it can be set to $1$, as a unit of $GDP$ ---
$Y$, $S$, $F$ then represent the percentage part of the potential
$GDP$), $m$ is the capital inflow/savings ratio and $r$ is is the
debt refund/output ratio. From the economic point of view, the
condition
\begin{equation}
\label{cond}
\frac{a}{vr}>1.
\end{equation}
guarantees the ability of the economy to refund the debt. Another
--- stronger condition --- can be $a > r$ ($v \in (0,1)$,
normally much lower then 1, see \cite{dadda}), which implies the
condition \eqref{cond}. The condition $a>r$ together with
\eqref{cond} was used in \cite{vosvrda}. We will see that the
equilibria are stable for $a > r$, although orbitally stable and
even unstable cycles can occur in this economically "normal" case.
In the section \ref{near zero} we let some of the parameters cross
the zero axis and accept their negativity to explain behaviour of
trajectories in the positive neighbourhood of zero.
>From the mathematical point of view, S. Bouali \cite{bouali}
introduced the same system as an idealized economic model of firm
profit in the form
\begin{equation}
\begin{gathered}
\dot R=aP+pR(\kappa-P^2),\\
\dot P=v(R+F),\\
\dot F=mR-rP,
\end{gathered}
\end{equation}
where $P$ is a firm profit, $R$ are reinvestments and $F$
represent debts, coefficients $a$, $p$, $v$, $m$ and $r$ are
corresponding rates or proportions.
\section{Equilibria and its stability}
The system \eqref{eq:system1} is antisymmetric. If
$(x(t),y(t),z(t))$ is a solution of \eqref{eq:system1}, than
$(-x(t),-y(t),-z(t))$ is also its solution. Solving the system
\begin{equation}
\label{eq:equil}
\begin{gathered}
0=ay+px(\kappa-y^2),\\
0=v(x+z),\\
0=mx-ry,
\end{gathered}
\end{equation}
we find, that the system \eqref{eq:system1} has three equilibria:
$E_0=(x_0,y_0,z_0)=(0,0,0)$,
\begin{gather*}
E_1=(x_1,y_1,z_1)=\Bigl (
-\sqrt{\frac{amr+p\kappa r^2}{pm^2}},\sqrt{\frac{amr+p\kappa
r^2}{pm^2}},\sqrt{\frac{am+p\kappa r}{pr}}
\Bigr ),\\
E_2=(x_2,y_2,z_2)=\Bigl ( \sqrt{\frac{amr+p\kappa
r^2}{pm^2}},-\sqrt{\frac{amr+p\kappa
r^2}{pm^2}},-\sqrt{\frac{am+p\kappa r}{pr}} \Bigr ).
\end{gather*}
We will study the equilibria $E_0$ and $E_1$, since results for $E_2$ are
analogous to the antisymmetric equilibrium $E_1$. Since the
Jacobian matrix has the form
\begin{equation}
\label{eq:jacobi}
J=
\begin{pmatrix}
p(\kappa-y^2)&a-2pxy&0\\
v&0&v\\
m&-r&0
\end{pmatrix},
\end{equation}
we can get the corresponding eigenvalues by solving cubic
equations. These formulas are very complicated, but we can use
them for numerical examples. On the other hand, we can get
information about stability of both equilibria right from the
characteristic polynomial
\begin{equation}
\label{eq:charpoly}
p(\lambda)=\lambda^3+a_1\lambda^2+a_2\lambda+a_3,
\end{equation}
where $a_1=p(y^2-\kappa)$, $a_2=v(r-a+2pxy)$ and
$a_3=v(pry^2+2mpxy-p\kappa r-am)$:
\noindent\textbf{1.} Equilibrium $E_0=(0,0,0)$:
Coefficients of the characteristic polynomial are
\begin{equation}
\label{polyE0} a_1=-p\kappa<0, \; a_2=v(r-a), \;
a_3=-v(rp\kappa+am)<0.
\end{equation}
From the Hurwitz criterion it yields that in the case $a \geq r$,
the characteristic polynomial has one positive root. In the case
that $a0, \;
a_2=v(r+a+\frac{2rp\kappa}{m})>0, \; a_3=2v(rp\kappa + am)>0.
\end{equation}
\begin{lemma}
\label{hurwitz} Characteristic polynomial is Hurwitzian (that is
all eigenvalues have negative real parts) if and only if $a>r$.
\end{lemma}
\begin{proof}
The characteristic polynomial is Hurwitzian if and only if the
determinant
\begin{equation*}
D_3=\begin{vmatrix} a_1 & 1 & 0\\
a_3 & a_2 & a_1\\
0 & 0 & a_3
\end{vmatrix}
\end{equation*}
and all its main subdeterminants are positive. The subdeterminant
$D_1=a_1=am/r>0$, subdeterminant $D_2 = \begin{vmatrix} a_1 & 1 \\
a_3 & a_2 \end{vmatrix} =a_1a_2 -a_3$ and the determinant
$D_3=a_3(a_1a_2-a_3)$ are positive if and only if
\begin{align*}
a_1a_2-a_3 &= \frac{am}{r}v(r+a+\frac{2rp\kappa}{m})-2v(rp\kappa +
am) \\
&= (a-r)\Bigl ( \frac{amv}{r} + 2vp \kappa \Bigr ) > 0,
\end{align*}
that is if and only if $a>r$.
\end{proof}
Both the equilibria cannot change its stability on the real axis
of the eigenvalues plane, since the characteristic polynomial
cannot have positive real root. It is possible, if the eigenvalues
crossed the imaginary axis. This type of bifurcation is called
Hopf, and it will be discussed in the next section.
\section{Hopf bifurcation and stability of limit cycles}
If we look for the Hopf bifurcation, we have to find two purely
imaginary eigenvalues $\lambda_{1,2}=i\omega$ of the
characteristic polynomial. Denote the third real eigenvalue
$\lambda_3$. Since substituting into \eqref{eq:charpoly} gives
$$
(\lambda-i\omega)(\lambda+i\omega)(\lambda-\lambda_3)=
\lambda^3-\lambda^2\lambda_3+\lambda\omega^2-\omega^2\lambda_3,
$$
the necessary condition for Hopf bifurcation is
\begin{equation}
\label{eq:hopf}
a_1a_2=a_3.
\end{equation}
Solving a system of five equations \eqref{eq:equil},
\eqref{eq:hopf} with substitution from \eqref{eq:charpoly} and
$a_2=\omega^2$ according to three variables $x, y, z$ and
two parameters $a, m$, we get two solutions. For the trivial
equilibrium $E_0$, we necessarily have
\begin{equation}
\label{mpk}
m=-p\kappa.
\end{equation}
We can exclude this case since all parameters are positive for now
(we will return to this bifurcation later). The next solution is
equilibrium
\begin{equation}
\label{hopfeq}
E_{1,2}=\Bigl (\pm
\frac{r}{m}\sqrt{\frac{m+p\kappa}{p}}, \pm
\sqrt{\frac{m+p\kappa}{p}},\mp
\frac{r}{m}\sqrt{\frac{m+p\kappa}{p}} \Bigr )
\end{equation}
with parameters
satisfying $a=r$ and $m=\frac{2vrp\kappa}{\omega^2-2vr}$, which
gives
\begin{equation}
\label{omega}
\omega=\sqrt{\frac{2vr(p\kappa+m)}{m}}.
\end{equation}
From the previous results we know that the third real eigenvalue
has to be negative and we get its value $\lambda_3=-m$.
We will show that the Hopf bifurcation appears for parameter $a$
while crossing the critical value $r$. Using the implicit function
theorem we can compute the derivative of the complex eigenvalue
$\lambda$ with respect to $a$ for the equilibrium $E_1$ (using
\eqref{polyE1}):
\begin{equation}
\label{dlambda}
\frac{d\lambda}{da}=-\frac{\frac{dp(\lambda)}{da}}{\frac{dp(\lambda)}{d\lambda}}=
-\frac{\lambda^2\frac{m}{r}+\lambda v +2vm}{3\lambda^2+ 2 \lambda
\frac{am}{r} + v(r+a)+\frac{2vrp\kappa}{m} }.
\end{equation}
Substituting $a=r$ and $\lambda_{1,2}=\pm i\omega$ into \eqref{dlambda},
we evaluate
\begin{equation}
\mathop{\rm Re} \frac{d\lambda}{da} |_{a=r}=-\frac{(2mv-\frac{m \omega^2}{r})(-3\omega^2+2vr+\frac{2vrp\kappa}{m}) + 2v\omega^2m}
{(-3\omega^2+2vr+\frac{2vrp\kappa}{m})^2+4\omega^2m^2}
\end{equation}
The denominator of this fraction is positive, substituting \eqref{omega}
into the nominator we have
\begin{equation}
\label{dlambda2}
\mathop{\rm sgn} \mathop{\rm Re} \frac{d\lambda}{da} |_{a=r} = \mathop{\rm sgn}
\frac{-4rv^2(m+p\kappa)(m+2p\kappa)}{m}<0.
\end{equation}
The transversality condition for Hopf bifurcation is fulfilled and
according to \cite{kuznecov} the Hopf bifurcation gives rise to a
limit cycle near the equilibria at the critical value $a=r$.
According to the Lemma \eqref{hurwitz}, the equilibrium $E_{1,2}$
change its unstability to stability as the parameter $a$ grows
and crosses the critical value $r$. Locally we can educe the same
result from \eqref{dlambda2}. Clearly, the following theorem
holds:
\begin{theorem}
\label{ar}
Equilibria $E_{1,2}$ of the system \eqref{eq:system1} are stable
for $a>r$ and lose their stability if a parameter crosses the Hopf
bifurcation manifold $a=r$.
\end{theorem}
Note that the condition \eqref{cond} is satisfied on the Hopf
bifurcation hyperplane, unlike $a>r$ used in \cite{vosvrda}.
In contrast to \cite{vosvrda}, we can prove that the Hopf
bifurcation (in both previous cases) can give rise to both the
stable and the unstable cycles. Following the projection method
for center manifold computation (see \cite{kuznecov}), we can get
formula for the first Lyapunov coefficient $l_1$. Negative sign of
this coefficient implies stability, positive sign unstability of
the arising limit cycle near the equilibrium. Symbolically this
formula is very complicated, so we computed the Hopf bifurcation
curves numerically using the continuation program Content.
The bifurcation border $GH$ of the generalized Hopf bifurcation
($l_1=0$) has typical form presented on figures \ref{fig1},
\ref{fig2} and \ref{fig3}. The parameter $\kappa$ can be set to 1
as it was explained earlier, parameters $v$, $p$, $m$ correspond
to ratios and the bifurcation diagrams are similar for all
parameter values in the whole interval $\langle 0, 1 \rangle$
(here and below, the concrete values of the parameters were chosen
arbitrarily, but to follow examples in \cite{vosvrda}, where
misleading results were made).
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1} %GHav.ps
\end{center}
\caption{$a=r$, $p=0.1$, $\kappa=1$, $m=0.19$.} \label{fig1}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2} % GHap.ps
\end{center}
\caption{$a=r$, $v=0.5$, $\kappa=1$, $m=0.19$.} \label{fig2}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig3} %GHam.ps
\end{center}
\caption{$a=r$, $v=0.5$, $\kappa=1$, $p=0.1$.} \label{fig3}
\end{figure}
For illustrating the diagrams, we present typical diagram \ref{fig4}
of the Hopf bifurcation in the plane $a,r$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig4} %ar.ps
\end{center}
\caption{The typical Hopf bifurcation curve.} \label{fig4}
\end{figure}
\section{Folding of cycles and period doubling}
In this section typical bifurcation diagrams will be presented.
The parameters in examples are chosen to correspond the previous
bifurcation diagrams \ref{fig1}, \ref{fig2}, \ref{fig3} and
\ref{fig4}. By a convention, solid lines represent stable
equilibria, dashed lines are unstable equilibria, solid circles
correspond to orbitally stable cycles and empty circles to the
unstable ones.
\subsection*{Example 1.}
Let us focus on the example with parameters $r=0.25$, $p=0.1$,
$v=0.5$, $\kappa=1$ and $m=0.19$. From \eqref{dlambda2} it yields
that $\mathop{\rm sgn} \frac{d \mathop{\rm Re} \lambda}{da} \doteq -0.01$ and as parameter
$a$ grows and cross $r=0.25$, $\mathop{\rm Re} \lambda$ decreases very
slowly. Since the Lyapunov coefficient is $l_1 \doteq -0.0318$, it
give rise to an orbitally stable cycle from the stable
equilibrium. The branches of the periodic solutions are presented
on the figure \ref{fig5}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig5} % doub05.ps
\end{center}
\caption{The typical Hopf bifurcation curve with arising stable
cycle and period doubling of the cycle for $a \doteq 0.09537$.}
\label{fig5}
\end{figure}
\subsection*{Example 2.}
On the other hand, for parameters $r=0.25$, $p=0.01$,
$v=0.031847$, $\kappa=1$ and $m=0.19$ (this example was already
shown in \cite{vosvrda}, example 3, but it was improperly
interpreted as supercritical Hopf bifurcation with a stable
mounting cycle) the Lyapunov coefficient is $l_1 \doteq 0.000278$
and so the subcritical Hopf bifurcation takes place here, an
unstable cycle mount from the unstable equilibrium. On the figure
\ref{fig6} the branches of the periodic solutions are presented.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig6} % doubling.ps
\end{center}
\caption{The typical Hopf bifurcation curve with arising unstable
cycle.} \label{fig6}
\end{figure}
You can see the Hopf bifurcation critical point ($HB$) for
$a=0.25$ (period $48.53$) and unstable cycle branch around a
stable equilibrium $E_1$ nearby. This branch folds ($LP$) at $a
\doteq 0.2695$ (period 51.92) and changes its stability. This is
the reason, why some solutions tend to a stable cycle and others
to a stable equilibrium. This coexistence of two attractors (a
stable fixed point and a stable cycle) is called hysteresis. When
we continue further, the cycle on the stable cycle branch
bifurcates by period doubling ($PD_1$) at $a \doteq 0.1693$
(period 61.78). The unstable part of the branch give rise to a
chaotic motion (a strange attractor). This phenomenon will be
discussed in the next section. A stable cycle (period 69.83) mount
again for $a \doteq 0.03342$ from the stable branch ($PD_2$) near
zero.
The stable limit cycle surrounding the unstable one near the
subcritical bifurcation point together with the trajectory tending
to the stable equilibrium $E_1$ is shown on the figure \ref{fig7}.
It demonstrates the above Example 2. for $a=0.265$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig7} % stablesub.ps
\end{center}
\caption{The stable cycle and a trajectory (converging to a stable
equilibrium) repelling from the unstable cycle near the
subcritical Hopf bifurcation due to folding. For a similar figure
you can use parameters from the Example 2. and $a \in
(0.25,0.2695)$ and initial conditions near $E_1$ for a trajectory
tending to the stable equilibrium $E_1$ and for the stable cycle
use some farther initial condition, for example $x(0)=10$,
$y(0)=10$, $z(0)=10$ and $t > 6000$.} \label{fig7}
\end{figure}
\section{Period doubling route to chaos}
Period doubling of the limit cycle can be shown in the previous
Example 2 (using continuation programs, for example Xppaut). First
branch bifurcates for $a_1 \doteq 0.1692534$ (see figure
\ref{fig6}), the next doubling takes place for $a_2 \doteq
0.1415250$ (see figure \ref{fig8}) and then for $a_3 \doteq
0.1346818$ (see figure \ref{fig9}), $a_4 \doteq 0.1331614$,
$a_5\doteq 0.1328327$, $a_6 \doteq 0.1327623$,
$a_7 \doteq 0.1327473$, \dots .
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig8} % doubling2.ps
\end{center}
\caption{Second branch of the period doubling.} \label{fig8}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig9} % doubling3.ps
\end{center}
\caption{Third branch of the period doubling.} \label{fig9}
\end{figure}
The ratio of distances between period doubling bifurcation points
is $\delta_1 = \frac{a_1-a_2}{a_2-a_3} \doteq
\frac{0.0277284}{0.0068432} \doteq 4.052$,
$\delta_2 = \frac{a_2-a_3}{a_3-a_4} \doteq \frac{0.0068432}{0.0015204} \doteq
4.501$, $\delta_3 \doteq 4.625$, $\delta_4 \doteq \delta_5 \doteq
4.669$ which is already pretty close to Feigenbaum constant.
This period doubling of the limit cycle lead to a chaotic motion
in a bounded region --- existence of a strange attractor. A
chaotic region is characterized by a sensitive dependence on
initial conditions, and arbitrarily close initial conditions lead
to evolutions that diverge exponentially fast with time. This
divergence is characterized by the maximal Lyapunov exponent (see
precise definition in \cite{alligood} for example). Briefly
speaking, Lyapunov exponents measure average rate of divergence of
nearby trajectories (for $\lambda > 0$) or convergence (for
$\lambda < 0$) respectively. In our case of a three dimensional
system, we have three Lyapunov exponents (for a trajectory or an
attractor respectively), for a fixed point the signs are
$(-,-,-)$, for a stable cycle $(0,0,-)$ and for a strange
attractor $(+,0,-)$. The maximal Lyapunov exponent is negative for
a fixed point, zero for a stable cycle and positive for chaotic
strange attractor.
Changing parameter $a$ from $0$ to $0.3$ (200 stps) in the Example
2, we computed the maximal Lyapunov exponent as a function of $a$.
We used program Xppaut (Xppaut computes the maximal Lyapunov
exponent along a computed trajectory by linearizing in each point
of the trajectory, advancing one time step using a normalized
vector, computing the expansion, and summing the log of the
expansion. The average of this over the trajectory is a rough
approximation of the maximal Lyapunov exponent. For further
details see methods presented in \cite{benettin}, \cite{wolf} or
\cite{ruelle}.) The maximal Lyapunov exponent is computed
numerically (using Xppaut) in finite number of points on
trajectory converging to the equilibrium, periodic orbit or a
strange attractor, starting at some initial point within the basin
of attraction with first amount of transient iterations being
discarded to converge to an attractor. The results are presented
on the next figures.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig10} % lyapexpa.ps
\end{center}
\caption{Estimation of the maximal Lyapunov exponent for $a \in
\langle 0, \, 0.3 \rangle$, $r=0.25, \, p=0.01, \, v=0.031847, \,
\kappa=1$ and $m=0.19$. Initial conditions for all computed
trajectories: $x(0)=0.45$, $y(0)=0.8$ and $z(0)=0.35$ for $t \in
\langle 3000, 6000 \rangle$} \label{fig10}
\end{figure}
On the figure \ref{fig10} you can see a range, where the maximal
Lyapunov exponent is positive and that implies existence of a
chaotic strange attractor for these values of parameters.
Changing the parameter $v$ from $0$ to $0.3$ for $a = 0.1$ we can
see a range, where the maximal Lyapunov exponent is positive again
on the figure \ref{fig11}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig11} % lyapexpv.ps
\end{center}
\caption{Estimation of the maximal Lyapunov exponent for $v \in
\langle 0, \, 0.3 \rangle$, $a=0.1, \, r=0.25, \, p=0.01, \,
\kappa=1$ and $m=0.19$. Initial conditions for all computed
trajectories: $x(0)=0.45$, $y(0)=0.8$ and $z(0)=0.35$ for $t \in
\langle 3000, 6000 \rangle$} \label{fig11}
\end{figure}
Changing parameter $p$ from $0$ to $1$ for fixed parameters
$a=0.1$, $r=0.25$, $v=0.031847$, $\kappa=1$ and $m=0.19$ and
initial conditions $x=0.45$, $y=0.8$, $z=0.35$ and $t \in \langle
3000, 6000 \rangle$, we got different results, since the maximal
Lyapunov exponent is always either positive or zero (numerical
programs of course give negative results, but it's due to the
finite number of computed points and the numerical methods), that
is either chaotic or periodic trajectories occur. For example, for
$p=0.2$ and $p=0.24$ the maximal Lyapunov exponent is positive,
but for $p \doteq 0.22$ zero (see figure \ref{fig12}) and the
trajectory converges to a stable limit cycle (see figures
\ref{fig13} and \ref{fig14}). This phenomenon is caused by
period-doubling route to chaos and a fractal structure of the
bifurcation diagram with stable areas (stable areas of the
Poincar\'e section more precisely). The maximal Lyapunov exponent
is tending to zero from $p \doteq 0.67$ (see figure \ref{fig15})
and periodic orbits occur. Maximal Lyapunov exponent takes
positive values for a big range of the $(0,1)$ interval (see also
figure \ref{fig21}) and that means this parameter cannot be used
for controlling the system.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig12} % lyapexpp2.ps
\end{center}
\caption{Estimation of the maximal Lyapunov exponent for
$p \in \langle 0.1, 0.3 \rangle$, $a=0.1$, $r=0.25$, $v=0.031847$,
$\kappa=1$ and $m=0.19$. Initial conditions for all computed
trajectories: $x(0)=0.45$, $y(0)=0.8$ and $z(0)=0.35$ for
$t \in \langle 3000, 6000 \rangle$} \label{fig12}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig13} % phase02.ps
\end{center}
\caption{Chaotic orbit for positive maximal Lyapunov exponent for
$p=0.2$.} \label{fig13}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig14} % phase022.ps
\end{center}
\caption{Periodic orbit for $p=0.22$} \label{fig14}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig15} % lyapexpp.ps
\end{center}
\caption{Estimation of the maximal Lyapunov exponent for
$p \in \langle 0.67, 0.8 \rangle$, $a=0.1$, $r=0.25$, $v=0.031847$,
$\kappa=1$ and $m=0.19$. Initial conditions for all computed
trajectories: $x(0)=0.45$, $y(0)=0.8$ and $z(0)=0.35$ for
$t \in \langle 3000, 6000 \rangle$} \label{fig15}
\end{figure}
The last figure \ref{fig16} presents dependence of the maximal
Lyapunov exponent on the parameter $m$. Maximal Lyapunov exponent
takes positive values as $m$ goes to 1 (see also figure
\ref{fig22}) and that means this parameter cannot be used for
controlling the system.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig16} % lyapexpm.ps
\end{center}
\caption{Estimation of the maximal Lyapunov exponent for
$m \in \langle 0, 0.3 \rangle$, $a=0.1$, $r=0.25$, $v=0.031847$,
$\kappa=1$ and $p=0.01$. Initial conditions for all computed
trajectories: $x=0.45$, $y=0.8$ and $z=0.35$ for
$t \in \langle 3000, 6000 \rangle$} \label{fig16}
\end{figure}
\section{Poincar\'e section and a period three orbit}
Let us take a look at the Poincar\'e section. A Poincar\'e map
consists of a discrete set of values picked whenever one of the
variables passes through some prescribed value. We chose to plot
successive maxima of the variable $y$.
We computed local maxima $y(t_n)=y_n$ of the trajectory from the
Example 2. with initial conditions $x(0)=0.45$, $y(0)=0.8$ and
$z(0)=0.35$ for $t \in \langle 3000, 30000 \rangle$ and got a
sequence ${y_n}$. When we plotted this sequence to $y_n$ vs.
$y_{n+1}$ plane (Ruelle plot in Xppaut), we got a black curve on
the figure \ref{fig17}. The intersection with the first quadrant
axes is a fixed point $y_n=y_{n+1}$ ($FP$) - a period one orbit.
The sequence ${y_n}$ plotted in $y_n$ vs. $y_{n+3}$ plane, gave
another curve (blue). The intersections with the first quadrant
axes different from $FP$ are fixed point $y_n=y_{n+3}$ -
corresponding to period three orbits. According to \cite{li}, the
trajectory is chaotic. An example of a trajectory with a period
three $(x(0) \doteq 1.9821$, $y(0) \doteq 5.8385$,
$z(0) \doteq -1.9821)$ found due to this Poincar\'e mapping is on
figure \ref{fig18}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig17} % poincare.ps
\end{center}
\caption{Fixed point of the Poincar\'e mapping corresponding to a
period one orbit and intersections corresponding to period three
orbits.} \label{fig17}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig18} % 3period.ps
\end{center}
\caption{A period three orbit.} \label{fig18}
\end{figure}
Compared to computing maximal Lyapunov exponent, this evidence of
chaos can be used for single trajectories only and we have no
image of the parameter dependence. On the other hand, maximal
Lyapunov exponent estimation has an error, caused by finiteness of
computed points and even may break down during reorthogonalization
of the matrix $Q$ (while using the standard $QR$ method of
computing submitted by Benettin et al. in \cite{benettin}), since
the computed matrix $Q$ may deviate from the origin one during the
Gram-Schmidt orthogonalization type procedure. This last mentioned
problem has been studied and overcame by Udwadia and Bremen using
Cayley transformation (see \cite{udwadia}).
\section{0-1 test for chaos}
Recently, a new test for determining chaos was introduced by
Gottwald and Melbourne. In contrast to the usual method of
computing the maximal Lyapunov exponent, their method is applied
directly to the time series data and does not require phase space
reconstruction. We computed time series corresponding to
solutions $x$ of the system \ref{eq:system1} for parameters $a, \;
v, \; p$ and $m$ by the fourth-order classical Runge-Kutta method
with timestep 0.25 and initial conditions $x=0.45$, $y=0.8$ and
$z=0.35$. We used sampling time $\tau_s=12$ (using $\tau_s=5$ the
data are oversampled still) and got time series $x(t_j)$,
$t_j =j \tau_s$ for $j=1, \dots, 2000$. According to \cite{gottwald}, we
computed 100 times $p_c(n) = \sum_{j=1}^n x(t_j) \cos{jc}$ and
$q_c(n)=\sum_{j=1}^n x(t_j) \sin{jc}$, $j=1, \dots, 2000$ for
arbitrary chosen $c \in (0, \pi)$ and computed median $K$ of the
asymptotic growth rates $K_c$ of the mean-square-displacement
$$
M_c(n) = \lim_{N \to \infty} \frac1N \sum_{j=1}^N [p_c(j + n) -
p_c(j)]^2 + [q_c(j + n) - q_c(j)]^2.
$$
Since $N=2000$ and for the
estimation of the limit we need $n << N$, we used the last
$n_{cut}=200$. For computing $K_c$ we used the correlation method.
The median $K \approx 0$ indicates regular dynamics, $K \approx 1$
indicates chaos. Results are displayed on figures \ref{fig19},
\ref{fig20}, \ref{fig21} and \ref{fig22} and correspond to the
estimations of the maximal Lyapunov exponent on figures
\ref{fig10}, \ref{fig11}, \ref{fig12} and \ref{fig16}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig19} % 01testa.eps
\caption{0-1 test for chaos: chaos range for the parameter $a$,
$r=0.25$, $p=0.01$, $v=0.031847$, $\kappa=1$ and $m=0.19$.}
\label{fig19}
\end{center}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig20} % 01testv.eps
\caption{0-1 test for chaos: chaos range for the parameter $v$,
$a=0.1$, $r=0.25$, $p=0.01$, $\kappa=1$ and $m=0.19$.}
\label{fig20}
\end{center}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig21} % 01testp.eps
\caption{0-1 test for chaos: chaos range for the parameter $p$,
$a=0.1$, $r=0.25$, $v=0.031847$, $\kappa=1$ and $m=0.19$.}
\label{fig21}
\end{center}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics*[width=0.9\textwidth]{fig22} % 01testm.eps
\caption{0-1 test for chaos: chaos range for the parameter $m$,
$a=0.1$, $r=0.25$, $v=0.031847$, $\kappa=1$ and $p=0.01$.}
\label{fig22}
\end{center}
\end{figure}
\section{Hopf bifurcation near zero}
\label{near zero}
Now we will turn our attention to the zero neighbourhood. It was
already mentioned above, that for the trivial equilibrium $E_0$
(always unstable), the necessary condition \eqref{mpk} for Hopf
bifurcation is $m=-p\kappa$, while $\lambda=\pm i\omega$, where
$\omega^2= v(r-a)>0$, that is for $r>a$ (since $v>0$). The
economic system can never reach the critical value, since we
assumed the parameters to be positive. On the other hand, local
behaviour of the trajectories near this critical value may
coincide with the zero neighbourhood and affect the properties of
the positive quadrant. Due to this we let some necessary
parameters reach and even cross zero to negative values in the
next part.
Using the implicit function theorem we can compute the derivative
of the complex eigenvalue $\lambda$ according to $m$ for the
equilibrium $E_0$ (using \eqref{polyE0}) $\omega^2=v(r-a)$:
\begin{equation}
\frac{d\lambda}{dm}=-\frac{\frac{dp(\lambda)}{dm}}{\frac{dp(\lambda)}{d\lambda}}=
-\frac{-av}{3\lambda^2+2m\lambda+\omega^2}.
\end{equation}
Substituting $\lambda_{1,2}=\pm i\omega$, we evaluate
\begin{equation}
\label{dlambdadm}
\mathop{\rm sgn} \mathop{\rm Re} \frac{d\lambda}{dm}= \mathop{\rm sgn} \frac{-av}{2(\omega^2+m^2)}<0.
\end{equation}
The transversality condition for Hopf bifurcation is fulfilled and
according to \cite{kuznecov} the Hopf bifurcation gives rise to a
limit cycle near the trivial equilibrium at the critical value
$m=-p\kappa$ for $a0$, $v>0$ and
$p>0$ (subcritical Hopf bifurcation). Generalized Hopf bifurcation
of the trivial equilibrium occurs on the parametric manifold
\begin{equation}
\label{GH}
M_{GH}=\{(a, p, \kappa, v, m, r) \in \mathbb{R}^6: m=-p\kappa, a=0,
v>0, r>0, p>0\}
\end{equation}
\end{theorem}
\begin{proof}
We already proved the Hopf bifurcation occures for $m=-p \kappa$
for $a0$ and $r>a$ (the
opposite signs are irrelevant for the economic model). For an
invariant expression for the first Lyapunov coefficient we need to
find eigenvectors $Q$ and $P$ such that $JQ=i \omega Q$, $J^TP=-i
\omega P$ and $\langle P,Q \rangle = 1$. One of the eigenvectors
corresponding to $\lambda=i\omega$ is
$$
Q=\bigl (a, m+i\omega, -r+i\frac{\omega m}v \bigr ).
$$
The corresponding normalized eigenvector is $P=\frac1c \bigl
(1,-\frac{i\omega}v,1 \bigr )$, where
$$
c= \langle
(1,-\frac{i\omega}v,1 \bigr ),Q \rangle = 2(a-r) + 2i\frac{\omega
m}{v} \neq 0$$
(as in \cite{kuznecov} we use an unusual scalar
product definition $\langle x,y \rangle= \sum \overline{x_i} y_i$
to keep the Kuznetsov's notation). According to \cite{kuznecov}
$$
\mathop{\rm sgn} l_1= \mathop{\rm sgn} \frac1{2\omega}\mathop{\rm Re} \langle P, C(Q,Q,\bar{Q})\rangle,
$$
where for $i=1, \dots 3$ $F_i$ stays for nonlinear part of the
right hand side of \eqref{eq:system1} and we define
$$
C_i(x,y,z)=\sum_{j,k,l=1}^3 \frac{\partial^3 F_i(\xi)}
{\partial \xi_j \partial \xi_k \partial
\xi_l}|_{\xi=0}x_jy_kz_l,
$$
that is,
$C_1(x,y,z)=-4p(x_1y_2z_2+x_2y_1z_2+x_2y_2z_1)$, $C_2(x,y,z)=0$
and $C_3(x,y,z)=0$. (The general formula for $l_1$ contains also
quadratic members of the nonlinear part of the right hand side of
\eqref{eq:system1}, but it vanishes in our case). From that we get
$$
\mathop{\rm sgn} l_1=\mathop{\rm sgn} \frac{apv}{ \omega}.
$$
The generalized Hopf bifurcation occures for $l_1=0,$ that is only
for $a=0$. We have to exclude cases $p=0$, when the system
\eqref{eq:system1} is strictly linear, and of course $v=0$, when
$\lambda_{1,2}=\pm i \omega=0$. Consequently, $l_1>0$ for $r>a>0$,
$v>0$ and $p>0$ and subcritical Hopf bifurcation occurs near zero
for $m=-p \kappa$.
\end{proof}
As a consequence of \eqref{dlambdadm} and the Theorem \ref{l1} we
see, that unstable cycle arises on the right hand side of $m$,
that is for $m>-p\kappa$, since the real part of the complex
conjugate eigenvalues decreases from positive to negative values
as $m$ increases and an unstable equilibrium $(0,0,0)$ become
stable on the center manifold with an unstable cycle in its local
neighbourhood. This is the reason, why we included study of the
economically impossible Hopf bifurcation of the zero equilibrium
into account - an unstable limit cycle may continue (and already
does) in the positive parametric space. You can see the
bifurcation diagram on the figure \ref{fig23}
\begin{figure}[ht]
\begin{center}
\includegraphics*[width=0.7\textwidth]{fig23} % zero.ps
\end{center}
\caption{The typical Hopf bifurcation diagram with arising
unstable cycle near zero for parameters $a=0.1$,
$r=0.25$, $p=0.01$, $v=0.031847$, $\kappa=1$ and $m \in (-0.02,0.1)$.}
\label{fig23}
\end{figure}
\section{Conclusions}
In this paper, we studied several types of bifurcations in the
system \eqref{eq:system1} that are connected with periodic and
non-periodic bounded trajectories that may represent economic
cycle and non-periodic chaotic oscillations in macroeconomic
quantities. From the economic point of view, condition
\eqref{cond} does not guarantee existence of a stable equilibrium.
Nor even stronger condition $a > r$, that according to Theorem
\ref{ar} guarantees existence of a stable equilibria, does not
guarantee asymptotic tending towards one of the stable fixed
points, since for a set of parameters the first Lyapunov
coefficient is positive and folding of the unstable cycle into a
stable one give rise to a large basin of attraction of a stable
periodic solution for $a>r$ close to $r$. This result is in the
contrary to the conclusions made in \cite{vosvrda}. If the
debt/output ratio is less than the variation of the marginal
propensity to savings then the economic equilibrium is locally
stable, but the economy near a stable state may be on a non-local
destabilizing cycle asymptotically tending to a non-local stable
cycle. As we can see from results described by the figure
\ref{fig1}, for normal values of capital-output ratio $v << 1$,
the first Lyapunov coefficient is positive for quite big range of
parameter $a$, even the less is the value of capital-output ratio
$v$, the bigger "unstable" range we get for $a$. For normal
economic parameters $a>>r$ the equilibria $E_{1,2}$ are always
stable and attracting, no cycles or chaos appear in the system.
But for low $v$, $p$ or large $m$, there exist trajectories
corresponding to unstable trade cycles that may even change into
non-periodic bounded chaotic unpredictable regime. Due to this we
cannot agree with the conclusion in \cite{vosvrda} that if the
capital inflow/savings ratio is less than double the ratio of
capitalized profit then the system is in a stable state. You can
see on the figures \ref{fig2} and \ref{fig3} ($m=0.19$, $p=0.1$)
that for low variation of the marginal propensity to savings this
is not true, unstable cycles occurs. As it can be seen from
figures \ref{fig16} and \ref{fig22}, parameter $m$, that is
capital inflow/savings ratio, cannot be used for controlling the
system at all. Similarly $p$, the ratio of the capitalized profit,
cannot be used for controlling (see figures \ref{fig12},
\ref{fig15} and \ref{fig21}). The first conclusion in
\cite{vosvrda} is true: an increasing of the capitalization of the
profits demonstrates well-known results in economics that the
capitalization of profits causes the stabilization of the system.
It should be taken into account, that the model of foreign
financing is highly simplified and therefore the conclusions may
represent possibilities of the real economic behaviour only. On
the other hand the model is "nearly linear" - and linear models
are often used in economics to show "the invisible hand of the
market" that leads the economy to the stable and predictable,
equalized quantities. The one cubic term included in the first
equation (also included in the model without foreign financing
with no bifurcations at all, for analysis of this model see
\cite{vosvrda}) give rise to a great deal of nonlinear dynamics -
to periodic and chaotic motion with no invisible hand to lead.
The mathematical results may be applied to the second model of
firm introduced by \cite{bouali} that is based on the Hunt's
hypothesis that the call for loans pushes the profit ratio of
stockholders’capital. Bouali shows various periodicity and chaos
in the system without deeper mathematical background. The results
derived here explain the dynamics in this model of firm. In order
to illustrate farther economic applicability of these results, we
cite the economic conclusions made by Bouali. ``\dots rules and
principles of finance governance built in static framework may
lose their validity. The findings of a well corporate debt policy
connected to a well dividend policy may lead to an unpredictable
and hazardous motion of the profit. In our 3D system, the rise of
the loss level is an endogenous outcome of the borrowing policy
and is not determined by a shock of economic recession. Against
the common sense, the profit motion is worsened by the braking of
dividend distribution!''
The existence of hysteresis and various types of bifurcations and
chaos open a question, whether the stabilization policy is
efficient. The policy advice is based mostly on linear models
while the economy is actually characterized by significant
nonlinearities. Linear modelling of a system with existence of
significant nonlinearity in the data may provide misleading
results. In the case of hysteresis, stabilization policy may lead
to destabilizing trade cycles instead of tending to an equilibrium
or to chaos. In the case of various cycle bifurcations the data
itself cannot be correctly analyzed by methods based on linear
modelling and in the case of chaos, the precision of a forecast
with a very small error in the initial conditions worsens
exponentially over time in addition to this previously mentioned
failings.
\begin{thebibliography}{99}
\bibitem{alligood} Alligood, K.T., Sauer, T.D., Yorke, J.A. (1996)
{\it Chaos -
an Introduction to Dynamical Systems}, Berlin, Heidelgerg, New
York, Springer-Verlag
\bibitem{benettin} Benettin, B., Galgani, L., Giorgilli, A. and
Strelcyn, J. M. (1980) { Lyapunov Characteristic Exponents for
Smooth Dynamical Systems and Hamiltonian Systems; a Method for
Computing All of Them, Part I. and Part II.}, {\it Meccanicca} 15,
9-30
\bibitem{bouali} Bouali, S. (2002)
{The Hunt Hypothesis and the Dividend Policy
of the Firm. The Chaotic Motion of the Profits.}, {\it 8th
International Conference of the Society for Computational
Economics Computing in Economics and Finance}, Aix-en-Provence,
France, June 27-29
\bibitem{dadda} D'adda, C., Scorcu, A.E. (2003) {On the time stability
of the output-capital ratio}, {\it Economic Modelling} 20, issue
6, 1175-1189
\bibitem{gottwald} Gottwald, A G., Melbourne, I. (2009)
{On the Implementation of the 0–1 Test for Chaos},
{\it SIAM J. Appl. Dyn. Syst.} 8, issue 1, pp. 129-145
\bibitem{kuznecov} Kuznetsov, Y.A. (1998) {\it Elements of Applied
Bifucation Theory}, Second Edition, Applied Mathematical Sciences
112, Berlin, Heidelgerg, New York, Springer-Verlag
\bibitem{li} Li, T. Y., Yorke, J.A. (1975) {Period Three Implies
Chaos}, {\it American Mathematical Monthly} 82, 985-992
\bibitem{ruelle} Eckmann, J. P., Ruelle, D. (1985)
{Ergodic Theory of Chaos and Strange
Attractors}, {\it Rev. Mod. Phys.} 57, 617-656.
\bibitem{vosvrda} Vo\v svrda, M. (2001) {Bifurcation Routes and Economic
Stability}, {\it Bulletin of the Czech Econometric Society} 14,
43-60
\bibitem{udwadia} Udwadia, F. E., von Bremen H. F. (2001) {Computation
of Lyapunov Characteristic Exponents for Continuous Dynamical
Systems}, {\it Math. Phys.}, Vol. 53, 23-146
\bibitem{wolf} Wolf A., Swift, J. B., Swinney, H. L. and Vastano, J.
A. (1985) {Determining Lyapunov Exponents from a Time Series},
{\it Physica D} 16, 285-317.
\end{thebibliography}
\end{document}