\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 209, pp. 1--11.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/209\hfil Computation of focal values]
{Computation of focal values and stability analysis
 of 4-dimensional systems}

\author[B. Sang,  Q.-L. Wang, W.-T. Huang \hfil EJDE-2015/209\hfilneg]
{Bo Sang, Qin-Long Wang, Wen-Tao Huang}

\address{Bo Sang (corresponding author) \newline
School of Mathematical Sciences, Liaocheng University,
Liaocheng 252059, China}
\email{sangbo\_76@163.com}

\address{Qing-Long Wang \newline
School  of Science, Hezhou University, Hezhou 542800, China}
\email{wqinlong@163.com}

\address{Wen-Tao Huang \newline
School  of Science, Hezhou University, Hezhou 542800, China}
\email{huangwentao@163.com}

\thanks{Submitted July 1, 2015. Published August 10, 2015.}
\subjclass[2010]{34C05, 34C07}
\keywords{Focal value; limit cycle; Hopf bifurcation}

\begin{abstract}
 This article presents a recursive formula for computing the $n$-th
 singular point values of a class of 4-dimensional  autonomous systems,
 and establishes the algebraic equivalence between focal values and singular
 point values. The formula is linear and then avoids complicated integrating
 operations, therefore the calculation can be carried out by computer algebra
 system such as Maple. As an application of the formula, bifurcation analysis
 is made for a quadratic system with a Hopf equilibrium, which can have three
 small limit cycles around an equilibrium point. The theory and  methodology
 developed in this paper can be used for higher-dimensional systems.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\allowdisplaybreaks

\section{Preliminaries}

Consider a $C^{k}$-smooth system ($k\geq{2}$)
\begin{equation}\label{sys1}
\begin{gathered}
\frac{dx}{dt}=Ax+f(x,y),\\
\frac{dy}{dt}=By+g(x,y),
\end{gathered}
\end{equation}
where $x\in \mathbb{R}^{n_{c}}, y\in \mathbb{R}^{n_{s}}$, $A$ and $B$
are constant matrices, and $f(x,y),g(x,y)$ are functions with
$$
f(0,0)=0,g(0,0)=0, Df(0,0)=0, Dg(0,0)=0.
$$
Suppose that $A$ has $n_{c}$ critical eigenvalues (i.e. eigenvalues with
$\operatorname{Re} \lambda$ = 0) and all $n_{s}$ eigenvalues of $B$ satisfy
$\operatorname{Re} \lambda< 0$.
According to the Center Manifold Theorem (see e.g. \cite{k04}), there exists
a (local) center manifold
$y=h(x)$ with  $h(0) = 0, Dh(0) = 0$,
and system \eqref{sys1} is  topologically equivalent near $(0, 0)$ to the system
\begin{equation}\label{sys2}
\begin{aligned}
 \frac{dx}{dt}&=Ax+f(x,h(x)),\\
\frac{dy}{dt}&=By.
\end{aligned}
\end{equation}
The first equation in \eqref{sys2} is called the restriction of system
\eqref{sys1} to its center manifold at the origin. Thus,
the dynamics of \eqref{sys1} near a non-hyperbolic equilibrium are determined
by this restriction, since the second equation in \eqref{sys2} is linear
and has exponentially decaying solutions.

If $A$ has a simple pair of purely imaginary eigenvalues
$\pm \omega{i}$ ($\omega>0$), system \eqref{sys1} undergoes a Hopf bifurcation
or a multiple Hopf bifurcation under proper perturbations of parameters.
The computation of focal values (Lyapunov coefficients) plays an important
 role in the study of small-amplitude limit cycles appeared in these
bifurcations (see \cite{bvl,gy09,han15,lv,mrs,t14,yu13,w10}
and reference therein). For system \eqref{sys1},  the projection method was
used for computing the first and the second focal values (see \cite{k04});
 a perturbation technique based on multiple time scales was used for computing
focal values (see \cite{yu98}). For a class of 3-dimensional systems,
a recursive formula was presented  for computing focal values (see \cite{w10}).
It should be noted that the theory  and methodology described in \cite{s14,w10}
can be applied to $N$-dimensional systems, where $N\geq {4}$.

\section{Focal values of a class of 4-dimensional systems}

Consider a class of analytic systems in $\mathbb{R}^4$,
\begin{equation}\label{sys3}
\begin{aligned} 
\frac{{dx_1}}{dt}
&=-x_2+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
\tilde{a}_{j_1,j_2,j_3,j_4}x_1^{j_1}x_2^{j_2}x_3^{j_3}x_4^{j_4}
=X_1(x_1,x_2,x_3,x_4),\\
\frac{{dx_2}}{dt}
&=x_1+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}\tilde{b}_{j_1,j_2,j_3,j_4}x_1^{j_1}
 x_2^{j_2}x_3^{j_3}x_4^{j_4}
