Mutiple Hopf Bifurcation on Center Manifold

Abstract

In this chapter, by researching the algorithm of the formal series, and deducing the recursion formula of computing the nondegenerate and degenerate singular point quantities on center manifold, we investigate the Hopf bifurcation of high‐dimensional nonlinear dynamic systems. And more as applications, the singular point quantities for two classes of typical three‐ or four‐dimensional polynomial systems are obtained, the corresponding multiple limit cycles or Hopf cyclicity restricted to the center manifold are discussed.

1. Introduction

This chapter is concerned with Hopf bifurcation restricted to the center manifold from the equilibrium for three‐, four‐, and more higher‐dimensional nonlinear dynamical systems.

Let us first consider the generic real systems which take the form

$\begin{array}{l}\dot{\mathbf{x}}=Ax+f(x)\hfill \end{array}$E1

where $x=(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_n)\in {\mathbb{R}}^n$, $A\in {\mathbb{R}}^{n\times n}\mathrm{,}n\in \mathbb{N}$, and $f(x)$ is sufficiently smooth with $f(0)=0$, $Df(0)=0$. Then the origin is an equilibrium. For dynamical analysis of systems (1), it is very important to discuss the asymptotic behavior and the existence of periodic orbits at the origin. When the Jacobi matrix $A$ has an eigenvalue with zero real part, the phase portraits in the vicinity of the origin is not easy to be determined. In particular, a system (1) has the following form

$\{\begin{array}{l}{\dot{\mathbf{x}}}_1=A_1{x}_1+f_1(x_1,x_2)\hfill \\ {\dot{\mathbf{x}}}_2=A_2{x}_2+f_2(x_1,x_2)\hfill \end{array}$E2

where $x_1=(x_1\mathrm{,}{x}_2\mathrm{,}\dots \mathrm{,}{x}_{n_c}{)}^T\in {\mathbb{R}}^{n_c}\mathrm{,}{x}_2=(x_{n_c+1}\mathrm{,}\dots \mathrm{,}{x}_n{)}^T\in {\mathbb{R}}^{n_s}$ with $n_c+n_s=n,$ $A_1$ and $A_2$ are constant matrices, and $f_1(x_1,x_2)$, $f_2(x_1,x_2)$ are functions with

$f_1(0,0)=0\mathrm{,}{f}_2(0,0)=0\mathrm{,}D{f}_1(0,0)=0\mathrm{,}D{f}_2(0,0)=0$E42

Suppose that $A_1$ has $n_c$ critical eigenvalues (i.e., eigenvalues with Re $\lambda$ = 0) and all $n_s$ eigenvalues of $A_2$ satisfy Re $\lambda <0$. According to the Center Manifold Theorem (see, e.g., [1, 2]), there exists a (local) center manifold $x_2=\mathbf{h}({\mathbf{x}}_1)$ with $h(0)=0\mathrm{,}Dh(0)=0,$ and system (2) is topologically equivalent near $(\mathbf{0},\mathbf{0})$ to the system

$\{\begin{array}{l}{\dot{\mathbf{x}}}_1=A_1{x}_1+f_1(x_1,\mathbf{h}(x_1))\hfill \\ {\dot{\mathbf{x}}}_2=A_2{x}_2\mathrm{.}\hfill \end{array}$E3

The first equation in Eq. (3) is called the restriction of system (2) to its center manifold at the origin. The local center manifold, which is tangent to the $(x_1\mathrm{,}{x}_2\mathrm{,}\dots \mathrm{,}{x}_{n_c})$‐plane (hyperplane) at the origin and which contains all the recurrent behavior of system (2) in a neighborhood of the origin, since the second equation in (3) is linear and has exponentially decaying solutions (see, e.g., [3]). Thus, the dynamics of Eq. (2) near a nonhyperbolic equilibrium are determined by this restriction. Generally, the local center manifold is not necessarily unique, but if the origin is a center restricted to a local center manifold for system (2), then the center manifold is unique and analytic, which is presented by the Lyapunov Center Theorem proved in Ref. [4].

If $A$ has a simple pair of purely imaginary eigenvalues $\pm \omega \mathrm{i}$ ($\omega >0$), system (1) undergoes a Hopf bifurcation or multiple Hopf bifurcation in a neighborhood of the origin on the local center manifold under proper perturbations of parameters. The computation of focal values (Lyapunov coefficients) plays an important role in the study of small‐amplitude limit cycles appearing in these bifurcations (see [5–14] and references therein). The projection method was used for computing the first and the second focal values (see [2]), and a perturbation technique based on multiple time scales was used for computing focal values (see [15]). For a class of three‐dimensional systems, the formal series method was presented with a recursive formula for computing singular point quantities (see [16]), here the theory and methodology described in Refs. [16, 17] can be applied to $n$‐dimensional systems, where $n\ge 4$.

If $A$ has some zero eigenvalues for system (1), the Hopf bifurcation problem at the origin on the local center manifold becomes generally more difficult in comparison to the nondegenerate case. Take the degenerate singular point with a zero linear part in planar system, for example, the investigation of Hopf bifurcation from the equilibrium has to involve detecting the monodromy and distinguishing between a center and a focus [18, 19]. For that matter, several available approaches and corresponding results can be seen in [18–25], and one can easily find that the results on the bifurcation of limit cycles are very less. Remarkably, the author of reference [26] in 2001 gave the formal series method of calculating the singular point quantities of the degenerate critical point, which made it possible to investigate multiple Hopf bifurcation of higher degree polynomial systems [27, 28]. Here we extend its application to the local center manifold of more higher‐dimensional system.

2. Case of the nondegenerate singular point

In this section, we consider Hopf bifurcation from the nondegenerate origin of system (1) restricted to the center manifold, in which the Jacobian matrix $A$ has a pair of pure imaginary eigenvalues and its other eigenvalues are all negative. As the particular case, for planar systems there exist some good computer algebra procedure to calculate the focal values (see survey article [29], monograph [30], and references therein), here the formal series method of computing singular point quantities on the local center manifold for high‐dimensional system originated from the work of [31–33] in planar systems.

2.1. The formal series method of computing nondegenerate singular point quantities on center manifold

Considering the Jacobian matrix $A$ at the origin of system (1) has a pair of purely imaginary eigenvalues and a negative one, then by certain nondegenerate transformation, the system (1) can be changed into the following system:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=-y+{\displaystyle \sum _{k+j+l=2}^{\infty }}{A}_{kjl}{x}^k{y}^j{u}^l=X(x\mathrm{,}y\mathrm{,}u),\\ \frac{\text{d}y}{\text{d}t}=x+{\displaystyle \sum _{k+j+l=2}^{\infty }}{B}_{kjl}{x}^k{y}^j{u}^l=Y(x\mathrm{,}y\mathrm{,}u),\\ \frac{\text{d}u}{\text{d}t}=-d_{0}u+{\displaystyle \sum _{k+j+l=2}^{\infty }}{\tilde{d}}_{kjl}{x}^k{y}^j{u}^l=\tilde{U}(x\mathrm{,}y\mathrm{,}u)\end{array}$E4

where $x\mathrm{,}y\mathrm{,}u\mathrm{,}{A}_{kjl}\mathrm{,}{B}_{kjl}\mathrm{,}{\tilde{d}}_{kjl}\in \mathbb{R}(k\mathrm{,}j\mathrm{,}l\in \mathbb{N})$ and $d_0>0$.

Here, we recall first the calculation method of the singular point quantities on center manifold for the above real three‐dimensional nonlinear dynamical systems. By means of transformation

$z=x+yi,w=x-yi,u=u,T=it,i=\sqrt{-1}$E5

system (4) is also transformed into the following complex system:

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}T}=z+{\displaystyle \sum _{k+j+l=2}^{\infty }}{a}_{kjl}{z}^k{w}^j{u}^l=Z(z\mathrm{,}w\mathrm{,}u),\\ \frac{\text{d}w}{\text{d}T}=-w-{\displaystyle \sum _{k+j+l=2}^{\infty }}{b}_{kjl}{w}^k{z}^j{u}^l=-W(z\mathrm{,}w\mathrm{,}u),\\ \frac{\text{d}u}{\text{d}T}=\text{i}{d}_{0}u+{\displaystyle \sum _{k+j+l=2}^{\infty }}{d}_{kjl}{z}^k{w}^j{u}^l=U(z\mathrm{,}w\mathrm{,}u)\end{array}$E6

