\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}