=X_2(x_1,x_2,x_3,x_4),\\
\frac{{dx_3}}{dt}
&=\mu x_3-\omega x_4+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
\tilde{c}_{j_1,j_2,j_3,j_4}x_1^{j_1}x_2^{j_2}x_3^{j_3}x_4^{j_4}\\
&=X_3(x_1,x_2,x_3,x_4),\\
\frac{{dx_4}}{dt}
&=\omega x_3+\mu x_4+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
\tilde{d}_{j_1,j_2,j_3,j_4}x_1^{j_1}x_2^{j_2}x_3^{j_3}x_4^{j_4}\\
&=X_4(x_1,x_2,x_3,x_4),
\end{aligned}
\end{equation}
with $\mu<0$, $\omega\geq {0}$,
$x_1,x_2,x_3,x_4,\tilde{a}_{j_1j_2j_3j_4},\tilde{b}_{j_1j_2j_3j_4},
\tilde{c}_{j_1j_2j_3j_4},
\tilde{d}_{j_1j_2j_3j_4}$ in $\mathbb{R}$, and
$j_1$, $j_2$, $j_3$, $j_4$ in $\mathbb{N}$. 
Our motivation for the
study of system \eqref{sys3} is that a physical model of airfoil
with cubic nonlinearity can be transformed into a special case of
system \eqref{sys3} via a linear transformation, see \cite{l92,z04}.

The Jacobian matrix of system \eqref{sys3} has eigenvalues
$\pm {i}, \mu\pm \omega{i}$
at the equilibrium $(x_1,x_2,x_3,x_4)=(0,0,0,0)$.
Then by the Center Manifold Theorem  system \eqref{sys3} has a (local)
center manifold tangent to the $(x_1, x_2)$ plane at the origin.
Moreover this center manifold can be represented as
\begin{equation}\label{eqn}
\begin{aligned}
x_3=g_1(x_1,x_2), \\
x_4=g_2(x_1,x_2),
\end{aligned}
\end{equation}
where
$g_1(0,0)=g_2(0,0)=0, Dg_1(0,0)=Dg_2(0,0)=0$, and the dynamics
of \eqref{sys3} restricted to the center manifold are given by
\begin{equation}\label{sys4}
\begin{gathered}
\frac{{dx_1}}{dt}
=-x_2+\sum_{j_1+j_2=2}^{\infty}\tilde{a}_{j_1,j_2}
 x_1^{j_1}x_2^{j_2}=X_1(x_1,x_2,g_1(x_1,x_2),g_2(x_1,x_2)),\\
\frac{{dx_2}}{dt}=x_1+\sum_{j_1+j_2=2}^{\infty}\tilde{b}_{j_1,j_2}
 x_1^{j_1}x_2^{j_2}=X_2(x_1,x_2,g_1(x_1,x_2),g_2(x_1,x_2)),
\end{gathered}
\end{equation}
with
\begin{gather*}
  \frac{d{g_1}(x_1,x_2)}{dt}\Big|_{\eqref{sys4}}
=X_3(x_1,x_2,g_1(x_1,x_2),g_2(x_1,x_2)),\\
\frac{d{g_2(x_1,x_2)}}{dt}\Big|_{\eqref{sys4}}
=X_4(x_1,x_2,g_1(x_1,x_2),g_2(x_1,x_2)).
\end{gather*}
From \cite{l08}, we know that for system \eqref{sys4} one can
derive successively and uniquely the terms of the following formal
series
\begin{equation}
\tilde{F}(x_1,x_2)=x_1^2+x_2^2+h.o.t.
\end{equation}
such that
\begin{equation}\label{stab1}
\begin{aligned}
\frac{d{\tilde{F}}}{dt}\Big|_{\eqref{sys4}}
&=\frac{{\partial \tilde{F}}}{{\partial x_1}}X_1(x_1,x_2,
g_1(x_1,x_2),g_2(x_1,x_2))\\
&\quad +\frac{{\partial \tilde{F}}}{{\partial x_2}}
X_2(x_1,x_2,g_1(x_1,x_2),g_2(x_1,x_2)),\\
&=\sum_{n=1}^{\infty} V_{n}(x_1^2+x_2^2)^{n+1},
\end{aligned}
\end{equation}
where $V_{n}$ are called the $n$th focal values of system
\eqref{sys4} and the original system \eqref{sys3}, and the acronym
h.o.t. stands for higher-order terms.


By change of variables
\begin{equation}\label{trans}
\begin{gathered}
x_1=\frac{1}{2}(x+y),\quad x_2=-\frac{i}{2}(x-y), \quad
x_3=\frac{1}{2}(z+u), \\
 x_4=-\frac{i}{2}(z-u), \quad t=-{i}T,
\end{gathered}
\end{equation}
the original system \eqref{sys3} can be transformed into the following complex system
\begin{equation}\label{sys5}
\begin{aligned} 
\frac{dx}{{dT}}&=x+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
 {a}_{j_1,j_2,j_3,j_4}x^{j_1}y^{j_2}z^{j_3}u^{j_4}=X(x,y,z,u),\\