where $z\mathrm{,}w\mathrm{,}T\mathrm{,}{a}_{kjl}\mathrm{,}{b}_{kjl}\mathrm{,}{d}_{kjl}\in \mathbb{C}(k\mathrm{,}j\mathrm{,}l\in \mathbb{N})$, the systems (4) and (6) are called concomitant.

Theorem 1 (see [16]). For system (6), using the program of term by term calculations, we can determine a formal power series:

$\begin{array}{l}F(z\mathrm{,}w\mathrm{,}u)=zw+{\displaystyle \sum _{\alpha +\beta +\gamma =3}^{\infty }}{c}_{\alpha \beta \gamma }{z}^{\alpha }{w}^{\beta }{u}^{\gamma }\hfill \end{array}$E7

such that

$\frac{\text{d}F}{\text{d}T}=\frac{\partial F}{\partial z}Z-\frac{\partial F}{\partial y}W+\frac{\partial F}{\partial u}U={\displaystyle \sum _{m=1}^{\infty }}{\mu }_m{(zw)}^{m+1}$E8

where $c_{110}=\mathrm{1,}{c}_{101}=c_{011}=c_{200}=c_{020}=\mathrm{0,}{c}_{kk0}=0$, $k=\mathrm{2,\; 3,}\cdots$.

Definition 1. The ${\mu }_m$ in the expression (8) is called the $m$th singular point quantity at the origin on center manifold of system (6) or (4), $m=\mathrm{1,2,}\cdots$.

Theorem 2 (see [16, 34]). For the $m$th singular point quantity and the $m$th focal value at the origin on center manifold of system (4), i.e., ${\mu }_m$ and $v_{2m+1}\mathrm{,}m=\mathrm{1,2,}\cdots$, we have the following relation:

$v_{2m+1}(2\pi )=i\pi {\mu }_m+i\pi {\displaystyle \sum _{k=1}^{m-1}}{\xi }_m^{(k)}{\mu }_k$E9

where ${\xi }_m^{(k)}(k=\mathrm{1,2,}\cdots \mathrm{,}m-1)$ are polynomial functions of coefficients of system (6). Usually, it is called algebraic equivalence and written as $v_{2m+1}\sim i\pi {\mu }_m$.

Based on the previous work in Ref. [16], we have developed the calculation method of the focal values on the center manifold for real four‐dimensional nonlinear dynamical systems in Ref. [35]. In fact, here Theorem 1 can be generalized in the $n$‐dimensional real systems as follows

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=-y+\text{h}\text{.o}\text{.t}\mathrm{.}=X(x\mathrm{,}y\mathrm{,}u),\hfill \\ \frac{\text{d}y}{\text{d}t}=x+\text{h}\text{.o}\text{.t}\mathrm{.}=Y(x\mathrm{,}y\mathrm{,}u),\hfill \\ \frac{\text{d}{u}_i}{\text{d}t}=-d_i{u}_i+\text{h}\mathrm{.}\text{o}\text{.t}\mathrm{.}={\tilde{U}}_i(x\mathrm{,}y\mathrm{,}u),i=\mathrm{1,2,}\cdots \mathrm{,}n-2\hfill \end{array}$E10

where $u=(u_1\mathrm{,}{u}_2\mathrm{,}\cdots \mathrm{,}{u}_{n-2})$, h.o.t denotes the terms in $x\mathrm{,}y\mathrm{,}{u}_1\mathrm{,}{u}_2\mathrm{,}\cdots \mathrm{,}{u}_{n-2}$ with orders greater than or equal to 2, and all $d_i>0$.

By means of transformation of Eq. (5), system (10) can be transformed into the following complex system

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}T}=z+{\displaystyle {\sum }_{{}_{k+j+1=2}}^{\infty }}{a}_{{}_{kj1}}{z}^k{w}^j{u}^1=Z(z\mathrm{,}w\mathrm{,}u),\hfill \\ \frac{\text{d}w}{\text{d}T}=-w-{\displaystyle {\sum }_{{}_{k+j+l=2}}^{\infty }}{b}_{{}_{kj1}}{w}^k{z}^j{u}^1=-W(z\mathrm{,}w\mathrm{,}u),\hfill \\ \frac{\text{d}{u}_i}{\text{d}T}=i{d}_i{u}_i+{\displaystyle {\sum }_{{}_{k+j+1=2}}^{\infty }}{d}_{{}_{kj1}}{z}^k{w}^j{u}^1=U_i(z\mathrm{,}w\mathrm{,}u),i=\mathrm{1,2,}\cdots \mathrm{,}n-2\hfill \end{array}$E11