\frac{dy}{{dT}}&=-y+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}{b}_{j_1,j_2,j_3,j_4}
 x^{j_1}y^{j_2}z^{j_3}u^{j_4}=Y(x,y,z,u),\\
\frac{{dz}}{{dT}}&=(\omega-\mu i)z+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
{c}_{j_1,j_2,j_3,j_4}x^{j_1}y^{j_2}z^{j_3}u^{j_4}=Z(x,y,z,u),\\
\frac{{du}}{{dT}}&=-(\omega+\mu i)u+\sum_{j_1+j_2+j_3+j_4=2}^{\infty}
{d}_{j_1,j_2,j_3,j_4}x^{j_1}y^{j_2}z^{j_3}u^{j_4}=U(x,y,z,u).
\end{aligned} 
\end{equation}
And, applying the transformation
\begin{equation}\label{ntrans}
x_1=\frac{1}{2}(x+y), x_2=-\frac{i}{2}(x-y), t=-{i}T,
\end{equation}
the restriction system \eqref{sys4} can be transformed into the following
complex system
\begin{equation}\label{res10}
\begin{gathered}
\frac{dx}{{dT}}=x+\sum_{j_1+j_2=2}^{\infty}{a}_{j_1,j_2}x^{j_1}y^{j_2}
 =X(x,y,z(x,y),u(x,y)),\\
\frac{dy}{{dT}}=-y+\sum_{j_1+j_2=2}^{\infty}{b}_{j_1,j_2}x^{j_1}y^{j_2}
 =Y(x,y,z(x,y),u(x,y)),
\end{gathered}
\end{equation}
where
\begin{gather}
\begin{aligned}
z(x,y)&= x_3+ix_4\\
      &= g_1(x_1,x_2)+ig_2(x_1,x_2)\\
      &= g_1(\frac{x+y}{2},-\frac{{i}}{2}(x-y))+ig_2(\frac{x+y}{2},
-\frac{i}{2}(x-y)),
\end{aligned} \label{cond11}\\
\begin{aligned}
u(x,y)&= x_3-ix_4\\
      &= g_1(x_1,x_2)-ig_2(x_1,x_2)\\
      &= g_1(\frac{x+y}{2},-\frac{{i}}{2}(x-y))-ig_2(\frac{x+y}{2},
-\frac{i}{2}(x-y)).
\end{aligned}\label{cond12}
\end{gather}
with
\begin{equation}\label{in13}
\begin{gathered}
\frac{dz(x,y)}{dT}\Big|_{\eqref{res10}}=Z(x,y,z(x,y),u(x,y)),\\
\frac{du(x,y)}{dT}\Big|_{\eqref{res10}}=U(x,y,z(x,y),u(x,y)).
\end{gathered}
\end{equation}

\begin{lemma}[\cite{l01}] \label{lem1}
For system \eqref{res10}, one can derive successively and uniquely the terms
 of the formal series
\[
G(x,y)= xy+\sum_{\alpha+\beta=3}^{\infty}{C_{\alpha,\beta}}x^{\alpha}y^{\beta},\\
\]
such that
\begin{equation}\label{eq10}
\begin{aligned}
\frac{dG}{dT}\Big|_{\eqref{res10}}
&=\frac{{\partial G}}{{\partial x}}{X(x,y,z(x,y),u(x,y))} +
\frac{{\partial G}}{{\partial y}}{Y(x,y,z(x,y),u(x,y))}\\
&=\sum_{n=1}^{\infty}W_{n}^{(2)}(xy)^{n+1},
\end{aligned}
\end{equation}
where ${C_{\alpha,\beta}}$ are determined by the recursive formula
(see \cite {l01}), and
$W_{n}^{(2)}$ are called the $n$th singular point values of system \eqref{res10}.
\end{lemma}

\begin{lemma}[\cite{l01}] \label{lem2}
If $W_1^{(2)}=W_2^{(2)}=\dots=W_{n-1}^{(2)}=0$, we have
\begin{equation}\label{res}
V_{n}=iW_{n}^{(2)},
\end{equation}
where $V_{n}$ is the $n$th focal value of system \eqref{sys4}, and
$W_{j}^{(2)}$ are the  $j$th singular point values of system  \eqref{res10},
 $j=1,2,\dots,n-1,n$.
\end{lemma}

\begin{theorem}\label{thm1}
For system \eqref{sys5}, we can derive successively and uniquely the terms
of the formal series
\begin{align*}
F(x,y,z,u)
&= xy+\sum_{\alpha+\beta+\gamma+\delta=3}^{\infty}
{C_{\alpha,\beta,\gamma,\delta}}x^{\alpha}y^{\beta}z^{\gamma}u^{\delta},\\
&\overset{\bigtriangleup}{=} \sum_{\alpha+\beta+\gamma+\delta=2}^{\infty}
{C_{\alpha,\beta,\gamma,\delta}}x^{\alpha}y^{\beta}z^{\gamma}u^{\delta},
\end{align*}
such that
\begin{equation}\label{eq10b}
\frac{dF}{dT}\Big|_{\eqref{sys5}}=\frac{{\partial F}}{{\partial x}}{X}
+\frac{{\partial F}}{{\partial y}}{Y}+\frac{{\partial F}}{{\partial z}}{Z}
 +\frac{{\partial F}}{{\partial u}}{U}=\sum_{n=1}^{\infty}W_{n}(xy)^{n+1},
\end{equation}
where ${C_{\alpha,\beta,\gamma,\delta}}$ are determined by the recursive formula
\begin{equation}\label{eq11}
\begin{aligned}
C_{\alpha,\beta,\gamma,\delta}
&= \frac{1}{\beta-\alpha-\gamma(\omega-\mu i)+\delta(\omega+\mu i)}\\
&\quad\times \sum_{j_1+j_2+j_3+j_4={3}}^{\alpha+\beta+\gamma+\delta+2}
\Big[(\alpha-j_1+1)a_{j_1,j_2-1,j_3,j_4}
+(\beta-j_2+1)b_{j_1-1,j_2,j_3,j_4}\\
&\quad +(\gamma-j_3)c_{j_1-1,j_2-1,j_3+1,j_4}
 +(\delta-j_4)d_{j_1-1,j_2-1,j_3,j_4+1}\Big]\\
&\quad\times C_{\alpha-j_1+1,\beta-j_2+1, \gamma-j_3,\delta-j_4}
\end{aligned}
\end{equation}
with ${C_{\alpha,\alpha,0,0}}=0$, $W_{n}$ are determined by the recursive formula
\begin{equation}\label{eq12}
W_{n}=\sum_{j_1+j_2=3}^{2n+4}[(n-j_1+2)a_{j_1,j_2-1,0,0}
+(n-j_2+2)b_{j_1-1,j_2,0,0}]C_{n-j_1+2,n-j_2+2,0,0}.
\end{equation}
\end{theorem}

\begin{proof}
Note that
\begin{align*}
\frac{dF}{dT}\Big|_{\eqref{sys5}}
&= \frac{{\partial F}}{{\partial x}}{X} 
 +\frac{{\partial F}}{{\partial y}}{Y}+\frac{{\partial F}}{{\partial z}}{Z}
 +\frac{{\partial F}}{{\partial u}}{U}\\
&= \sum_{\alpha+\beta+\gamma+\delta\geq {2}}{[\alpha-\beta
 +\gamma(\omega-\mu i)-\delta(\omega+\mu i)
 ]C_{\alpha,\beta,\gamma,\delta}}x^{\alpha}y^{\beta}z^{\gamma}u^{\delta}\\
&\quad +\sum_{\alpha+\beta+\gamma+\delta\geq {2}}
 \sum_{j_1+j_2+j_3+j_4\geq {3}} 
\Big(\alpha a_{j_1,j_2-1,j_3,j_4}+\beta b_{j_1-1,j_2,j_3,j_4}\\
&\quad +\gamma c_{j_1-1,j_2-1,j_3+1,j_4}
 +\delta d_{j_1-1,j_2-1,j_3,j_4+1}\Big)\\
&\quad\times  C_{\alpha,\beta,\gamma,\delta}x^{\alpha+j_1-1}y^{\beta+j_2-1}
 z^{\gamma+j_3}u^{\delta+j_4}\\
&= \sum_{\alpha+\beta+\gamma+\delta \geq{2}}x^{\alpha}y^{\beta}z^{\gamma}u^{\delta}
\Big\{[\alpha-\beta+\gamma(\omega-\mu i)-\delta(\omega+\mu i)]
 C_{\alpha,\beta,\gamma,\delta}\\
&\quad +\sum_{j_1+j_2+j_3+j_4\geq {3}}\Big[(\alpha-j_1+1)a_{j_1,j_2-1,j_3,j_4}
 +(\beta-j_2+1)b_{j_1-1,j_2,j_3,j_4}\\
&\quad +(\gamma-j_3)c_{j_1-1,j_2-1,j_3+1,j_4}
+(\delta-j_4)d_{j_1-1,j_2-1,j_3,j_4+1}\Big]\\
&\quad\times C_{\alpha-j_1+1,\beta-j_2+1, \gamma-j_3,\delta-j_4}\Big\}.
\end{align*}
Comparing the above power series with the right side of \eqref{eq10b},
we can obtain the recursive
formulas \eqref{eq11} and \eqref{eq12}.
\end{proof}

\begin{theorem}\label{thm2}
If $W_1=W_2=\dots=W_{n-1}=0$, then
\begin{equation}\label{res2}
W_{n}^{(2)}=W_{n},
\end{equation}
where $W_{j}$ are the  $j$th singular point values of system  \eqref{sys5}, 
and $W_{n}^{(2)}$ is the $n$th singular point value of system\eqref{res10}.
\end{theorem}