where the subscript “$kj1$” denotes “$kj{l}_1\cdots l_{n-2}$”, $u^1=u_1^{l_1}{u}_2^{l_2}\cdots u_{n-2}^{l_{n-2}}$, and $l={\displaystyle \sum _{i=1}^{n-2}}{l}_i$, all $u_i\in \mathbb{R}$, $z\mathrm{,}w\mathrm{,}T\mathrm{,}{a}_{kj1}\mathrm{,}{b}_{kj1}\mathrm{,}{d}_{kj1}\in \mathbb{C}(k\mathrm{,}j\mathrm{,}{l}_i\in \mathbb{N})$, we call that system (10) and system (11) are concomitant.

Theorem 3. For system (11), using the program of term by term calculations, we can determine a formal power series:

$\begin{array}{l}F(z\mathrm{,}w\mathrm{,}\mathit{u})=zw+{\displaystyle \sum _{\alpha +\beta +\ell =3}^{\infty }}{c}_{\alpha \beta \ell }{z}^{\alpha }{w}^{\beta }{\mathit{u}}^{\ell }\hfill \end{array}$E12

such that

$\frac{\text{d}F}{\text{d}T}=\frac{\partial F}{\partial z}Z-\frac{\partial F}{\partial y}W+{\displaystyle \sum _{i=1}^{n-2}}\frac{\partial F}{\partial u_i}{U}_i={\displaystyle \sum _{m=1}^{\infty }}{\mu }_m{(zw)}^{m+1}$E13

where the subscript “$\alpha \beta \ell$” denotes “$\alpha \beta {\gamma }_1\cdots {\gamma }_{n-2}$”, ${\mathit{u}}^{\ell }=u_1^{{\gamma }_1}{u}_2^{{\gamma }_2}\cdots u_{n-2}^{{\gamma }_{n-2}}$, and $\ell ={\displaystyle \sum _{i=1}^{n-2}}{\gamma }_i$, and more setting $c_{\alpha \beta \ell }=0$ with $0\le \alpha +\beta +\ell \le 2$ except for $c_{110}=1$, and $c_{kk0}=0$ with $k\ge 2$.

Proof. It is very similar to the proving course of Theorem 1.3.1 in [16], by computing carefully and comparing the above power series with the two sides of (13), we can obtain the expression of ${\mu }_m$.

Definition 2. The ${\mu }_m$ in the expression (13) is called the $m$th singular point quantity at the origin on center manifold of system (11) or (10), $m=\mathrm{1,2,}\cdots$.

Remark 1. Similar to Theorem 2, there exists a equivalence between ${\mu }_m$ and $v_{2m+1}$, namely, if ${\mu }_1={\mu }_2=\cdots ={\mu }_{m-1}=\mathrm{0,}{\mu }_m\ne 0$, then $v_3=v_5=\cdots =v_{2m-1}=\mathrm{0,}{v}_{2m+1}=i\pi {\mu }_m\mathrm{,}m=\mathrm{1,2,}\cdots$, and vice versa.

Corollary 1. The origin of system (10) or (11) is a center restricted to the center manifold if and only if ${\mu }_m=0$ for all $m$.

Remark 2. From the relation given by Remark 1 and Corollary 1, the center‐focus problem and Hopf bifurcation of equilibrium point restricted to the center manifold can be figured out by the calculation of singular point quantities for system (10).

2.2. An example of four‐dimensional system

Recently, the study of chaos has become a hot research topic, and the attention of many researchers is turning to 4D systems from 3D dynamical systems, for example, the authors of Ref. [36] investigated Hopf bifurcation of a 4D‐hyoerchaotic system by applying the normal form theory in 2012, but its multiple Hopf bifurcation on the center manifold have not been considered. Here, we will investigate the system further by computing the singular point quantities of its equilibrium point, which takes the following form

$\{\begin{array}{l}{\dot{x}}_1=a(x_2-x_1)\\ {\dot{x}}_2=c{x}_1-x_2+x_4-x_1{x}_3\\ {\dot{x}}_3=x_1{x}_2-b{x}_3+e{x}_1^2\\ {\dot{x}}_4=-K{x}_2\end{array}$E14

where $a\mathrm{,}b\mathrm{,}c\mathrm{,}e\mathrm{,}K\in \mathbb{R}$. Obviously, system (14) has only one isolated equilibrium: $O(\mathrm{0,0,0,0})$ when $K\ne 0$. Therefore, we only need to consider $O$. The Jacobian matrix of system (14) at $O$ is

$A=\left(\begin{array}{cccc}-a& a& 0& 0\\ c& -1& 0& 1\\ 0& 0& -b& 0\\ 0& -K& 0& 0\end{array}\right)$E43

with the characteristic equation:

$(\lambda +b)({\lambda }^3+(a+1){\lambda }^2+(a-ac+K)\lambda +aK=0.$E15

To guarantee that $A$ has a pair of purely imaginary eigenvalues $\pm i\omega (\omega >0)$ and two negative real eigenvalues ${\lambda }_1\mathrm{,}{\lambda }_2$, we let its characteristic equation take the form

$({\lambda }^2+{\omega }^2)(\lambda -{\lambda }_1)(\lambda -{\lambda }_2)=0.$E44

Thus, we obtain the critical condition of Hopf bifurcation at $O$:

$\begin{array}{l}{a}^2(c-1)={\omega }^2\mathrm{,}K=a(a+1)(c-1),{\lambda }_1=-b\mathrm{,}{\lambda }_2=-a-1\hfill \end{array}$E16

where $a>-\mathrm{1,}b>\mathrm{0,}c>1$, namely, $c=\frac{a^2+{\omega }^2}{a^2}\mathrm{,}K=\frac{(a+1){\omega }^2}{a}$. Under the conditions (16), one can find a nondegenerate matrix

$P=\left(\begin{array}{cccc}-\frac{i{a}^2}{(a+1)(a+i\omega )\omega }& \frac{i{a}^2}{(a+1)(a-i\omega )\omega }& 0& -\frac{a^2}{{\omega }^2}\\ -\frac{ia}{(a+1)\omega }& \frac{ia}{a\omega +\omega }& 0& \frac{a}{{\omega }^2}\\ 0& 0& 1& 0\\ 1& 1& 0& 1\end{array}\right)$E45

such that

$P^{-1}AP=\left(\begin{array}{cccc}\omega i& 0& 0& 0\\ 0& -\omega i& 0& 0\\ 0& 0& -b& 0\\ 0& 0& 0& -a-1\end{array}\right)$E17

Namely, we can use the nondegenerate transformation and the time rescaling: $T=it\omega$ to make the system (14) become the following same form as the complex system (11) with $n=4$:

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}T}=z+{\displaystyle {\sum }_{{}_{k+j+l+n=2}}^2}{a}_{{}_{kjln}}{z}^k{w}^j{u}^l{v}^n=Z(z\mathrm{,}w\mathrm{,}u\mathrm{,}v),\hfill \\ \frac{\text{d}w}{\text{d}T}=-w-{\displaystyle {\sum }_{{}_{k+j+l+n=2}}^2}{b}_{{}_{kjln}}{w}^k{z}^j{u}^l{v}^n=-W(z\mathrm{,}w\mathrm{,}u\mathrm{,}v),\hfill \\ \frac{\text{d}u}{\text{d}T}=\frac{bi}{\omega }u+{\displaystyle {\sum }_{{}_{k+j+l+n=2}}^2}{d}_{{}_{kjln}}{z}^k{w}^j{u}^l{v}^n=U(z\mathrm{,}w\mathrm{,}u\mathrm{,}v),\hfill \\ \frac{\text{d}v}{\text{d}T}=\frac{(a+1)i}{\omega }v+{\displaystyle {\sum }_{{}_{k+j+l+n=2}}^2}{e}_{{}_{kjln}}{z}^k{w}^j{u}^l{v}^n=V(z\mathrm{,}w\mathrm{,}u\mathrm{,}v)\hfill \end{array}$E18