\begin{proof}
Suppose that $W_1=W_2=\dots=W_{n-1}=0$, there exists a unique polynomial
\begin{equation}\label{eq15}
F_{2n+2}(x,y,z,u)=xy+\sum_{\alpha+\beta+\gamma+\delta=3}^{2n+2}
{C_{\alpha,\beta,\gamma,\delta}}x^{\alpha}y^{\beta}z^{\gamma}u^{\delta},
\end{equation}
such that
\begin{equation}\label{eq16}
\frac{dF_{2n+2}}{dT}\Big|_{\eqref{sys5}}
=\frac{{\partial F_{2n+2}}}{{\partial x}}{X} +
\frac{{\partial F_{2n+2}}}{{\partial y}}{Y}+\frac{{\partial F_{2n+2}}}{{\partial z}}{Z}
+\frac{{\partial F_{2n+2}}}{{\partial u}}{U}
=W_{n}(xy)^{n+1}+\dots.
\end{equation}
Let
\begin{equation}\label{eq17}
G_{2n+2}(x,y)=F_{2n+2}(x,y,z(x,y),u(x,y)),
\end{equation}
where $z(x,y),u(x,y)$ are determined by conditions \eqref{cond11} 
and \eqref{cond12}.
We have
\begin{align*} %\label{eq18}
\frac{dG_{2n+2}}{dT}\Big|_{\eqref{res10}}
&= \frac{{\partial G_{2n+2}}}{{\partial x}}{X(x,y,z(x,y),u(x,y))} +
\frac{{\partial G_{2n+2}}}{{\partial y}}{Y(x,y,z(x,y),u(x,y))}\\
&= \Big(\frac{{\partial F_{2n+2}}}{{\partial x}} 
 +\frac{{\partial F_{2n+2}}}{{\partial z}}\frac{\partial{z}}{\partial x}
+\frac{{\partial F_{2n+2}}}{{\partial u}}\frac{\partial{u}}{\partial x}
 \Big){X(x,y,z(x,y),u(x,y))}\\
&\quad +\Big(\frac{{\partial F_{2n+2}}}{{\partial y}} +\frac{{\partial F_{2n+2}}}
 {{\partial z}}\frac{\partial{z}}{\partial y}
+\frac{{\partial F_{2n+2}}}{{\partial u}}\frac{\partial{u}}{\partial y}\Big)
 {Y(x,y,z(x,y),u(x,y))}\\
&= \frac{{\partial F_{2n+2}}}{{\partial x}}{X(x,y,z(x,y),u(x,y))} 
 +\frac{{\partial F_{2n+2}}}{{\partial y}}{Y(x,y,z(x,y),u(x,y))}\\
&\quad +\frac{{\partial F_{2n+2}}}{{\partial z}}\Big(\frac{\partial{z}}{\partial x}
 X(x,y,z(x,y),u(x,y))
+\frac{\partial{z}}{\partial y}Y(x,y,z(x,y),u(x,y))\Big)\\
&\quad +\frac{{\partial F_{2n+2}}}{{\partial u}}\Big(\frac{\partial{u}}{\partial x}
 X(x,y,z(x,y),u(x,y))
+\frac{\partial{u}}{\partial y}Y(x,y,z(x,y),u(x,y))\Big)\\
&= \frac{{\partial F_{2n+2}}}{{\partial x}}{X(x,y,z(x,y),u(x,y))}
  +\frac{{\partial F_{2n+2}}}{{\partial y}}{Y(x,y,z(x,y),u(x,y))}\\
&\quad +\frac{{\partial F_{2n+2}}}{{\partial z}}\frac{dz(x,y)}{dT}
 \Big|_{\eqref{res10}}++\frac{{\partial F_{2n+2}}}{{\partial u}}
 \frac{du(x,y)}{dT}\Big|_{\eqref{res10}}\\
&= \frac{{\partial F_{2n+2}}}{{\partial x}}{X(x,y,z(x,y),u(x,y))} 
 +\frac{{\partial F_{2n+2}}}{{\partial y}}{Y(x,y,z(x,y),u(x,y))}\\
&\quad +\frac{{\partial F_{2n+2}}}{{\partial z}}Z(x,y,z(x,y),u(x,y))
 +\frac{{\partial F_{2n+2}}}{{\partial u}}U(x,y,z(x,y),u(x,y))\\
&= \frac{dF_{2n+2}}{dT}\Big|_{\eqref{sys5},z=z(x,y),u=u(x,y)}\\
&= W_{n}(xy)^{n+1}+\dots.
\end{align*}
Hence, the $n$th singular point value $W_{n}^{(2)}$ of system \eqref{res10}  
satisfies $W_{n}^{(2)}=W_{n}$.
\end{proof}

Summarizing the results of Lemma \ref{lem2} and Theorem \ref{thm2} 
we obtain the following Corollary.

\begin{corollary}\label{coro1}
If $W_1=W_2=\dots=W_{n-1}=0$, then
\begin{equation}\label{res3}
V_{n}=iW_{n},
\end{equation}
where $V_{n}$ is the $n$th focal value of system \eqref{sys3}, 
and $W_{j}$ are the  $j$th singular point values of system  \eqref{sys5}, 
$j=1,2,\dots,n-1,n$.
\end{corollary}