where $u\in \mathbb{R}$, $z\mathrm{,}w\mathrm{,}T\in \mathbb{C}$, and all $a_{kjln}=b_{kjln}=d_{kjln}=e_{kjln}=0$ except the following coefficients

$\begin{array}{l}{a}_{0011}=\frac{a^3+a^2(1+i\omega )+ia\omega }{2{\omega }^2(a+i\omega +1)}\mathrm{,}{a}_{0110}=\frac{a(\omega -ia)}{2\omega \left(a^2+a+\omega (\omega -i)\right)}\mathrm{,}\hfill \\ b_{kjln}={\overline{a}}_{kjln}(ikjl=\mathrm{0011,0110}),\hfill \\ d_{0002}=\frac{i{a}^3(1-a)e}{{\omega }^5}\mathrm{,}{d}_{0101}=-\frac{a^4(2e+1)-a^3(1+i\omega )}{(a+1){\omega }^4(a-i\omega )}\mathrm{,}\hfill \\ d_{0200}=\frac{a^3\omega +i{a}^4(e+1)}{{(a+1)}^2{\omega }^3{(a-i\omega )}^2}\mathrm{,}{d}_{1001}=\frac{a^4(2e+1)+a^3(i\omega -1)}{(a+1){\omega }^4(a+i\omega )}\mathrm{,}\hfill \\ d_{1100}=-\frac{2i{a}^4(e+1)}{{(a+1)}^2{\omega }^3\left(a^2+{\omega }^2\right)}\mathrm{,}{d}_{2000}=-\frac{a^3\omega -i{a}^4(e+1)}{{(a+1)}^2{\omega }^3{(a+i\omega )}^2}\mathrm{,}\hfill \\ e_{0011}=-\frac{ia(a+1)}{\omega \left(a^2+2a+{\omega }^2+1\right)}\mathrm{,}{e}_{0110}=-\frac{a}{(a-i\omega )\left(a^2+2a+{\omega }^2+1\right)}\mathrm{,}\hfill \\ e_{1010}=\frac{a}{(a+i\omega )\left(a^2+2a+{\omega }^2+1\right)}\hfill \end{array}$E46

where ${\overline{a}}_{kjln}$ denotes the conjugate complex number of $a_{kjln}$.

According to Theorem 3, we obtain the recursive formulas of $c_{\alpha \beta \gamma }$ and ${\mu }_m$.

Theorem 5. For system (18), setting $c_{\alpha \beta \gamma \lambda }=0$ with $0\le \alpha +\beta +\gamma +\lambda \le 2$ except for $c_{1100}=1$, and $c_{kk00}=0$ with $k\ge 2$, we can derive successively and uniquely the terms of the following formal series (12) with $n=4$, such that (13) with $n=4$ holds and if $\alpha \ne \beta$ or $\alpha =\beta \mathrm{,}{\lambda }^2+{\gamma }^2\ne 0$, $c_{\alpha \beta \gamma \lambda }$ is determined by following recursive formula:

$\begin{array}{l}{c}_{{}_{\alpha \beta \gamma \lambda }}=\frac{\omega }{\omega (\alpha -\beta )+i(b\gamma +(a+1)\lambda )}\hfill \\ \{-d_{2000}(1+\gamma )c[\alpha -\mathrm{2,}\beta \mathrm{,}\gamma +\mathrm{1,}\lambda ]-d_{1100}(\gamma +1)c[\alpha -\mathrm{1,}\beta -\mathrm{1,}\gamma +\mathrm{1,}\lambda ]-\hfill \\ e_{1010}(\lambda +1)c[\alpha -\mathrm{1,}\beta \mathrm{,}\gamma -\mathrm{1,}\lambda +1]-d_{1001}(\gamma +1)c[\alpha -\mathrm{1,}\beta \mathrm{,}\gamma +\mathrm{1,}\lambda -1]+\hfill \\ b_{0110}(\beta +1)c[\alpha -\mathrm{1,}\beta +\mathrm{1,}\gamma -\mathrm{1,}\lambda ]-d_{0200}(\gamma +1)c[\alpha \mathrm{,}\beta -\mathrm{2,}\gamma +\mathrm{1,}\lambda ]-\hfill \\ e_{0110}(\lambda +1)c[\alpha \mathrm{,}\beta -\mathrm{1,}\gamma -\mathrm{1,}\lambda +1]-d_{0101}(\gamma +1)c[\alpha \mathrm{,}\beta -\mathrm{1,}\gamma +\mathrm{1,}\lambda -1]-\hfill \\ e_{0011}\lambda c[\alpha \mathrm{,}\beta \mathrm{,}\gamma -\mathrm{1,}\lambda ]-d_{0002}(\gamma +1)c[\alpha \mathrm{,}\beta \mathrm{,}\gamma +\mathrm{1,}\lambda -2]+\hfill \\ b_{0011}(\beta +1)c[\alpha \mathrm{,}\beta +\mathrm{1,}\gamma -\mathrm{1,}\lambda -1]-a_{0110}(\alpha +1)c[\alpha +\mathrm{1,}\beta -\mathrm{1,}\gamma -\mathrm{1,}\lambda ]-\hfill \\ a_{0011}(\alpha +1)c[\alpha +\mathrm{1,}\beta \mathrm{,}\gamma -\mathrm{1,}\lambda -1]\}\hfill \end{array}$E19