\section{Stability analysis and degenerate Hopf bifurcation of 4-dimensional 
systems}

As an application of the results obtained in Section 2,  we consider a class 
of 4-dimensional quadratic systems
\begin{equation}\label{sys27} 
\begin{gathered} 
\frac{{dx_1}}{dt}=-x_2-a_1x_1x_{{3}}+a_1{x_1}^2,\\
\frac{{dx_2}}{dt}=x_1+a_2{x_2}^2+a_{{3}}x_1x_4,\\
\frac{{dx_3}}{dt}=-x_{{3}}-3\,x_4+a_2x_1x_2,\\
\frac{{dx_4}}{dt}=3\,x_{{3}}-x_4+a_4x_2x_{{3}}.
\end{gathered}
\end{equation}
Applying the recursive formulas \eqref{eq11}, \eqref{eq12} and the 
relation \eqref{res3}, we obtain the first four focal values 
$V_1, V_2, V_3, V_4$ of system \eqref{sys27}, where
\begin{gather*}
V_1= -{\frac {a_2 ( -9\,a_{{3}}+4\,a_1) }{104}},
\\
\begin{aligned}
V_2&= \Big(a_2a_{{3}} \big( 16288820\,{a_2}^2-609030\,a_2a
_{{3}}-342086\,a_2a_4-95543550\,{a_{{3}}}^2\\
&\quad -890367\,a_{{3}}a
_4+235464\,{a_4}^2 \big) \Big)/ 26166400,\\
\end{aligned}
\\
\begin{aligned}
V_3&= -\frac {a_2a_1}{509581408485615245796556800}
 \Big( 224716160498681163828769344{a_1}^4 \\
&\quad +23391560818512386519330256a_2{a_1}^{3}\\
&\quad +108275873757375918099514632a_4{a_1}^{3} \\
&\quad -519476645118198963455620296{a_2}^2{a_1}^2\\
&\quad +30897799320198541896741988{a_1}^2a_2a_4 \\
&\quad -2608855859269938585435348{a_2}^{3}a_1\\
&\quad -107827566965381888322448434a_1{a_2}^2a_4\\
&\quad +303786319994459587929292974{a_2}^4\\
&\quad -22649516190250502017478481{a_2}^{3}a_4 \Big),
\end{aligned}\\  
\begin{aligned}
V_4&= \big({a_1a_2^3}/123387308399787108515183370468179366315299358650342667\\
&\quad 6059415722103342538117098506095377725888053245618196607199839\\
&\quad 4040115200000\big) V_{4,1}
\end{aligned}
\end{gather*}
and 
\begin{align*}
V_{4,1}
&= 151090339549302141420998778418569305630254453586345260252801201492\\
&\quad 19567253160806434635187840970326748821238775751413433121260736a_1^4\\
&\quad + 6902376529367432579800252656210544307 8560901105903707085784563088\\
&\quad 5126271156862905782616567832231124991999486811410920350778464
 a_2{a_1}^{3}\\
&\quad - 4852869311179757471922805489901407470548086994648890234025102002\\
&\quad 3568501430203732065339603894449313619014039201735819511942218584\\
&\quad\times a_2^2{a_1}^2\\
&\quad + 221525274767927679154097633424938981478458627350721696040337613\\
&\quad 53821 383571195398731665935703556310724308392281864424192874143532\\
&\quad\times a_1^2a_2a_4\\
&\quad + 434975337903447724813720015538420471541695549645684668496869670\\
&\quad 2641954333545244084235631897431187963845500346767147747172307568\\
&\quad\times {a_2}^{3}a_1\\
&\quad +139925156878486203251650275215133948778485624374385690191244381\\
&\quad 7572961470824125252612724121465334679579965759273450126018896624\\
&\quad\times a_1{a_2}^2a_4\\
&\quad + 34130451364797077586851329711670405855176881838233683394612381\\
&\quad 837487804749465866800870958247599456165586110349460109943436651626\\
&\quad \times a_2^4\\
&\quad -21736399966566178674551590653380371038422117740968675610514295\\
&\quad 706203182388856113401624201196405347456215992463433041450052098371\\
&\quad \times a_2^3a_4.
\end{align*}
Here $V_{j}$ is reduced modulo the
Gr\"{o}bner basis of $\{V_{s}: 1 \leq s \leq j-1\}$ for $j\geq 2$.

From \eqref{stab1} and the computation of focal values, we get the following
result on the stability of equilibrium for system \eqref{sys27}.

\begin{theorem} \label{thm3}
Let $\mathfrak{X}$ be the vector field determined on $\mathbb{R}^4$
by system \eqref{sys27}.
For any center manifold $W^{c}$ of \eqref{sys27} at the origin of 
$\mathbb{R}^4$, with regard to $\mathfrak{X}|W^{c}$:
\begin{enumerate}
\item  if $V_1\neq {0}$ then the origin is a first order fine focus,
 whose stability is determined by $\operatorname{sgn}V_1$ (i.e., is asymptotically stable 
if and only if $V_1 < 0$);

\item  if $V_1={0}, V_2 \neq {0}$ then the origin is a second order fine focus,
 whose stability is determined by $\operatorname{sgn}V_2$;

\item  if $V_1=V_2=0, V_3 \neq {0}$ then the origin is a third order fine focus, 
whose stability is determined by $\operatorname{sgn}V_3$;

\item  if $V_1=V_2=V_3=0, V_4 \neq {0}$ then the origin is a fourth order 
fine focus, whose stability is determined by $\operatorname{sgn}V_4$ .
\end{enumerate}
\end{theorem}


\begin{theorem} \label{thm4}
Suppose that 
\begin{equation}\label{cond}
\begin{gathered}
a_1=  \mu\,a_2,\\
a_{{3}}= {\frac{4}{9}}\mu \,a_2,\\
a_4= \frac{F_1(a_2)}{F_2},
\end{gathered}
\end{equation}
where $\mu$ is one of four real roots of the polynomial
\begin{equation}
\begin{aligned}
P(\lambda)
&= 10799985775088040041585439500413156070699520{\lambda}^{8}\\
&\quad +6212790925817703957911789339676019539872928{\lambda}^{7}\\
&\quad -28584647405659759515776062165338760876812544{\lambda}^{6}\\
&\quad -15713629519470248682246891010001873803911248{\lambda}^{5}\\
&\quad +23373024697622801653163411830546195419879000{\lambda}^4\\
&\quad +13015318562221586811114336007971951317356918{\lambda}^{3}\\
&\quad -4091058129972655937988952508114319926195742{\lambda}^2\\
&\quad -3522262300428215202057982497109271042377107\lambda\\
&\quad -1499503145083805883162951654763171839112800,
\end{aligned}
\end{equation}
and
\begin{gather}
\begin{aligned}
F_1(a_2)
&= -6a_2 \Big( 37452693416446860638128224{\mu}^4+
3898593469752064419888376{\mu}^{3} \\
&\quad -86579440853033160575936716{\mu}
^2-434809309878323097572558\mu \\
&\quad +50631053332409931321548829 \Big),
\end{aligned}\\
\begin{aligned}
F_2&= 108275873757375918099514632{\mu}^{3}+30897799320198541896741988{\mu}^2\\
&\quad -107827566965381888322448434\mu-22649516190250502017478481,\\
\end{aligned}\\
a_2 \neq {0},
\end{gather}
then the origin of system \eqref{sys27} is a fourth order fine focus.
\end{theorem}

\begin{proof}
Let condition \eqref{cond} be satisfied, then the first four focal values 
of system \eqref{sys27} are as follows:
$V_1= V_2=V_3=0$,
\begin{align*}
V_4&= -{a_2}^{8}Q{(\mu)}/
36887588669248862905116889599301849674634285420053378493\\
&\quad 1286682303008351322973629074374337069105705808904284987022190533816\\
&\quad 3182367071153612906739577075384387480378163200000,
\end{align*}
where
\begin{align*}
Q(\mu)
&= 5061565304752756281507283445466696498397806982562902502862568707\\
&\quad 12912529775726964251246968509875875991228610701216698005689607424\\
&\quad 73087731876897344650076374604360705836087040{\mu}^{7}\\
&\quad -340812912483222809928339721244205680969631979595064645753790303\\
&\quad 68058287137415281782406391986081459447638560645736575391889270596\\
&\quad 577079487175148071589735417325410968915137744{\mu}^{6}\\
&\quad - 824956771772875920524768605310963288557474695808927379648635146\\
&\quad 89934278817997190177837875064854601374617126333206151559276636649\\
&\quad 98789597343778156617 9221118658577582087242584{\mu}^{5}\\
&\quad + 81584483608286871615772633434382546877277513771408741974445245\\
&\quad 4845368677016347807819376306567104009462092441916792312570408246\\
&\quad 02372031012019997815662870970633944875629927764\,{\mu}^4\\
&\quad + 25896513600610435149766065028213468848105899037674273723965973\\
&\quad 85292517796469005762508262417287664215744015726820682298861228047\\
&\quad 9466387807108601473511477377951770327152251930{\mu}^{3}\\
&\quad -63871971031084648777832890365823291014772224724716750214834580\\
&\quad 9852821152694191443403735777254252142582767175529895793775396585\\
&\quad 05625608488378226940819951947411084516792539928{\mu}^2\\
&\quad +63320363852535951596180397565233937233465160852321916775662978\\
&\quad 7014918692824778554364353862099348531370990541173801942371682611\\
&\quad 4200270000631493834542751507656379882694394067\mu\\
&\quad +15864826572799943856438557624408837065205007313716175374949797\\
&\quad 7723864079669654897429092964389267887437504915634854709787668372\\
&\quad 90103030523917020740911403903049474302071125600.
\end{align*}
Computing the resultant $R(P, Q)$ between two polynomials $P(\lambda),Q(\lambda)$, 
we obtain
$R(P,Q)\neq{0}$.
Thus $V_4\neq {0}$ and  the origin is a fourth order fine focus for system 
\eqref{sys27}.
\end{proof}