and for any positive integer $m\mathrm{,}{\mu }_m$ is determined by following recursive formula:

$\begin{array}{l}{\mu }_m=d_{2000}c[-2+m\mathrm{,}m\mathrm{,1,0}]\hfill \\ +d_{1100}c[-1+m\mathrm{,}-1+m\mathrm{,1,0}]+d_{0200}c[m\mathrm{,}-2+m\mathrm{,1,0}]\hfill \end{array}$E20

and when $\alpha <0$ or $\beta <0$ or $\gamma <0$ or $\lambda <0$ or $\alpha =\beta \mathrm{,}\gamma =\lambda =0$, we have let $c_{\alpha \beta \gamma \lambda }=0$, and where each $c[\alpha \mathrm{,}\beta \mathrm{,}\gamma \mathrm{,}\lambda ]$ denotes $c_{{}_{\alpha \beta \gamma \lambda }}$.

By applying the above formulas in the Mathematica symbolic computation system, we figure out easily the first two singular point quantities of the origin of system (18):

$\begin{array}{l}{\mu }_1=ia{f}_1{[|a|bc{(a+1)}^2{d}_0]}^{-1}\mathrm{,}\hfill \\ {\mu }_2=108i{a}^3{b}^4{f}_2{f}_3^2{f}_4{[|a|c^2{d}_0{d}_1{}^2{d}_2{}^4{d}_3]}^{-1}\hfill \end{array}$E21

where

$\begin{array}{l}{f}_1=8a^{3}ce+8a^{3}c-8a^{3}e-8a^3-2a^{2}bce+2a^{2}be+8a^{2}ce+8a^{2}c\hfill \\ -8a^{2}e-8a^2+a{b}^{2}c+3a{b}^{2}e+2a{b}^2+2abc-2ab+3b^{2}e+3b^2\mathrm{,}\hfill \\ f_2=(2a+b+{2)}^3(2ae+2a-b)(e+1),\hfill \\ f_3=4a^{2}e+4a^2-3abe-2ab+4ae+4a+b\mathrm{,}\hfill \\ f_4=8a^5{c}^2-16a^{5}c+8a^5-2a^{4}b{c}^2+2a^{4}bc+8a^4{c}^2-16a^{4}c+8a^4+2a^3{b}^{2}c\hfill \\ -2a^3{b}^2-4a^{3}bc+4a^{3}b-5a^2{b}^{3}c+4a^2{b}^3+2a^2{b}^{2}c\hfill \\ \begin{array}{l}-2a^2{b}^2-2a^{2}bc+2a^{2}b-2a{b}^3-b^3\mathrm{,}\\ \begin{array}{l}{d}_0=(a^{2}c+2a+1)(4a^{2}c-4a^2+b^2)(c-{1)}^{3/2}\mathrm{,}\hfill \\ d_1=8a^{3}c-8a^3-2a^{2}bc+2a^{2}b+8a^{2}c-8a^2+3a{b}^2+3b^2\mathrm{,}\hfill \\ d_2=8a^{2}e+8a^2-2abe+8ae+8a+b^2+2b\mathrm{,}\hfill \\ d_3=9a^{2}c-8a^2+2a+\mathrm{1,}\hfill \end{array}\end{array}\hfill \end{array}$E47

and the above expression of ${\mu }_2$ is obtained under the condition of ${\mu }_1=0$.

From Remark 1 and the singular point quantities (21), we have

Theorem 6. For the flow on center manifold of the system (14), the first 2 focal values of the origin are as follow

$v_3=\mathit{i}\mathit{\pi }{\mathit{\mu }}_1\mathrm{,}{v}_5=\mathit{i}\mathit{\pi }{\mathit{\mu }}_2$E22

where the expression of $v_5$ is obtained under the condition of $v_3=0$.

Remark 3. In contrast to the result and process in [36], one can easily see that our first quantity is basically consistent with its characteristic exponent of bifurcating periodic solutions, and our algorithm is easy to realize with computer algebra system due to the linear recursion formulas, and more convenient to investigate the multiple Hopf bifurcation on center manifold.

Considering its Hopf bifurcation form of Theorem 6, we have the following:

Theorem 7. At least two small limit cycles can be bifurcated from the origin of the 4D‐hyoerchaotic system (14), which lie in the neighborhood of the origin restricted to the center manifold.

The rigorous proof of the above theorem is very similar to the previous ones in [14, 16], namely, by calculating the Jacobian determinant with respect to the functions $v_3\mathrm{,}{v}_5$ and its variables, which will not be given here.

3. Case of the degenerate singular point

Up till now, study on bifurcation of limit cycles from the degenerate singularity of higher dimensional nonlinear systems (1) is hardly seen in published references. Here, we will investigate the Hopf bifurcation problem from the high‐order critical point on the center manifold.

3.1. The formal series method of computing degenerate singular point quantities on center manifold