\begin{theorem} \label{thm5}
If  condition \eqref{cond} holds, then there are perturbations of system
\eqref{sys27} yielding  three  small-amplitude limit cycles bifurcating
from the origin.
\end{theorem}

\begin{proof}
Under  condition \eqref{cond}, the Jacobian matrix of
focal values $V_1,V_2,V_3$  of system \eqref{sys27} with respect to
$a_1,a_3,a_4$ has full rank, i.e.,
\[
\operatorname{rank} \big[\frac{\partial
{(V_1,V_2,V_3})}{\partial
(a_1,a_3,a_4)}\big]_{\eqref{cond}}=3,
\]
hence by \cite[Theorem 2.3.2]{Han13} the claim follows.
\end{proof}


\subsection*{Acknowledgements}
 This work was supported by a grant from the National Natural Science 
Foundation of China (No. 11461021), a grant from the Foundation 
for Research in Experimental Techniques of Liaocheng University (No. LDSY2014110),
and a grant from the Scientific Research Foundation of Guangxi Education 
Department (No. ZD2014131). 
The author is grateful to the referee for the valuable remarks which helped 
to improve the manuscript.

\begin{thebibliography}{99}

\bibitem{bvl} L. Barreira,  C. Valls,  J. Llibre;
\emph{Integrability and limit cycles of the Moon-Rand
system}, Int. J. Nonlin. Mech., \textbf{69}(2015), 129--136.

\bibitem{gy09} M. Gyllenberg, P. Yan;
 \emph{Four limit cycles for a three-dimensional competitive Lotka-Volterra
system with a heteroclinic cycle}, Comput. Math. Appl., \textbf{
{58}}(2009), 649--669.

\bibitem{Han13} M. A. Han;
\emph{Bifurcation Theory of Limit Cycles}, Mathematics Monograph Series, 
vol. 25, Science Press, Beijing, 2013.

\bibitem{han15}M. A. Han, P. Yu;
\emph{Ten limit cycles around a center-type singular point in a 3-d
quadratic system with quadratic perturbation}, Appl. Math. Lett.,
{\textbf{44}}(2015), 17--20.

\bibitem{k04} Y. A. Kuznetsov;
\emph{Elements of Applied Bifurcation Theory}, Springer-Verlag, New York, 2004.

\bibitem{l92} J. K. Liu, L. C. Zhao;
\emph{Bifurcation analysis of airfoils in incompressible flow},  
Journal of Sound and Vibration, \textbf{{154}}(1992), 117--124.

\bibitem{l01} Y. R. Liu;
\emph{Theory of center-focus for a class of higher-degree critical points and 
infinite  points}, Sci. China Ser. A, \textbf{{44}}(2001), 37--48.

\bibitem{l08} Y. R. Liu, J. B. Li, W. T. Huang;
\emph{Singular point values, center problem and bifurcations of
limit cycles of two dimensional differential autonomous systems},
Nonlinear Dynamics Series, Vol. 6, Science Press, Beijing, 2008.

\bibitem{lv} J. Llibre,  C. Valls;
\emph{Hopf bifurcation of a generalized Moon-Rand system}, 
Commun. Nonlinear Sci. Numer. Simulat., \textbf{20}(2015), 1070--1077.

\bibitem{mrs} A. Mahdi, V. G. Romanovski, D. S. Shafer;
 \emph{Stability and periodic oscillations in the Moon-Rand systems}, 
Nonlinear Anal. Real World Appl., \textbf{14}(2013), 294--313.


\bibitem{s14} B. Sang;
\emph{Center problem for a class of degenerate
quartic systems}, Electron. J. Qual. Theory  Differ. Equ., no.
74(2014), 1--17.

\bibitem{t14} Y. Tian, P. Yu;
\emph{Seven limit cycles around a focus point in a simple
three-dimensional quadratic vector field}, Int. J. Bifurcation and
Chaos, \textbf{{24}}(2014),  1450083--1450092.

\bibitem{yu98} P. Yu;
\emph{Computation of normal forms via a perturbation technique}, 
Journal of Sound and Vibration, \textbf{{211}}(1998), 19--38.

\bibitem{yu13} P. Yu, M. A. Han;
\emph{Eight limit cycles around a center
in quadratic hamiltonian system with third-order perturbation},
Int. J. Bifurcation and Chaos, \textbf{{23}}(2013),
1350005--1350022.

\bibitem{w10} Q. L. Wang,  Y. R. Liu,  H. B. Chen;
\emph{Hopf bifurcation for a class of three-dimensional nonlinear
dynamic systems}, Bull. Sci. Math., \textbf{{134}}(2010), 786--798.

\bibitem{z04} Q. C. Zhang, H. Y. Liu, A. D. Ren;
\emph{The study of limit cycle flutter for airfoil with
nonlinearity}, Acta Aerodynamica Sinica, \textbf{22}(2004), 
332--336 (in Chinese).

\end{thebibliography}

\end{document}