Let us consider the real $n$‐dimensional systems with two zero eigenvalues and zero linear part as follows

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=(\delta x-y)(x^2+y^2{)}^q+{\displaystyle \sum _{k+j+1=2q+2}^{\infty }}{A}_{kj1}{x}^k{y}^j{u}^1=X(x\mathrm{,}y\mathrm{,}u),\hfill \\ \frac{\text{d}y}{\text{d}t}=(x-\delta y)(x^2+y^2{)}^q+{\displaystyle \sum _{k+j+1=2q+2}^{\infty }}{B}_{kj1}{x}^k{y}^j{u}^1=Y(x\mathrm{,}y\mathrm{,}u),\hfill \\ \frac{\text{d}{u}_i}{\text{d}t}=-d_i{u}_i+{\displaystyle \sum _{k+j+1=2}^{\infty }}{d}_{kjl}{z}^k{w}^j{u}^1=U_i(x\mathrm{,}y\mathrm{,}u),i=\mathrm{1,2,}\cdots \mathrm{,}n-2\hfill \end{array}$E23

where the subscript “$kj1$” denotes “$kj{l}_1\cdots l_{n-2}$”, $u^1=u_1^{l_1}{u}_2^{l_2}\cdots u_{n-2}^{l_{n-2}}$, and $l={\displaystyle \sum _{i=1}^{n-2}}{l}_i$, all $d_i>0$, $x\mathrm{,}y\mathrm{,}{u}_i\mathrm{,}t\mathrm{,}\delta \mathrm{,}{A}_{kjl}\mathrm{,}{B}_{kjl}\mathrm{,}{d}_{kjl}\in \mathbb{R}$, $q\mathrm{,}k\mathrm{,}j\mathrm{,}{l}_i\in \mathbb{N}$. Obviously, the origin of system (23) is a high‐order degenerate singular point with two zero eigenvalues and $n-2$ negative ones.

In order to discuss the calculation method of the focal values on center manifold of the system (23), from the center manifold theorem [1], we take an approximation to the center manifold:

$u=u(x\mathrm{,}y)=u_2(x\mathrm{,}y)+h\mathrm{.}o\mathrm{.}t\mathrm{.}$E24

where $u=(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_{n-2}{)}^T$, $u_2$ is a quadratic homogeneous polynomial vector in $x$ and $y$, and h.o.t. denotes the terms with orders greater than or equal to 3. Substituting $u=u(x\mathrm{,}y)$ into the equations of system (23), we obtain a real planar polynomial differential system as follows

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=(\delta x-y)(x^2+y^2{)}^q+{\displaystyle \sum _{k=2q+2}^{\infty }}{X}_k(x\mathrm{,}y)=\tilde{X}(x\mathrm{,}y),\hfill \\ \frac{\text{d}y}{\text{d}t}=(x-\delta y)(x^2+y^2{)}^q+{\displaystyle \sum _{k=2q+2}^{\infty }}{Y}_k(x\mathrm{,}y)=\tilde{Y}(x\mathrm{,}y)\hfill \end{array}$E25

where $X_k(x\mathrm{,}y),Y_k(x\mathrm{,}y)$ are homogeneous polynomials of degree $k$, and the origin is degenerate with a zero linear part.

For system (25), some significant works have been done in Refs. [26] and [27]. Let us recall the related notions and results.

By means of transformation (5)

$\begin{array}{l}z=x+yi\mathrm{,}w=x-yi\mathrm{,}u=u\mathrm{,}T=it\mathrm{,}i=\sqrt{-1}\mathrm{,}\hfill \end{array}$

system (25) is transformed into following system:

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}T}=(1-i\delta )z^{q+1}{w}^q+{\displaystyle \sum _{k+j=2q+2}^{\infty }}{a}_{kj}{z}^k{w}^j=Z(z\mathrm{,}w),\hfill \\ \frac{\text{d}w}{\text{d}T}=-(1+i\delta )z^q{w}^{q+1}-{\displaystyle \sum _{k+j=2q+2}^{\infty }}{b}_{kj}{z}^k{w}^j=-W(z\mathrm{,}w)\hfill \end{array}$E26

where $z\mathrm{,}w\mathrm{,}T$ are complex variables and for any positive integer $k\mathrm{,}j$, we have $a_{kj}={\overline{b}}_{kj}$, then systems (25) and (26) are called concomitant.

For any positive integer $k$, we denote

$f_k(z\mathrm{,}w)={\displaystyle \sum _{\alpha +\beta =k}}{c}_{\alpha \beta }{z}^{\alpha }{w}^{\beta }$E48

a homogeneous polynomial of degree $k$ with $c_{00}=\mathrm{1,}{c}_{kk}=\mathrm{0,}k=\mathrm{1,\; 2,}\cdots$.

Theorem 8 ([26, 27]). For system (26) with $\delta =0$, we can derive successively the terms of the following formal series:

$F(z\mathrm{,}w)=zw[1+{\displaystyle \sum _{m=1}^{\infty }}\frac{f_{m(2q+3)}(z\mathrm{,}w)}{{(zw)}^{m(q+1)}}]$E27

such that

$\frac{\text{d}F}{\text{d}T}=\frac{\partial F}{\partial z}Z-\frac{\partial F}{\partial w}W=(zw{)}^q{\displaystyle \sum _{m=1}^{\infty }}{\mu }_m{(zw)}^{m+1}\mathrm{.}$E28

Definition 3. If $\delta =0$ holds, ${\mu }_m$ in expression (28) is called the $m$th singular point quantity at the degenerate singular point for system (26) or (1.3.26) is also called the $m$th singular point quantity of the origin on the center manifold of system (23), where $m=\mathrm{1,2,}\cdots \mathrm{.}$

Similar to Theorem 2, there also exists a equivalence between the $m$th singular point quantity and the $m$th focal value $v_{2m+1}(2\pi )$ at the origin on center manifold of system (23).

Theorem 9. For system (23) with $\delta =0$, and any positive integer m, the following assertion holds: $v_{2m+1}(2\pi )\sim i\pi {\mu }_m$, namely

$v_{2m+1}(2\pi )=i\pi ({\mu }_m+{\displaystyle \sum _{k=1}^{m-1}}{\xi }_m^{(k)}{\mu }_k),$E29

where ${\xi }_m^{(k)}(k=\mathrm{1,2,}\cdots \mathrm{,}m-1)$ are polynomial functions of coefficients of system (26). Then, the relation between $v_{2m+1}(2\pi )$ and ${\mu }_m$ is called the algebraic equivalence.

Remark 4. In fact, from Theorem 2, for any positive integer $m=\mathrm{2,3,}\cdots$, if ${\mu }_1={\mu }_2=\cdots ={\mu }_{m-1}=0$ and $v_1(2\pi )=v_3(2\pi )\cdots =v_{2m-1}(2\pi )=0$ hold, and vice versa. And more the stability and bifurcation of the origin of system (23) can be figured out by calculating the singular point quantities.

Corollary 2. The origin of system (23) is a center restricted to the center manifold if and only if ${\mu }_m=0$ for all $m$.

3.2. An example of three‐dimensional system

Now we consider an example for system (23) with $n=3$, it can be put in its concomitant form as follows

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}T}=(1-i\delta )z^{2}w+uz(a_{20}{z}^2+a_{11}zw+a_{02}{w}^2)=Z\mathrm{,}\\ \frac{\text{d}w}{\text{d}T}=-(1+i\delta )z{w}^2-uw(b_{20}{w}^2+b_{11}wz+b_{02}{z}^2)=-W\mathrm{,}\\ \frac{\text{d}u}{\text{d}T}=iu+i{d}_{1}zw=U\mathrm{,}\end{array}$E30

where $d_1\ne 0$ and

$a_{ij}=A_i+i{B}_i\mathrm{,}{b}_{ij}=A_i-i{B}_i\mathrm{,}{A}_i\mathrm{,}{B}_i\in \mathbb{R}\mathrm{,}i\mathrm{,}j=\mathrm{0,1,2,}$E31

namely, $a_{ij}={\overline{b}}_{ij}$. Then for the center manifold of system (30), from the transformation (5), we can determine the formal expression (24): $u=u(x\mathrm{,}y)=\tilde{u}(z\mathrm{,}w)$, thus obtain

$\{\begin{array}{c}\frac{\text{d}z}{\text{d}T}=(1-i\delta )z^{2}w+\tilde{u}z(a_{20}{z}^2+a_{11}zw+a_{02}{w}^2)=\tilde{Z},\hfill \\ \frac{\text{d}w}{\text{d}T}=-(1+i\delta )z{w}^2-\tilde{u}w(b_{20}{w}^2+b_{11}wz+b_{02}{z}^2)=-\tilde{W}.\hfill \end{array}$E32

Remark 5. For system (32), the corresponding $n=1$ in (27) and (28) of Theorem 8, we figure out that each ${\mu }_m$ is related to only the coefficients of the first $2m+3$ order terms of system (32), $m=\mathrm{1,2,}\cdots$. Here, we determine the above $\tilde{u}$ just to the sixth‐order term as follows

$\tilde{u}(z\mathrm{,}w)={\displaystyle \sum _{k=2}^6}{\tilde{u}}_k(z\mathrm{,}w)$E33

where ${\tilde{u}}_k$ is a homogeneous polynomial in $z\mathrm{,}w$ of degree $k$ and

$\begin{array}{l}{\tilde{u}}_2=-d_{1}zw\mathrm{,}{\tilde{u}}_4=2\delta d_1{z}^2{w}^2\mathrm{,}{\tilde{u}}_3={\tilde{u}}_4={\tilde{u}}_5=\mathrm{0,}\hfill \\ {\tilde{u}}_6=-i{d}_{1}wz((a_{02}-b_{20})d_1{w}^{3}z+(a_{11}{d}_1-b_{11}{d}_1-8i{\delta }^2)w^2{z}^2\hfill \\ +(a_{20}-b_{02})d_{1}w{z}^3).\hfill \end{array}$E34

Hence, $\tilde{Z}$ and $\tilde{W}$ in system (32) are two polynomials with degree 9.

Theorem 10. For system (32) with $\delta =0$, we can derive successively the terms of the formal series (27), such that (28) holds ($c_{\alpha \beta }$, ${\mu }_m$ in Appendix A).

Applying the powerful symbolic computation function of the Mathematica system and the recursive formulas in Theorem 10, and from Remark 5, we obtain the first three singular point quantities as follows

$\begin{array}{l}{\mu }_1=-d_1(a_{11}-b_{11}),\hfill \\ {\mu }_2=d_1^2(b_{20}{b}_{02}-a_{20}{a}_{02}),\hfill \\ {\mu }_3=-2i{d}_1^2(a_{02}{a}_{20}+b_{02}{b}_{20}-a_{02}{b}_{02}-a_{20}{b}_{20})\hfill \end{array}$E35

In the above expression of each ${\mu }_k\mathrm{,}k=\mathrm{2,3}$, we have already let ${\mu }_1=\cdots ={\mu }_{k-1}=0$.

Thus, from Theorem 9 and Eqs. (35) and (31), we have

Theorem 11. For the flow on center manifold of system (30),${}_{\delta =0}$, the first three focal values $v_{2i+1}(2\pi )(i=\mathrm{1,2,3})$ of the origin are as follows

$\begin{array}{l}{v}_3=2\pi d_1{B}_1\mathrm{,}\\ v_5=2\pi d_1^2(A_2{B}_0+A_0{B}_2),\\ v_7=2\pi d_1^2[(A_0-A_2{)}^2+(B_0+B_2{)}^2]\end{array}$E36

Theorem 12. For the flow on center manifold of (30),${}_{\delta =0}$, the origin is a three‐order weak focus, i.e., $v_3=v_5=\mathrm{0,}{v}_7\ne 0$ if and only if

$B_1=\mathrm{0,}{A}_2{B}_0+A_0{B}_2=0\text{and}{(A_0-A_2)}^2+{(B_0+B_2)}^2\ne 0$E37

Remark 6. For the coefficients of system (30)${}_{\delta =0}$, there exists necessarily a group of critical values: $A_i=A_i^{*}\mathrm{,}{B}_i=B_i^{*}(i=\mathrm{0,1,2})$ such that the conditions (37) hold, for example:

$A_1^{*}=B_1^{*}=\mathrm{0,}{A}_0^{*}=B_0^{*}=\mathrm{1,}{B}_2^{*}=-A_2^{*}=13$E38

Now we consider Hopf bifurcation of limit cycles from the origin for perturbed system (30).

Theorem 13. At least three limit cycles can be bifurcated from the origin of system (30) restricted to the center manifold, which lie in the neighborhood of the origin.

Proof. From Theorem 11, one can easily calculate the Jacobian determinant with respect to the functions $v_3\mathrm{,}{v}_5\mathrm{,}{v}_7$ and variables $B_1\mathrm{,}{B}_0\mathrm{,}{A}_0$,

$J=\frac{\partial (v_3\mathrm{,}{v}_5\mathrm{,}{v}_7)}{\partial (B_1\mathrm{,}{B}_0\mathrm{,}{A}_0)}=-2{\pi }^3{d}_1^5[8(A_0{A}_2-A_2^2-B_0{B}_2-B_2^2)]$E39

Considering the conditions (37) of Theorem 12 and substituting the group of critical values of Eq. (38) into Eq. (39), we obtain $J=649{\pi }^3{d}_1^5\ne 0$. Thus, we take some appropriate perturbations for the coefficients of system (32) to make the following two conditions:

$(v_1(2\pi )-1)v_3<\mathrm{0,}{v}_3{v}_5<\mathrm{0,}{v}_5{v}_7<0$E40

and

$\left|e^{2\pi \delta }-1\right|\ll \left|v_3\right|\ll \left|v_5\right|\ll \left|v_7\right|$E41

hold, one must obtain that the succession function on the center manifold has three small real positive roots, just the system (30) has at least three limit cycles in the neighborhood of the origin. We can refer to references [16, 26, 27] for more details about the construction of limit cycles.

Remark 7. In general, in order to find more limit cycles in the neighborhood of the origin of system (30), we should add more higher order terms of $\tilde{u}(z\mathrm{,}w)$ determined in Eq. (33). Here we propose a conjecture that system (30) has at most three limit cycles through Hopf bifurcation restricted to a center manifold from the origin. However, the center conditions or integrability at the degenerate singularity will need further study.

4. Conclusion and discussion

The two classes of methods for computing the nondegenerate and degenerate singular point quantities on center manifold of the three‐, four‐, and more higher dimensional polynomial systems are discussed here, and more as the applications of them, the multiple limit cycles or Hopf cyclicity of two typical nonlinear dynamic systems restricted to the corresponding center manifolds are investigated.

Appendix A

$\begin{array}{l}c[\alpha \mathrm{,}\beta ]=\\ \frac{1}{5(\alpha -\beta )}{d}_1\mathrm{\{}{b}_{02}^2(3\beta -2\alpha )+a_{20}{b}_{02}(20-\beta -\alpha )-a_{20}^2(20+2\beta -3\alpha ))\\ \times d_{1}c[\alpha -\mathrm{17,}\beta -13]+((a_{11}{b}_{02}+a_{20}{b}_{11})(20-\beta -\alpha )-2b_{02}{b}_{11}(5-3\beta +2\alpha )-\\ 2a_{11}{a}_{20}(15+2\beta -3\alpha ))d_{1}c[\alpha -\mathrm{16,}\beta -14]+((a_{02}{b}_{02}+a_{11}{b}_{11}+a_{20}{b}_{20})(20-\\ \beta -\alpha )-(a_{11}^2+2a_{02}{a}_{20})(10+2\beta -3\alpha )-(b_{11}^2+2b_{02}{b}_{20})(10-3\beta +2\alpha ))d_{1}c[\alpha -\\ \mathrm{15,}\beta -15]+((a_{02}{b}_{11}+a_{11}{b}_{20})(20-\beta -\alpha )-2b_{11}{b}_{20}(15-3\beta +2\alpha )-2a_{02}{a}_{11}(5+\\ 2\beta -3\alpha ))d_{1}c[\alpha -\mathrm{14,}\beta -16]+(a_{02}{b}_{20}(20-\beta -\alpha )-b_{20}^2(20-3\beta +2\alpha )-\\ a_{02}^2(2\beta -3\alpha ))d_{1}c[\alpha -\mathrm{13,}\beta -17]-b_{02}(5+3\beta -2\alpha )+a_{20}(5+2\beta -3\alpha ))ic[\alpha -\\ \mathrm{6,}\beta -4]-(b_{11}(3\beta -2\alpha )+a_{11}(2\beta -3\alpha ))ic[\alpha -\mathrm{5,}\beta -5]\\ +(b_{20}(5-3\beta +2\alpha )+a_{02}(5-2\beta +3\alpha ))ic[\alpha -\mathrm{4,}\beta -6]\end{array}$

$\begin{array}{l}\tilde{\mu }[\alpha ]=-\frac{d_1}{5}\mathrm{\{}(a_{20}^2(\alpha -20)+2a_{20}{b}_{02}(10-\alpha )+b_{02}^2\alpha )d_{1}c[\alpha -\mathrm{17,}\alpha -13]\\ +(2a_{11}{a}_{20}(\alpha -15)-2(a_{11}{b}_{02}+a_{20}{b}_{11})(\alpha -10)+2b_{02}{b}_{11}(\alpha -5))d_{1}c[\alpha -\\ \mathrm{16,}\alpha -14]+((a_{11}^2+2a_{02}{a}_{20}-2a_{02}{b}_{02}-2a_{11}{b}_{11}+b_{11}^2-2a_{20}{b}_{20}+2b_{02}{b}_{20})(\alpha -\\ 10))d_{1}c[\alpha -\mathrm{15,}\alpha -15]+2((a_{02}{b}_{11}+a_{11}{b}_{20})(10-\alpha )-b_{11}{b}_{20}(15-\alpha )-a_{02}{a}_{11}(5-\\ \alpha ))d_{1}c[\alpha -\mathrm{14,}\alpha -16]+(b_{20}^2(\alpha -20)-2(a_{02}{b}_{20})(\alpha -10)+a_{02}^2\alpha )d_{1}c[\alpha -\mathrm{13,}\alpha -\\ 17]+(a_{20}(\alpha -5)-b_{02}(5+\alpha ))ic[\alpha -\mathrm{6,}\alpha -4]+(a_{11}-b_{11})\alpha iic[\alpha -\mathrm{5,}\alpha -\\ 5]-(b_{20}(\alpha -5)-a_{02}(5+\alpha ))ic[\alpha -\mathrm{4,}\alpha -6]\mathrm{\}},\\ {\mu }_m=\tilde{\mu }[5m],\end{array}$

where $c[k\mathrm{,}j]=c_{kj}$.

Acknowledgments

This work was supported by Natural Science Foundation of China grants (11461021, 11261013), Nature Science Foundation of Guangxi (2015GXNSFAA139011), Research Foundation of Hezhou University (No.HZUBS201302), and Guangxi Education Department Key Laboratory of Symbolic Computation and Engineering Data Processing.

MR(2000) Subject Classification: 34C05, 37C07