Lyapunov coefficients for degenerate Hopf bifurcations and an application in a model of competing populations

Accepted Manuscript Lyapunov coefficients for degenerate Hopf bifurcations and an application in a model of competing populations

Jocirei D. Ferreira, Aida P. González Nieva, Wilmer Molina Yepez

PII: DOI: Reference:

S0022-247X(17)30500-0 http://dx.doi.org/10.1016/j.jmaa.2017.05.040 YJMAA 21401

To appear in:

Journal of Mathematical Analysis and Applications

Received date:

30 August 2016

Please cite this article in press as: J.D. Ferreira et al., Lyapunov coefficients for degenerate Hopf bifurcations and an application in a model of competing populations, J. Math. Anal. Appl. (2017), http://dx.doi.org/10.1016/j.jmaa.2017.05.040

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Lyapunov Coefficients for Degenerate Hopf Bifurcations and an Application in a Model of Competing Populations Jocirei D. Ferreiraa,∗∗, Aida P. Gonz´alez Nievab , Wilmer Molina Yepezb a Federal

University of Mato Grosso, Institute of Exact and Earth Science, Barra do Gar¸cas-MT, Brazil b Departmento de Matem´ aticas, Universidad del Cauca, Popay´an, Colombia

Abstract This paper focus on the study of codimension one and two Hopf bifurcations and the pertinent Lyapunov stability coefficients for a general reaction-diffusion system. We display algebraic expressions for the first and second Lyapunov coefficients for the infinite dimensional system subject to Neumann boundary conditions. As an application, a special subset of a three-dimensional Lotka-Volterra dynamical system with diffusions subject to Neumann boundary conditions is analyzed. The main goal is to perform a detailed local stability analysis for the proposed predator-prey model to show the existence of multiple spatially homogeneous and non-homogeneous periodic orbits, due to the occurrence of codimension one Hopf bifurcation. Keywords: Reaction diffusion system, Lyapunov coefficients, degenerate Hopf bifurcation, Predator-prey system, spatially homogeneous bifurcation, spatially non-homogeneous bifurcation.

1. Introduction Mathematical ecology has its roots in population ecology, which treats the increase and fluctuation of population. Since Lotka-Volterra’s pioneering work on population dynamics, numerous mathematical models have been proposed to study the relation between predator and prey populations. This type of models is generally known as predator-prey models. Due to the importance of such models in ecology, predator-prey models have been intensively studied and will continue to be one of the dominant themes in the area. The goal of the Lotka-Volterra model is the quantification of the interaction within and between species and their inorganic environment, and the investigation of the temporal variation of groups of individuals of various species. However, it is true that time and space are inseparable “coordinates,” and only when populations of organisms are considered in both time and space can the ecological situations be understood. As it is well known, the classical predator-prey model reflects population changes due to predation in a situation where predator and prey densities are not spatially dependent. It does not take into account neither the fact that population is not homogeneously distributed nor the fact that predators and preys naturally develop strategies for survival. Both of these considerations involve diffusion processes which can be quite intricate as different concentration levels of preys and predators cause different population movements. The role of diffusion in the modelling of many physical, chemical and biological processes has been extensively studied. Starting with Turing’s paper [4], diffusion have been observed as causes of the spontaneous emergence of ordered structures, called patterns, in a variety of non-equilibrium situations. It must be realized that spatial ecology addresses the fundamental effects of space on the dynamics of individual species and on the structure, dynamics, diversity, and stability of multispecies communities. Essentially, this subject is designed to highlight the importance ∗ The ordering of the authors is alphabetical since all authors contributed equally to the planning and execution of this research, and the writing of the paper. ∗∗ Corresponding author. Email addresses: [email protected] (Jocirei D. Ferreira), [email protected] ( Aida P. Gonz´alez Nieva), [email protected] (Wilmer Molina Yepez)

Preprint submitted to Elsevier

May 22, 2017

of space in the areas of stability, patterns of diversity, invasions, coexistence, and pattern generation. The abovementioned patterns (fronts, spirals, targets, hexagons, stripes, and dissipative solitons) can be found in various types of reaction-diffusion systems in spite of large discrepancies in the local reaction terms. It has also been argued that reaction-diffusion processes are an essential basis for those connected with animal coats and skin pigmentation [14, 15]. In the case that the population changes are not spatially dependent it is tacitly assumed that members of the same or different populations may meet any other member with equal probability. This assumption plays an important role in problems concerning reproduction and also in predator-prey and other relations. This means that populations or ecological systems are considered either to be concentrated in a point of space or they are “well stirred”. This assumption is acceptable if the system is that of a small lake, a small wood, or a test tube in a laboratory. However, one generally has to take into consideration that populations must exploit a smaller or larger territory on the surface of the Earth or the sea, and the distance between the members may count in possible interactions. Thus, it is more realist to consider models which describe the dynamics of populations that live in a 1, 2 or 3 dimensional space, for instance fish in the sea or in a lake, animals on the savannah, or worms moving up and down on a single tree. In these cases the dynamics is described, for example, by a system of reaction-diffusion equations, that is, a nonlinear or rather quasilinear system of parabolic partial differential equations. All beings, including different kinds of populations live in a spatial world and it is a natural phenomenon that a substance goes from high-density regions to low-density regions. As a result, more and more scholars use spatial artifacts to model the interactions between the prey and the predator species. In a natural setting when populations diffuse/migrate from one environment to another, they are often subject to unfavorable conditions. Naturally, predation will obligate the individuals of populations to follow defense strategies as, for example, to diffuse in the direction of lower concentrations of its predators. As a consequence, the predators also will develop some strategies to take advantage of the defense switching of the prey, looking at the migration of the prey and diffusing in that direction. For more interesting account of various aspects and examples in population ecology and more details of the spatial variation of groups of individuals of various species we refer the readers to [5, 7]. In response to this need of considering and analyzing temporal and spatial variations in mathematical models describing the dynamics interaction among groups of individuals, as a start point in this exposition we analyze a general reaction-diffusion system consisting of three equations and subject to Neumann boundary conditions on spatial domain Ω = (0, π), with  ∈ R+ . Our analyzes focus on the study of degenerate Hopf bifurcation and the pertinent Lyapunov stability coefficients, displaying expressions for the first and second Lyapunov coefficients. The expressions for these coefficients will be of interest in the study of the codimension one and two Hopf bifurcation for a class of reaction-diffusion systems subject to Neumann boundary conditions modeling real world problems. It is important to observe that a similar analyzes was carried out in [17] for n−dimensional ODE systems, where the authors displayed expressions for the first, second, third and fourth Lyapunov coefficients. Following these ideas, we carry out a review of the method found in [6] (see also [11]) for the calculation of the first Lyapunov coefficient for the nonlinear autonomous reaction-diffusion system. It is worth observing that the procedure considered in this paper to get the second Lyapunov coefficient requiring no explicit formal evaluation of the center manifold, has not been found by the authors in the current literature and is better adapted to the needs of long calculations of the Lyapunov coefficients. Furthermore, this method can be adapted to obtain the expressions for the third and fourth Lyapunov coefficients with few modifications. Finally, we observe that a similar analysis for the calculation of the Lyapunov coefficients can be performed by considering a domain Ω ⊂ RN (N = 1, 2, 3), which is interesting to study the effects of high dimensional domain in reaction-diffusion models since this consideration provides scope for the development of increasingly realistic models. We defer our work in this direction to a subsequent exposition since this analysis is beyond the scope of the current paper. As an application of the Hopf bifurcation formulas, we investigate the occurrence of codimension 1 Hopf bifurcation in a diffusive predator-prey system represented by the following PDE system ut = d1 u xx + γ1 [1 + K − u − Kw]u, vt = d2 v xx + γ2 [u − v]v, wt = d3 w xx + γ3 [v − w]w,

2

x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0.

(1)

We will assume that the functions u, v and w satisfy the following Neumann boundary conditions and initial conditions u x (0, t) = v x (0, t) = w x (0, t) = 0 u x (π, t) = v x (π, t) = w x (π, t) = 0 u (x, 0) = u0 (x) , v (x, 0) = v0 (x) , w (x, 0) = w0 (x)

t > 0, t > 0, x ∈ (0, π) ,

(2)

where Ω = (0, π) is an open connected bounded interval of R with  ∈ R+ ; In system (1), u = u(x, t), v = v(x, t) and w = w(x, t) represent the population densities; d1 , d2 , d3 > 0 are the diffusion rates; the parameters γ1 , γ2 , γ3 and K are assumed to be strictly positive, where the parameter K will be considered as our bifurcation parameter. It is worth noting that the related ordinary differential equation system of model (1)-(2), which describes the temporal variation of species u˙ = γ1 [1 + K − u − Kw]u, v˙ = γ2 [u − v]v, w˙ = γ3 [v − w]w,

x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0,

(3)

was considered in [8] (see also [26, 27]). The authors carried out a bifurcation analysis and investigated the existence and stability of limit cycles for system (3). We are studying a special subset of three-dimensional Lotka-Volterra dynamical system with diffusion which we can interpret biologically as follow. In the system (1)-(2), the first equation shows that w is a predator and u is a prey. In the second equation , v is a predator and u is a prey. Finally, in the third equation, w is a predator and v is a prey. We are particularly interested in the study of the stability of isolated equilibrium populations, wherein we assume that the prey and predators are diffusing in the domain (0, π) ⊂ R. So as to ensure that the equilibrium of system (1)-(2) is isolated and corresponds to that of the related ODE model we invoke the Neumann boundary conditions. Other studies for this kind of model can be find in [16, 19, 20, 25]. Our main purpose in the application presented in this exposition is to generalize the result given in [8] for the reaction diffusion system (1)-(2), that is, we show that the populations coexist in a periodic manner due to the occurrence of an Andronov-Hopf bifurcation even when diffusion is considered in the system (3). It is important to observe that, to the best of our knowledge there is no result in the existing literature regarding to the direction of the Hopf bifurcation and stability of the bifurcating periodic solutions for the system (1)-(2). Pursuing the bifurcation analysis carried out in [8] for the system (3), we perform a detailed bifurcation analysis for the constant coexistence equilibrium solution of system (1)-(2), on spatial domain Ω = (0, π), with  ∈ R+ . We show that under certain conditions on the parameters values, the system (1)-(2) undergoes a codimension 1 Hopf bifurcation and the bifurcating periodic orbits are either spatially homogeneous, which coincides with the periodic solution of the corresponding ODE system, or spatially non-homogeneous. Furthermore, we show that there exist infinity Hopf bifurcation points where spatially non-homogeneous periodic orbits bifurcate. These periodic orbits correspond to the spatial eigen-modes cos( jx) ( j ≥ 1) where π is the length of the spatial domain. The present paper is organized as follows. In Section 2 we present a review of the method found in [6] (see also [11]) for the calculation of the first Lyapunov coefficient for a general reaction-diffusion systems. Furthermore, we present the procedure to calculate the expression for the second Lyapunov coefficient, pushing forward the method found in [6, 11]. In Section 3, we defer a study of its homogeneous equilibrium points (or homogeneous solutions) and investigate their asymptotic behavior. Furthermore, as an application of the Hopf bifurcation formulas, we consider the related codimension 1 Hopf bifurcation for the system (1)-(2) with the spatial domain Ω = (0, π),  ∈ R+ . Specifically, we show the occurrence of spatially homogeneous and non-homogeneous codimension 1 Hopf bifurcation for the system (1)-(2). In section 4, we present some simulations carried in MATLAB to verify the veracity of the analytical results obtained for system (1)-(2) related to the non-homogeneous Hopf bifurcation. A discussion follows in the final section. 2. Explicit expressions for the first and second Lyapunov coefficients In this section we present a review of the method described in [6] (see also [11]) for the calculation of the first Lyapunov coefficient for a general reaction-diffusion system subject to Neumann boundary conditions. The calculation of the second Lyapunov coefficient has not been found by the authors in the current literature. Thus, we display the 3

expression for the second Lyapunov coefficient pushing forward the method presented in [6, 11]. It is worth observing that the expressions for the first and second Lyapunov coefficients, denoted by l1 and l2 respectively, will be of interest in the study of the codimension one and two Hopf bifurcation for a whole class of reaction-diffusion systems subject to Neumann boundary conditions modeling real world problems as, for example, in the study of models describing the dynamics interaction within and between species and their environment. Accordingly, we consider the reaction diffusion model represented by the following PDE system ⎧ ⎪ x ∈ Ω, t > 0, u = d1 Δu + f (K, u, v, w) ⎪ ⎪ ⎨ t v = d Δv + g(K, u, v, w) x ∈ Ω, t > 0, (4) ⎪ t 2 ⎪ ⎪ ⎩ w = d Δw + h(K, u, v, w) x ∈ Ω, t > 0, t

3

subject to Neumann boundary conditions ∂u ∂ν

(x, t) =

∂v ∂ν

(x, t) =

∂w ∂ν

(x, t) = 0 on ∂Ω × [0, ∞),

(5)

and initial conditions u (x, 0) = u0 (x) , v (x, 0) = v0 (x) , w (x, 0) = w0 (x) on Ω,

(6)

 where Ω ⊂ RN (N ≥ 1), ν = ν(x) denotes the outer unit normal to ∂Ω, Δ = Nj=1 ∂2 /∂x j 2 is the Laplacian operator, + di > 0, i = 1, 2, 3 are the diffusion rates and K ∈ R is a parameter. As a starting point in this exposition, we considered a one dimensional domain Ω = (0, π), with  ∈ R+ ; the functions f, g, h : R × R3 → R are C k (k ≥ 5) with f (K, 0, 0, 0) = g(K, 0, 0, 0) = h(K, 0, 0, 0) = 0. We define the real-valued Sobolev space   X := (u, v, w) ∈ H 2 (0, π) × H 2 (0, π) × H 2 (0, π)/(u x , v x , w x )| x=0,π = 0 and its complexification XC := X ⊕ iX = {x1 + ix2 : x1 , x2 ∈ X}. If we consider the Taylor series expansion of the functions f , g and h at (K, 0, 0, 0), the linearized system associated with system (4) at (K, 0, 0, 0) is given by ⎞ ⎛ ⎛ ⎜⎜⎜ ut ⎟⎟⎟ ⎜⎜⎜ u ⎟ ⎜⎜⎜ ⎜ ⎜⎜⎝ vt ⎟⎟⎟⎟⎠ = L(K) ⎜⎜⎜⎝⎜ v wt w

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ,

where

⎛ ∂2 ⎜⎜⎜ d1 ∂x 2 + A(K) ⎜⎜⎜ L(K) := ⎜⎜⎜ D(K) ⎝ G(K)

B(K) ∂ d2 ∂x 2 + E(K) H(K) 2

C(K) F(K) ∂2 d3 ∂x 2 + I(K)

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ ⎠

(7)

is the linearized operator with the domain DL(K) = XC and A(K) = fu (K, 0), D(K) = gu (K, 0), G(K) = hu (K, 0),

B(K) = fv (K, 0), C(K) = fw (K, 0), E(K) = gv (K, 0), F(K) = gw (K, 0), H(K) = hv (K, 0), I(K) = hw (K, 0).

It is well known that the following hypothesis is a necessary condition for system (4)-(6) to exhibits Hopf bifurcation. Hopf hypothesis: There exists a neighborhood B(K0 , δ) of K0 such that for K ∈ B(K0 , δ), L(K) has a pair of complex, simple, conjugate eigenvalues Γ(K) = α(K) ± iω(K), continuously differentiable in K, such that α(K0 ) = 0, ω(K0 ) = ω0 > 0, and α (K0 )  0; Furthermore, all other eigenvalues of L(K) have no null real parts for K ∈ B(K0 , δ). Hence, to carry out a Hopf bifurcation analysis for the system (4)-(6), we study the eigenvalues of the differential operator L(K). It is well known that the eigenvalue problem  −ϕ = μϕ; x ∈ (0, π), (8) ϕ (0) = ϕ (π) = 0. 4

has eigenvalues μm =

m2 , 2

m = 0, 1, 2, 3, . . ., with corresponding eigenfunctions ϕm (x) = cos ⎛ ⎛ ⎞ ⎜⎜⎜ ϕ ⎟⎟⎟  ⎜ a ∞ mx ⎜⎜⎜⎜ m ⎜⎜⎜ ⎟⎟⎟ ⎜⎜ bm cos ⎜⎜⎝ ψ ⎟⎟⎠ =  ⎜⎝ c m=0 φ m



mx 



. Let

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

(9)

be an eigenfunction of L(K) with eigenvalue β(K), that is L(K)(ϕ, ψ, φ)T = β(K)(ϕ, ψ, φ)T . The following proposition follows immediately by taking into account (9) and therefore the proof is omitted. ⎞ ⎛ ⎜⎜⎜ am ⎟⎟⎟ ⎟ ⎜ Proposition 1. For each m = 0, 1, 2, 3, . . . we have that ⎜⎜⎜⎜ bm ⎟⎟⎟⎟ is an eigenvector of the operator Lm (K) defined by ⎠ ⎝ cm ⎛ 2 ⎜⎜⎜ A(K) − d1m2 ⎜⎜⎜ Lm (K) = ⎜⎜⎜ D(K) ⎜⎝ G(K)

B(K) 2 E(K) − d2m2 H(K)

CK) F(K) 2 I(K) − d3m2

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ , ⎟⎠

with eigenvalue β(K), that is ⎛ ⎜⎜⎜ am ⎜ Lm (K) ⎜⎜⎜⎜ bm ⎝ cm

⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ am ⎟⎟⎟ ⎜ ⎟⎟⎠ = β(K) ⎜⎜⎜⎜⎝ bm cm

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ , m = 0, 1, 2, 3, . . .

where am , bm and cm are the Fourier coefficients of ϕ, ψ and φ, respectively. The Proposition 1 states that the eigenvalues of L(K) are given by the eigenvalues of Lm (K) for m = 0, 1, 2, . . .. Therefore, we consider the characteristic equation of Lm (K) which is given by −λ3 + Mm λ2 + Nm λ + Pm = 0, m = 0, 1, 2, ..., where 2

Mm (K) = Nm (K) = Pm (K)

(E + I + A) − (d1 + d2 + d3 ) m2 − (BD + CG + FH) − (AI + EI + EA) 2 4 + (Ed1 + Id1 + Ad2 + Id2 + Ed3 + Ad3 ) m2 − (d1 d2 + d1 d3 + d2 d3 ) m4 = − (AFH + BDI + CEG) + (BFG + CDH + AEI) 2 + [(FH − EI) d1 + (CG − IA) d2 + (BD − EA) d3 ] m2 4 6 + (Ed3 d1 + Ad2 d3 + Id1 d2 ) m4 − d3 d1 d2 m6

Following the framework of [6] (see also [11]), the system (4) can be written in the abstract form dU = L (K) U + F (K, U) dt

(10)

with U = (u, v, w)T ∈ X and ⎛ ⎜⎜⎜ f (K, u, v, w) − A(K)u − B(K)v − C(K)w ⎜ F (K, U) = ⎜⎜⎜⎜ g (K, u, v, w) − D(K)u − E(K)v − F(K)w ⎝ h (K, u, v, w) − G(K)u − H(K)v − I(K)w

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ;

in particular, at K = K0 the system (10) takes the form dU = L (K0 ) U + F0 (U) where F0 (U) = F (K0 , U) . dt 5

(11)

We introduce on XC the complex-valued L2 (0, π) inner product defined by π

U1 , U2  =

(¯u1 u2 + v¯ 1 v2 + w¯ 1 w2 ) dx

(12)

0

with Ui = (ui , vi , wi ) ∈ XC (i = 1, 2) . Furthermore, we denote the adjoint operator of L (K) by L∗ (K) which satisfy T

U, L(K)V = L∗ (K)U, V for all U, V ∈ XC , with DL∗ (K) = XC .

(13)

The next proposition follows immediately using (13) and the proof will be omitted. Proposition 2. Let U, V be elements of XC and L (K0 ) the differential operator associated with the system (4) at K = K0 . Then the adjoint operator of L (K0 ) is given by ⎛ ∂2 ⎜⎜⎜ d1 ∂x D(K0 ) G(K0 ) 2 + A(K0 ) ⎜⎜⎜ ∂2 L (K0 ) = ⎜⎜⎜ B(K0 ) d2 ∂x2 + E(K0 ) H(K0 ) ⎝ ∂2 C(K0 ) F(K0 ) d3 ∂x 2 + I(K0 ) ∗

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ . ⎠

∗ (K), m = 1, 2, 3, ..., with eigenvalue As in Proposition 1 it can be proved that (a∗m , b∗m , c∗m )T is an eigenvector of Lm β(K) where

⎛ 2 ⎜⎜⎜ A(K) − d1m2 ⎜ ⎜ ∗ (K) := ⎜⎜⎜⎜ Lm B(K) ⎜⎝ C(K)

D(K) 2 E(K) − d2m2 F(K)

G(K) H(K) 2 I(K) − d3m2

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ . ⎟⎠

If the Hopf hypothesis holds at K = K0 then L (K0 ) has a pair of simple purely imaginary eigenvalues ±iω0 if and only if there exists a unique m ∈ N such that ±iω0 are purely imaginary eigenvalues of Lm (K0 ). We denote the T associated eigenfunction by q = cos mx  (am , bm , cm ) , with am , bm , cm ∈ C, such that L(K0 )q = iω0 q. Furthermore, we m ∗ ∗ ∗ ∗ T can choose q = cos  x(am , bm , cm ) ∈ XC as an eigenfunction of L∗ (K0 ) such that ¯ = 0 and q∗ , q = 1. L∗ (K0 ) q∗ = −iω0 q∗ , q∗ , q We introduce the spaces X c := {zq + zq : z ∈ C} and X su := {u ∈ X : q∗ , u = 0} and, for sake of completeness, we prove the following lemma that allows to decompose the space X as X = X c ⊕ X su . Lemma 1. Let X c be the generalized eigenspace of L(K) corresponding to the eigenvalues with null real parts and X su the space of functions u ∈ X such that < q∗ , u >= 0. The space X can be decomposed as X = X c ⊕ X su . ¯ z ∈ C, and < U, q∗ >= 0 i.e, Proof. It is easy to see that X c ∩ X su = {0}. In fact, if U ∈ X c ∩ X su then U = zq + z¯q, ∗ ∗ ∗ < zq + z¯q, ¯ q >= z < q, q > +¯z < q, ¯ q >= 0. So, z = 0 which implies U = 0. To show that X c ⊕ X su = X observe that if q ∈ XC = X + iX, then q = q1 + iq2 where q1 , q2 ∈ X. Furthermore, zq + z¯q¯ = z(q1 + iq2 ) + z¯(q1 − iq2 ) = q1 (z + z¯) + iq2 (z − z¯) = 2q1 Re(z) − 2q2 Im(z) ∈ X. So, if U ∈ X c ⊕ X su then U = zq + z¯q¯ + u where zq + z¯q¯ ∈ X c and u ∈ X su . From the previous observation it follows that zq + z¯q¯ ∈ X and from the definition of X su we have u ∈ X. Hence, U ∈ X. On the other hand, if U ∈ X and z =< q∗ , U > then < q∗ , U − zq − z¯q¯ >=< q∗ , U > −z = 0. Therefore,there exists z ∈ C such that U − zq − z¯q¯ ∈ X su . Furthermore, since U = (zq + z¯q) ¯ + (U − zq − z¯q) ¯ ∈ X c ⊕ X su it follows that X ⊂ X c ⊕ X su . Remark 1. The result in Lemma 1 can be obtained if we replace X by XC , which allows to decompose the space XC as XC = X c ⊕ X su . 6

From Lemma 1, if U = (u, v, w)T ∈ X, there exists r = (r1 , r2 , r3 )T ∈ X su and z ∈ C such that ⎛ ⎞ ⎞ ⎛ ⎧ mx ⎪ ⎜⎜⎜ u ⎟⎟⎟ ⎜⎜ r1 ⎟⎟ ¯ m cos mx ⎪  + z¯a  + r1 , ⎪ u = zam cos mx ⎜⎜⎜⎜ v ⎟⎟⎟⎟ = zq + z¯q¯ + ⎜⎜⎜⎜⎜ r2 ⎟⎟⎟⎟⎟ , or equivalently ⎨ ¯ m cos mx + r2 , cos + z ¯ b v = zb ⎪ m   ⎪ ⎜⎝ ⎟⎠ ⎟⎠ ⎜⎝ ⎪ ⎩ w = zc cos mx r3 ¯ m cos mx w m  + z¯c  + r3 . Proposition 3. In the coordinate of (14) the variables z and r can be expressed as  z = q∗ , U ¯ r = U − q∗ , Uq − q¯ ∗ , Uq. Furthermore, in the coordinates of (15) the system (11) assumes the form  dz ∗ (U) dt = iω0 z + q , F 0 dr ) (K r + H(z, z¯, r), = L 0 dt where

¯ H(z, z¯, r) = F0 (U) − q∗ , F0 (U)q − q¯ ∗ , F0 (U)q.

(14)

(15)

(16)

(17)

Remark 2. Observe that H(z, z¯, r) is a smooth function of z, z¯ ∈ C1 and r ∈ C3 . Moreover, the Taylor expansion of H (z, z¯, r) around (0, 0, r) is given by  1 (0, 0, r) zk z¯ j H (z, z¯, r) = k! j! Hk j =

k+ j≥2 H20 2 z + H202 z¯2 + H630 z3 + H221 z2 z¯ + H212 z¯z2 + H603 z¯3 + H2440 z4 + H631 z3 z¯ 2 z + H11 z¯ H13 3 H50 5 H22 2 2 + 4 z z¯ + 6 z¯z + H2404 z¯4 + 120 z + H2441 z4 z¯ + H1232 z3 z¯2 + H1223 z2 z¯3 + H2414 z¯z4 H05 5 6 + 120 z¯ + o(|z| ) + o(|z| |r|),

(18)

where Hk j (0, 0, r) =

∂k+ j ∂zk ∂z j

H(z, z, r)|z=0 ,

and

q∗ , F0 (U)

=

1 2 z + 12 G02 z¯2 + 3!1 G30 z3 + 12 G21 z2 z¯ + 12 G12 z¯z2 2 G 20 z + G 11 z¯ 1 1 3 + 3! G03 z¯ + 4! G40 z4 + 3!1 G31 z3 z¯ + 14 G22 z2 z¯2 + 3!1 G13 z¯z3 + 4!1 G04 z¯4 1 1 + 5!1 G50 z5 + 4!1 G41 z4 z¯ + 12 G32 z3 z¯2 + 12 G23 z2 z¯3 + 4!1 G14 z¯z4 + 5!1 G05 z¯5 1 2 + G10 , r z + G01 , r z¯ + 2 G110 , r z + 12 [ G210 , r + G120 , r] z¯z + 12 G220 , r z¯2 + 16 z3 G1110 , r + 16 z2 z¯ [ G1120 , r + G1210 , r + G2110 , r] + 16 z¯z2 [ G1220 , r + G2120 , r + G2210 , r] + 16 z¯3 G2220 , r 1 4 1 3 z G11110 , r + 24 z z¯ [ G11120 , r + G11210 , r + G12110 , r + G21110 , r] + 24 1 2 2 z z¯  [ G11220 , r + G12120 , r + G21120 , r + G12210 , r + G 21210 , r + G22110 , r] + 24 1 1 1 1 4 z¯z3 G12220 , r + 24 z¯ G22220 , r + 24

G21220 , r + 24

G22120 , r + G22210 , r + 24 1 5 + 120 z G111110 , r 1 4 z z¯ [ G111120 , r + G111210 , r + G112110 , r + G121110 , r + G211110 , r] + 120 1 3 2 z z¯ [ G111220 , r + G112210 , r + G121210 , r + G211210 , r + G112120 , r] + 120 1 3 2 z z¯ [ G121120 , r + G211120 , r + G122110 , r + G221110 , r + G212110 , r] + 120 1 2 3 z z¯ [ G112220 , r + G121220 , r + G122120 , r + G122210 , r + G211220 , r] + 120 1 2 3 z z¯ [ G212120 , r + G212210 , r + G221120 , r + G222110 , r + G221210 , r] + 120 1 z¯z4 [ G122220 , r + G212220 , r + G221220 , r + G222120 , r + G222210 , r] + 120 1 5 + 120 z¯ G222220 , r + · · ·

(19)

with Gi j ∈ C1 , 2 ≤ i + j ≤ 5, and the standard scalar product in L2 (0, π) is considered. The scalar product now means 3  G¯ i ri . Complex numbers and vectors involved in (19) can be computed by the formulas given in that < G, r >= i=1

Appendix 1, which are not presented here since these expressions are too long. 7

Assume that the function F0 (U) = F(K0 , U) can be written as F0 (U) =

1 2 B(U, U)

+ 16 C (U, U, U) +

1 24 D(U, U, U, U)

where BXY = B(X, Y), CXYZ = C(X, Y, Z), DXYZW = multilinear functions with ⎞ ⎞ ⎛ ⎛ ⎛ ⎜⎜⎜ x1 ⎟⎟⎟ ⎜⎜⎜ y1 ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜ ⎜⎜⎜ ⎜ X = ⎜⎜ x2 ⎟⎟ , Y = ⎜⎜ y2 ⎟⎟ , Z = ⎜⎜⎜⎜ ⎠ ⎠ ⎝ ⎝ ⎝ x3 y3

+

1 120 E(U, U, U, U, U)

  + o U6 ,

(20)

D(X, Y, Z, W) and EXYZWV = E(X, Y, Z, W, V) are symmetric z1 z2 z3

⎞ ⎞ ⎛ ⎛ ⎟⎟⎟ ⎜⎜⎜ w1 ⎟⎟⎟ ⎜⎜⎜ v1 ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜ ⎜ ⎟⎟⎠ , W = ⎜⎜⎝ w2 ⎟⎟⎠ , V = ⎜⎜⎜⎜⎝ v2 w3 v3

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

In coordinates, we have Bi (X, Y) =

 3  2  ∂ Fi (K0 , η) j,k=1

Di (X, Y, Z, W) =

∂η j ∂ηk

η=0

x j yk , Ci (X, Y, Z) =

 4  3  ∂ Fi (K0 , η) x j y k z l wr , ∂η j ∂ηk ∂ηl ∂ηr η=0 j,k,l,r=1

 3  3  ∂ Fi (K0 , η) x j yk z l , ∂η j ∂ηk ηl η=0 j,k,l=1

Ei (X, Y, Z, W, V) =

 3  j,k,l,r,p=1

∂5 Fi (K0 , η) ∂η j ∂ηk ∂ηl ∂ηr ∂η p

(21)  η=0

x j yk zl wr v p , (22)

where i = 1, 2, 3, and F1 = f, F2 = g, F3 = h. The next results will be of interest for later uses on the calculation of the first and second Lyapunov coefficients. Proposition 4. From (17), (18) and (20) the coefficients of the Taylor expansion of the function H (z, z¯, r), with all the partial derivatives appearing hereafter evaluated at (K0 , 0, 0, 0), are given by ⎧ ⎪ H31 = Dqqqq¯ − q∗ , Dqqqq¯ q − q¯ ∗ , Dqqqq¯ q¯ H20 = Bqq − q∗ , Bqq q − q¯ ∗ , Bqq q¯ ⎪ ⎪ ⎪ ⎪ ∗ ∗ ⎪ H13 = Dqq¯ q¯ q¯ − q∗ , Dqq¯ q¯ q¯ q − q¯ ∗ , Dqq¯ q¯ q¯ q¯ H22 = Dqqq¯ q¯ − q , Dqqq¯ q¯ q − q¯ , Dqqq¯ q¯ q¯ ⎪ ⎪ ⎪ ⎪ ∗ ∗ ⎪ ⎪ H04 = Dq¯ q¯ q¯ q¯ − q∗ , Dq¯ q¯ q¯ q¯ q − q¯ ∗ , Dq¯ q¯ q¯ q¯ q¯ H11 = Bqq¯ − q , Bqq¯ q − q¯ , Bqq¯ q¯ ⎪ ⎪ ⎪ ∗ ∗ ⎪ ⎪ = B −

q , B q −

q ¯ , B  q ¯ H50 = Eqqqqq − q∗ , Eqqqqq q − q¯ ∗ , Eqqqqq q¯ H q¯ q¯ q¯ q¯ q¯ q¯ ⎪ ⎪ ⎨ 02 ∗ ∗ = C −

q , C q −

q ¯ , C  q ¯ H41 = Eqqqqq¯ − q∗ , Eqqqqq¯ q − q¯ ∗ , Eqqqqq¯ q¯ H (23) ⎪ 30 qqq qqq qqq ⎪ ⎪ ⎪ ∗ ∗ ∗ ∗ ⎪ = C −

q , C q −

q ¯ , C  q ¯ H = E −

q , E q −

q ¯ , E  q ¯ H ⎪ 21 qq q ¯ qq q ¯ qq q ¯ 32 qqq q ¯ q ¯ qqq q ¯ q ¯ qqq q ¯ q ¯ ⎪ ⎪ ⎪ ⎪ H23 = Eqqq¯ q¯ q¯ − q∗ , Eqqq¯ q¯ q¯ q − q¯ ∗ , Eqqq¯ q¯ q¯ q¯ H12 = Cqq¯ q¯ − q∗ , Cqq¯ q¯ q − q¯ ∗ , Cqq¯ q¯ q¯ ⎪ ⎪ ⎪ ⎪ ∗ ∗ ⎪ ⎪ H14 = Eqq¯ q¯ q¯ q¯ − q∗ , Eqq¯ q¯ q¯ q¯ q − q¯ ∗ , Eq¯ q¯ q¯ q¯ q¯ H03 = Cq¯ q¯ q¯ − q , Cq¯ q¯ q¯ q − q¯ , Cq¯ q¯ q¯ q¯ ⎪ ⎪ ⎪ ⎩ H40 = Dqqqq − q∗ , Dqqqq q − q¯ ∗ , Dqqqq q¯ H05 = Eq¯ q¯ q¯ q¯ q¯ − q∗ , Eq¯ q¯ q¯ q¯ q¯ q − q¯ ∗ , Eq¯ q¯ q¯ q¯ q¯ q. ¯ Proof. The result follows by plugging U = zq + z¯q¯ + r in F0 (U) given by (20), replacing the obtained result in (17) and finally equating coefficients of like power in z and z¯ between the right hand side of (17) and the right hand side of (18), after respective substitutions. Remark 3. Since the expressions for Bqq , Bqq¯ , Bq¯ q¯ , Cqqq , Cqqq¯ , Cqq¯ q¯ , Cq¯ q¯ q¯ , Dqqqq , Dqqqq¯ , Dqqq¯ q¯ , Dqq¯ q¯ q¯ , Dq¯ q¯ q¯ q¯ , Eqqqqq , Eqqqqq¯ , Eqqqq¯ q¯ , Eqqq¯ q¯ q¯ , Eqq¯ q¯ q¯ q¯ and Eq¯ q¯ q¯ q¯ q¯ given in Proposition 4 are too long, we will present them in the Appendix 2. From appendix A of [6] the system (16) has a center manifold which can be represented locally as the graphic of a smooth function W c (0) = {(z, z, r) : r = V(z, z)} with V : C2 → X su . Thus, we can write r in the form r=

 2≤ j+k≤5

1 r jk z j z¯k + o(|z|6 ) j!k!

8

(24)

or r = V(z, z)

=

r20 2 z + r202 z¯2 + 3!1 r30 z3 + 12 r21 z2 z¯ + 12 r12 z¯z2 + 3!1 r03 z¯3 2 z + r11 z¯ 1 + 4!1 r40 z4 + 3!1 r31 z3 z¯ + 2!2! r22 z2 z¯2 + 3!1 r13 z¯z3 + 4!1 r04 z¯4   1 1 1 1 5 4 + 5! r50 z + 4! r41 z z¯ + 3!2! r32 z3 z¯2 + 2!3! r23 z2 z¯3 + 4!1 r14 z¯z4 + 5!1 r05 z¯5 + o |z|6

(25)

  with q∗ , ri j = 0 and r jk = r¯k j ∈ C3 , 2 ≤ j + k ≤ 5.

Lemma 2. From (17), (18), (25) and observing that dr ∂r dz ∂r d¯z = + = L(K0 )r + H(z, z¯, r) dt ∂z dt ∂¯z dt we obtain

    Bqq − q∗ , Bqq q − q¯ ∗ , Bqq q¯     Bqq¯ − q∗ , Bqq¯ q − q¯ ∗ , Bqq¯ q¯     ¯ Bq¯ q¯ − q∗ , Bq¯ q¯ q − q¯ ∗ , Bq¯ q¯ q,

(27)

    −3r20 q∗ , Bqq − 3r11 q¯ ∗ , Bqq             +3 Bqr20 − q∗ , Bqr20 q − q¯ ∗ , Bqr20 q¯ + Cqqq − q∗ , Cqqq q − q¯ ∗ , Cqqq q¯         −2r20 q∗ , Bqq¯ − r11 q∗ , Bqq − r02 q¯ ∗ , Bqq − 2r11 q¯ ∗ , Bqq¯             + Bqr − q∗ , Bqr ¯ ∗ , Bqr ¯ + 2 Bqr11 − q∗ , Bqr11 q − q¯ ∗ , Bqr11 q¯ ¯ 20 q − q ¯ 20 q      ¯ 20  + Cqqq¯ − q∗ , Cqqq¯ q − q¯ ∗ , Cqqq¯ q¯         −r20 q∗ , Bq¯ q¯ − 2r11 q∗ , Bqq¯ − 2r02 q¯ ∗ , Bqq¯ − r11 q¯ ∗ , Bq¯ q¯             ∗ ¯ ∗ , Bqr ¯ + Bqr02 − q∗ , Bqr02 q − q¯ ∗ , Bqr02 q¯ + 2 Bqr ¯ 11 − q , Bqr ¯ 11 q − q ¯ 11 q       ∗ ∗ + Cqq¯ q¯ − q , Cqq¯ q¯ q − q¯ , Cqq¯ q¯ q¯     −3r02 q¯ ∗ , Bq¯ q¯ − 3r11 q∗ , Bq¯ q¯             ∗ ¯ ∗ , Bqr ¯ + Cq¯ q¯ q¯ − q∗ , Cq¯ q¯ q¯ q − q¯ ∗ , Cq¯ q¯ q¯ q¯ , +3 Bqr ¯ 02 − q , Bqr ¯ 02 q − q ¯ 02 q

(28)

− [L (K0 ) − 2iω0 I] r20 −L (K0 ) r11 − [L (K0 ) + 2iω0 I] r02 − [L (K0 ) − 3iω0 I] r30

=

− [L (K0 ) − iω0 I] r21

=

− [L (K0 ) + iω0 I] r12

=

− [L (K0 ) + 3iω0 I] r03

=

− [L (K0 ) − 2iω0 I] r31

−L (K0 ) r22

=

=

(26)

= = =

        −6r20 q∗ , Bqr11 − 6r11 q¯ ∗ , Bqr11 − 3r11 q∗ , Bqr20 − 3r02 q¯ ∗ , Bqr20           ∗ ∗ ∗ ¯ ∗ , Bqq −3r11 q¯ ∗ , Bqr ¯ 20 − 3r20 q , Bqr ¯ 20 − 3r21 q , Bqq − 3r30 q , Bqq¯ − 3r12 q           −3r21 q¯ ∗ , Bqq¯ − r11 q∗ , Cqqq − 3r20 q∗ , Cqqq¯ − r02 q¯ ∗ , Cqqq − 3r11 q¯ ∗ , Cqqq¯             3 Br20 r11 − q∗ , Br20 r11 q − q¯ ∗ , Br20 r11 q¯ + 3 Bqr21 − q∗ , Bqr21 q − q¯ ∗ , Bqr21 q¯             ∗ ¯ ∗ , Bqr ¯ + 3 Cqr20 r11 − q∗ , Cqr20 r11 q − q¯ ∗ , Cqr20 r11 q¯ + Bqr ¯ 30 − q , Bqr ¯ 30 q − q ¯ 30 q       ¯ ∗ , Cqr ¯ − q∗ , Cqr +3 Cqr ¯ 20 r11 q − q ¯ 20 r11 q      ¯ 20 r11  +3 Cqqr11 − q∗ , Cqqr11 q − q¯ ∗ , Cqqr11 q¯             ∗ ¯ ∗ , Cqqr ¯ + Dqqqq¯ − q∗ , Dqqqq¯ q − q¯ ∗ , Dqqqq¯ q¯ , +3 Cqqr ¯ 20 − q , Cqqr ¯ 20 q − q ¯ 20 q

        ¯ ∗ , Bqr −2r11 q¯ ∗ , Bqr02 − 4r02 q¯ ∗ , Bqr11 − 2r02 q¯ ∗ , Bqr ¯ 20 − 4r11 q ¯ 11         − 4r q∗ , B ¯ 11 − 2r20 q∗ , Bqr02 −4r11 q∗ , Bqr11 − 2r11 q∗ , Bqr     ¯ 20  20  qr       −r12 q∗ , Bqq − r03 q¯ ∗ , Bqq − 4r21 q∗ , Bqq¯ − r30 q∗ , Bq¯ q¯ − 4r12 q¯ ∗ , Bqq¯ − r21 q¯ ∗ , Bq¯ q¯         −2r11 q∗ , Cqqq¯ − 2r20 q∗ , Cqq¯ q¯ − 2r11 q¯ ∗ , Cqq¯ q¯ − 2r02 q¯ ∗ , Cqqq¯             +2 Br11 r11 − q∗ , Br11r11 q − q¯ ∗ , Br11 r11 q¯ + Br02 r20 − q∗ , Br02 r20 q − q¯ ∗ , Br02r20 q¯ (29) +2 Bqr12 − q∗ , Bqr12 q − q¯ ∗ , Bqr12 q¯ + 2 Bqr q − q¯ ∗ , Bqr q¯ − q∗ , Bqr     ¯ 21  ¯ 21    ¯ 21  − q∗ , Cqr + Cqr ¯ ∗ , Cqr ¯ + Cqr02 r20 − q∗ , Cqr02 r20 q − q¯ ∗ , Cqr02 r20 q¯ ¯ 02 r20 q − q ¯ 02 r20 q             ¯ 02 r20 ∗ ¯ ∗ , Cqr ¯ +2 Cqr11 r11 − q∗ , Cqr11 r11 q − q¯ ∗ , Cqr11 r11 q¯ + 2 Cqr ¯ 11 r11 − q , Cqr ¯ 11 r11 q − q ¯ 11 r11 q             ∗ ¯ ∗ , Cqqr ¯ + Cqqr02 − q∗ , Cqqr02 q − q¯ ∗ , Cqqr02 q¯ + 4 Cqqr ¯ 11 − q , Cqqr ¯ 11 q − q ¯ 11 q             ∗ ¯ ∗ , Cq¯ qr ¯ + Dqqq¯ q¯ − q∗ , Dqqq¯ q¯ q − q¯ ∗ , Dqqq¯ q¯ q¯ , + Cq¯ qr ¯ 20 − q , Cq¯ qr ¯ 20 q − q ¯ 20 q 9

− [L (K0 ) − 4iω0 I] r40

=

− [L (K0 ) + 4iω0 I] r04

=

        −6r30 q∗ , Bqq − 12r20 q∗ , Bqr20 − 4r20 q∗ , Cqqq − 12r11 q¯ ∗ , Bqr20     −6r21 q¯ ∗ , Bqq − 4r11 q¯ ∗ , Cqqq             +3 Br20 r20 − q∗ , Br20 r20 q − q¯ ∗ , Br20 r20 q¯ + 4 Bqr30 − q∗ , Bqr30 q − q¯ ∗ , Bqr30 q¯       +3 Cqr20 r20 − q∗ , Cqr20 r20 q − q¯ ∗ , Cqr20 r20 q¯       ¯ ∗ , Cqr ¯ − q∗ , Cqr +3 Cqr ¯ 20 r20 q − q ¯ 20 r20 q            ¯ 20 r20  ∗ ∗ +6 Cqqr20 − q , Cqqr20 q − q¯ , Cqqr20 q¯ + Dqqqq − q∗ , Dqqqq q − q¯ ∗ , Dqqqq q¯ ,       ∗ ∗ −12r02 q¯ ∗ , Bqr ¯ 02 − 12r11 q , Bqr ¯ 02 − 6r12 q , Bq¯ q¯       −6r03 q¯ ∗ , Bq¯ q¯ − 4r02 q¯ ∗ , Cq¯ q¯ q¯ − 4r11 q∗ , Cq¯ q¯ q¯             ∗ +3 Br02 r02 − q∗ , Br02 r02 q − q¯ ∗ , Br02 r02 q¯ + 4 Bqr ¯ ∗ , Bqr ¯ ¯ 03 − q , Bqr ¯ 03 q − q ¯ 03 q       ∗ ∗ +3 Cqr02 r02 − q , Cqr02 r02 q − q¯ , Cqr02 r02 q¯       ¯ ∗ , Cqr ¯ − q∗ , Cqr +3 Cqr ¯ 02 r02 q − q ¯ 02 r02 q            ¯ 02 r02  ∗ ∗ q − q ¯ q ¯ + Dq¯ q¯ q¯ q¯ − q∗ , Dq¯ q¯ q¯ q¯ q − q¯ ∗ , Dq¯ q¯ q¯ q¯ q¯ , − q , C , C +6 Cq¯ qr ¯ 02 q¯ qr ¯ 02 q¯ qr ¯ 02

where I is the unit 3 × 3 matrix. Proof. To prove the Lemma, plug in the right hand side of (26) the Taylor expansion of H(z, z¯, r) given by (18). From (25) one obtains the derivatives of r = V(z, z) with respect to z and z¯; on the other hand, from (16) one obtains the dz ∂r d¯z derivatives of z and z¯ with respect to t. Next, replace these derivatives in the expression ∂r ∂z dt + ∂¯z dt . Finally, the result ∂r dz ∂r d¯z follows by equating coefficients of like power in z and z¯ in the equation ∂z dt + ∂¯z dt = L(K0 )r + H(z, z¯, r).

2.1. Calculation of the Lyapunov Coefficients In this subsection we obtain expressions for the first and second Lyapunov coefficients l1 and l2 , respectively. Equivalent definitions and algorithmic procedures to write the expressions for the Lyapunov coefficients for two dimensional ODE systems can be found in [1, 3], among others. These procedures apply also to n-dimensional ODE systems, if properly restricted to the center manifold. In [17, 29], the authors outlined a method for n-dimensional ODE systems requiring no explicit formal evaluation of the center manifold. It is worth observing that the method outlined in this paper pursues forward the algorithmic procedure presented in [17, 29] since we obtain expressions for the Lyapunov coefficients for a reaction diffusion systems subject to Newmann boundary conditions. Some of the results presented in this subsection are generalizations, for a infinite dimensional reaction diffusion systems, of the results presented in [17, 29]. To obtain expressions for the Lyapunov coefficients l1 and l2 we need to calculate some of the terms rk j , 2 ≤ j+k ≤ 5, given in Lemma 2. In [11] p. 1951, the authors presented the procedure to calculate r20 and r11 . To this end, they observed that in the equations given by (27) the terms ±2iω0 and 0 are not eigenvalues of L (K0 ), and the spectrum of L (K0 ) consists only of eigenvalues. Furthermore, the calculation of [2iω0 I − L (K0 )]−1 and [L (K0 )]−1 are restricted to the subspaces spanned by the eigen-modes 1 and cos( 2mx  ). In the following, we present the procedure to calculate r21 . To this purpose, we need to calculate the inner products involved in the equation for r21 given in (28). Noticing that mx mx (am , bm , cm )T , q¯ = cos (¯am , b¯ m , c¯ m )T and − [L (K0 ) − 2iω0 I] r20 = Bqq     + r20,1 where = r20,c cos 2mx  q = cos

we get r20

r20,c

⎛⎛ ⎜⎜⎜⎜⎜⎜ am 1 ⎜⎜ = [2iω0 I − L2m (Km )]−1 B ⎜⎜⎜⎜⎜⎜⎜⎜ bm ⎝⎝ 2 cm

⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ am ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎠ , ⎜⎜⎝ bm cm

⎛⎛ ⎞⎞ ⎜⎜⎜⎜ am ⎟⎟⎟⎟⎟⎟ 1 ⎟⎟⎟⎟⎟⎟⎟⎟ and r20,1 = [2iω0 I − L0 (Km )]−1 B ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜ bm ⎜⎝⎜⎝ ⎟⎠⎟⎠ 2 cm 10

⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ am ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎠ , ⎜⎜⎝ bm cm

⎞⎞ ⎟⎟⎟⎟⎟⎟ ⎟⎟⎟⎟⎟⎟⎟⎟ . ⎟⎠⎟⎠

thus, we calculate Bqr ¯ 20 Analogously Bqr11 where

⎛⎛ ⎜⎜⎜⎜⎜⎜ a¯ m ⎜⎜ = B ⎜⎜⎜⎜⎜⎜⎜⎜ b¯ m ⎝⎝ c¯ m

⎞ ⎞ ⎛⎛ ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜⎜⎜⎜ a¯ m mx 2mx ⎟⎟⎟ ⎟⎟⎟ ⎜⎜ cos + B ⎜⎜⎜⎜⎜⎜⎜⎜ b¯ m ⎟⎟⎠ , r20,c ⎟⎟⎠ cos ⎝⎝   c¯ m

⎞ ⎞ ⎟⎟⎟ ⎟⎟⎟ mx ⎟⎟⎟ ⎟ . ⎟⎟⎠ , r20,1 ⎟⎟⎟⎟⎠ cos 

⎛⎛ ⎜⎜⎜⎜⎜⎜ am ⎜⎜ = B ⎜⎜⎜⎜⎜⎜⎜⎜ bm ⎝⎝ cm

⎞ ⎞ ⎛⎛ ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜⎜⎜⎜ am mx 2mx ⎟⎟⎟ ⎟⎟⎟ ⎜⎜ cos + B ⎜⎜⎜⎜⎜⎜⎜⎜ bm ⎟⎟⎠ , r11,c ⎟⎟⎠ cos ⎝⎝   cm

⎞ ⎞ ⎟⎟⎟ ⎟⎟⎟ mx ⎟⎟⎟ ⎟ , ⎟⎟⎠ , r11,1 ⎟⎟⎟⎟⎠ cos 



r11,c

1 = − L2m (Km ) 2

Finally, we can calculate

−1

⎛⎛ ⎜⎜⎜⎜⎜⎜ am ⎜⎜ B ⎜⎜⎜⎜⎜⎜⎜⎜ bm ⎝⎝ cm

⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ a¯ m ⎟⎟⎟ ⎜⎜⎜ ¯ ⎟⎟⎠ , ⎜⎜⎝ bm c¯ m

⎞⎞ ⎛⎛  −1 ⎜⎜⎜⎜ am ⎟⎟⎟⎟⎟⎟ ⎜⎜⎜⎜ 1 ⎟⎟⎟⎟⎟⎟ ⎟⎟⎠⎟⎟⎠ and r20,1 = − L0 (Km ) B ⎜⎜⎜⎜⎝⎜⎜⎜⎜⎝ bm 2 cm

⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ a¯ m ⎟⎟⎟ ⎜⎜⎜ ¯ ⎟⎟⎠ , ⎜⎜⎝ bm c¯ m

⎞⎞ ⎟⎟⎟⎟⎟⎟ ⎟⎟⎟⎟⎟⎟ ⎟⎟⎠⎟⎟⎠ .

         ∗ ∗ ¯ ∗ , Bqr11 , q∗ , Bqqq¯ . q¯ ∗ , Bqr ¯ 20 , q , Bqr ¯ 20 , q , Bqr11 , q           ∗ ∗ ¯ ∗ , Bqr11 , q∗ , Bqqq¯ are long to present here, we refer the As the expressions for q¯ ∗ , Bqr ¯ 20 , q , Bqr ¯ 20 , q , Bqr11 , q readers to Appendix 2. To carry out the above calculations, it is worth observing that for m  0 we have π

cos2



0

mx 





dx =

π 2

,

π

cos



0

2mx 



cos2



mx 



dx =

π 4

,

π

cos4

0



mx 



dx =

3π 8 ;

on the other hand, for m = 0 we have π

cos2

0



mx 



dx =

π 0

cos



2mx 



cos2



mx 



dx =

π

cos4

0



mx 



dx = π.

From (28) we obtain a singular system for r21 − [L (K0 ) − iω0 I] r21

=

        −2r20 q∗ , Bqq¯ − r11 q∗ , Bqq − r02 q¯ ∗ , Bqq − 2r11 q¯ ∗ , Bqq¯             ¯ ∗ , Bqr ¯ + 2 Bqr11 − q∗ , Bqr11 q − q¯ ∗ , Bqr11 q¯ + Bqr − q∗ , Bqr ¯ 20 q − q ¯ 20 q      ¯ 20  + Cqqq¯ − q∗ , Cqqq¯ q − q¯ ∗ , Cqqq¯ q¯

(30)

which has a solution if and only if q∗ , a = 0, where         a = −2r20 q∗ , Bqq¯ − r11 q∗ , Bqq − r02 q¯ ∗ , Bqq − 2r11 q¯ ∗ , Bqq¯             ¯ ∗ , Bqr ¯ + 2 Bqr11 − q∗ , Bqr11 q − q¯ ∗ , Bqr11 q¯ − q∗ , Bqr + Bqr ¯ 20 q − q ¯ 20 q      ¯ 20  + Cqqq¯ − q∗ , Cqqq¯ q − q¯ ∗ , Cqqq¯ q¯ . We can find r21 by solving the nonsingular system [iω0 I − L (K0 )] r21 + sq

q∗ , r21  + 0s

= a = 0

(31)

with the condition q∗ , r21  = 0. For the sake of completeness, in the next proposition we prove that the homogeneous system associated with system (31) has a unique solution and if (r21 , s) is a solution of (31) with the condition q∗ , r21  = 0, then r21 is a solution of (30). Proposition 5. The non-homogeneous system (31) is non-singular. Furthermore, if (r21 , s) is a solution of (31) with the condition q∗ , a = 0, then r21 is a solution of (30). Reciprocally, if r21 is a solution of (30) then q∗ , a = 0. 11

Proof. Write X = X c ⊕ X su , where X c , X su are invariant spaces by L(K0 ). From the definition of X su , r21 ∈ X su if and only if q∗ , r21  = 0. Let (r21 , s) be a solution of the homogeneous equation obtained from (31), i. e., [iω0 I − L(K0 )] r21 + sq = 0 = 0.

q∗ , r21 

(32)

From the second equation of (32) it follows that r21 ∈ X su thus, (iω0 I − L (K0 ))r21 ∈ X su which implies that

q∗ , (iω0 I − L (K0 ))r21  = 0. Taking the inner product between q∗ and both sides of the first equation in (32) one obtains q∗ , (iω0 I − L (K0 ))r21  + s q∗ , q = 0, which implies s q∗ , q = 0 . Since q∗ , q = 1 it follows that s = 0. Consequently, [iω0 I − L (K0 )] r21 = 0 or, equivalently L (K0 ) r21 = iω0 r21 ; therefore, there exists α ∈ C such that r21 = αq. Additionally, from the second equation in (32) we have 0 = q∗ , r21  = q∗ , αq = α q∗ , q = α, i.e., (r21 , s) = (0, 0). Consequently, the non-homogeneous system (31) is non-singular, i.e., it has an unique solution. Let (r21 , s) solution of (31) and assumes that q∗ , a = 0. From the second equation of (31) it follows that r21 ∈ X su thus, (iω0 I − L (K0 ))r21 ∈ X su which implies q∗ , (iω(0 I − L (K0 ))r21  = 0. Taking the inner product of q∗ with both sides of the first equation in (31) one obtains q∗ , (iω0 I − L (K0 ))r21  + s q∗ , q = q∗ , a = 0 or, equivalently s q∗ , q = 0. Since q∗ , q = 1 it follows that s = 0. Substituting s = 0 into the first equation of (31) results [iω0 I − L (K0 )] r21 = a. Therefore r21 is a solution of (30). Reciprocally, let r21 solution of (30). Since q∗ , r21  = 0, i.e., r21 ∈ X su and the space X su is invariant by L(K0 ) then, a = [iω0 I − L (K0 )] r21 ∈ X su . Hence, q∗ , a = 0. Remark 4. The procedure described above can be adapted in connection with the determination of r12 , r32 and further coefficients in order to obtain the expressions for the third and fourth Lyapunov coefficients with few modifications. We defer our work in this direction to a subsequent exposition since this analysis is beyond the scope of the current paper. Once we get the coefficients of the Taylor expansion of the center manifold up to terms of order five , it follows from Appendix A of [6] that the system (19) restricted to the center manifold is given by the first equation of (19) with r = V(z, z). Therefore, the restriction of the reaction-diffusion system (19) to the center manifold up to terms of order five, can be written as  gi j i j 6 dz = iω0 z + q∗ , F0 (z, z¯, r) = iω0 z + (33) dt i! j! z z¯ + o(|z| ). 2≤i+ j≤5 Consequently, the dynamics of (19) can be determined by the dynamics of (33) where           g20 = q∗ , Bqq g11 = q∗ , Bqq¯ ; g02 = q∗ , Bq¯ q¯ ; g30 = 3 q∗ , Bqr20 + q∗ , Cqqq ;             ∗ ∗ ∗ ∗ g21 = 2 q∗ , Bqr11 + q∗ , Bqr ¯ 20 + q , Cqqq¯ ; g12 = q , Bqr02 + 2 q , Bqr ¯ 11 + q , Cqq¯ q¯ ;            ∗  ∗ ∗ ∗ ∗ g03 = 3 q∗ , Bqr ¯ 02 ; + q , Cq¯ q¯ q¯ ; g40 = 4 q , Bqr30 + 3 q , Br20 r20 + 6 q , Cqqr20 + q , Dqqqq ;             g31 = 3 q∗ , Bqr21 + q∗ , Bqr + 3 q∗ , Br20 r11 + q∗ , Dqqqq¯ + 3 q∗ , Cqqr11 + 3 q∗ , Cqqr ¯ 20 ;   ¯ 30    ∗   ∗  g22 = 2 q∗ , Bqr12 + 2 q∗ , Bqr ¯ 21 + 2 q , Br11 r11 + q , Br20 r02         ∗ + q∗ , Cq¯ qr + q∗ , Cqqr02 + 4 q∗ , Cqqr ¯ 20 + q , Dqqq¯ q¯ ;    ¯ 11          ∗ ∗ g13 = q∗ , Bqr03 + 3 q∗ , Bqr + 3 q∗ , Br11 r02 + 3 q∗ , Cqqr ¯ 02 + 3 q , Cq¯ qr ¯ 11 + q , Dqq¯ q¯ q¯ ;   ¯ 12         + 3 q∗ , Br20 r12 + q∗ , Br02 r30 + 6 q∗ , Br21 r11 g32 = 3 q∗ , Bqr22 + 2 q∗ , Bqr   ¯ 31           ∗ ∗ +3 q∗ , Cqr20 r02 + 3 q∗ , Cqqr12 + 6 q∗ , Cqr11 r11 + 6 q∗ , Cqr ¯ 20 r11 + q , Cq¯ qr ¯ 30 + 6 q , Cqqr ¯ 21         ∗ ∗ + q∗ , Dqqqr02 + 6 q∗ , Dqqqr ¯ 11 + 3 q , Dqq¯ qr ¯ 20 + q , Eqqqq¯ q¯ .

(34)

Remark 5. The quantities gi j , 2 ≤ i + j ≤ 5, given in (34) are obtained plugging equation (35) in (33) and then

12

equating coefficients of like power in z and z¯ in the equation (33) after substitution.            iω0 z + q∗ , F0 (z, z¯, r) = iω0 z + 12 q∗ , Bqq z2 + q∗ , Bqq¯ z¯z + 12 q∗ , Bq¯ q¯ z¯2 + 12 q∗ , Bqr20 + 16 q∗ , Cqqq z3       1 ∗ 2 + q∗ , Bqr11 + 12 q∗ , Bqr ¯ 20 + 2 q , Cqqq¯ z z¯             1 1 1 ∗ ∗ ∗ ∗ 3 + 2 q , Bqr02 + q , Bqr z2 + 12 q∗ , Bqr ¯ 11 + 2 q , Cqq¯ q¯ z¯ ¯ 02 + 6 q , Cq¯ q¯ q¯ z¯          1 q∗ , Dqqqq z4 + 16 q∗ , Bqr30 + 18 q∗ , Br20 r20 + 14 q∗ , Cqqr20 + 24        3 1 1 1  ∗ ∗ ∗ + 2 q , Bqr21 + 6 q , Bqr ¯ 30 + 2 q , Br20 r11 + z z¯        (35) 3 + 16 q∗ , Dqqqq¯ + 12 q∗ , Cqqr11 + 12 q∗ , Cqqr ¯ 20 z z¯            1 1 1 ∗ ∗ ∗ 2 2 + 12 q∗ , Bqr12 + 12 q∗ , Bqr ¯ 21 + 2 q , Br11 r11 + 4 q , Br20 r02 + 4 q , Cqqr02 z z¯       1 1 ∗ ∗ ∗ 2 2 + q , Cqqr ¯ 11 + 4 q , Cq¯ qr ¯ 20 + 4 q , Dqqq¯ q¯ z z¯       3 1 1 1  ∗ ∗ ∗ + 6 q , Bqr03 + 2 q , Bqr z ¯ 12 + 2 q , Br11 r02 z¯        1 1 ∗ ∗ z3 + 12 q∗ , Cqqr ¯ 02 + 2 q , Cq¯ qr ¯ 11 + 6 q , Dqq¯ q¯ q¯ z¯         ∗ 1 ∗ 1 4 + q∗ , 18 Br02 r02 + q∗ , 16 Bqr ¯ 03 + q , 4 Cq¯ qr ¯ 02 + q , 24 Dq¯ q¯ q¯ q¯ z¯      1  ∗  1 ∗ 1 1  ∗ ∗ + 24 q , Bqr40 + 24 q , Br30 r20 + 24 q , Br20 r30 + 8 q , Cqr20 r20 z5        1 1 1 q∗ , Cqqr30 + 12 q∗ , Dqqqr20 + 120 q∗ , Eqqqqq z5 + 12        1 ∗  4 1 1 ∗ q∗ , Bqr + 16 q∗ , Bqr31 + 24 ¯ 40 + 8 q , Br20 r21 + 8 q , Br21 r20 z z¯        1  ∗  1 4 q∗ , Br11 r30 + 12 q , Br30 r11 + 12 q∗ , Cqr20 r11 + 18 q∗ , Cqr + 12 ¯ 20 r20 z z¯            1 1 1 1 1 ∗ ∗ ∗ 4 + 6 q∗ , Dqqqr11 + 4 q∗ , Dqqqr ¯ 20 + 4 q , Cqqr21 + 6 q , Cqqr ¯ 30 + 24 q , Eqqqqq¯ z z¯       1 ∗  1  ∗  3 2 1 1 1  ∗ ∗ ∗ + 4 q , Bqr22 + 6 q , Bqr ¯ 31 + 8 q , Br20 r12 + 8 q , Br12 r20 + 24 q , Br02 r30 z z¯            1 + 24 q∗ , Br30 r02 + 14 q∗ , Br21 r11 + 14 q∗ , Br11 r21 + 14 q∗ , Cqr20 r02 + 14 q∗ , Cqqr12 z3 z¯2          1 1 ∗ ∗ 3 2 + 12 q∗ , Cqr11 r11 + 12 q∗ , Cqr ¯ 20 r11 + 12 q , Cq¯ qr ¯ 30 + 2 q , Cqqr ¯ 21 z z¯          1 1 1 1 ∗ ∗ ∗ ∗ 3 2 + q , Dqq¯ qr + 12 q , Dqqqr02 + 2 q , Dqqqr ¯ 20 + 12 q , Eqqqq¯ q¯ z z¯    ¯ 11  4        2 3 1 1 1 ∗ ∗ ∗ + 16 q∗ , Bqr13 + 14 q∗ , Bqr ¯ 22 + 24 q , Br03 r20 + 24 q , Br20 r03 + 4 q , Br11 r12 z z¯            2 3 + 14 q∗ , Br12 r11 + 18 q∗ , Br21 r02 + 18 q∗ , Br02 r21 + 12 q∗ , Cqr11 r02 + 14 q∗ , Cqr ¯ 20 r02 z z¯          1 1 1 1 ∗ ∗ ∗ ∗ 2 3 + q , Cqqr03 + 4 q , Cq¯ qr z z¯ + 2 q , Cqr ¯ 21 + 2 q , Cqqr     ¯ 12    ¯ 11 r11  12  1 1 1 1 ∗ ∗ ∗ ∗ 2 3 + + 4 q , Dqqqr ¯ 02 + 2 q , Dqq¯ qr ¯ 11 + 12 q , Dq¯ q¯ qr ¯ 20 + 12 q , Eqqq¯ q¯ q¯ z z¯       1  ∗  1 ∗  4 1 1 1  ∗ ∗ ∗ z + 24 q , Bqr04 + 6 q , Bqr ¯ 13 + 12 q , Br03 r11 + 12 q , Br11 r03 + 8 q , Br12 r02 z¯          1 ∗ 4 + q z¯ z + + 18 q∗ , Br02 r12 + 18 q∗ , Cqr02 r02 + 12 q∗ , Cqr , C ¯ 11 r02 qqr ¯ 03      6    1 1 1 ∗ ∗ ∗ z4 + + 14 q∗ , Cq¯ qr ¯ 12 + 4 q , Dqq¯ qr ¯ 02 + 6 q , Dq¯ q¯ qr ¯ 11 + 24 q , Eqq¯ q¯ q¯ q¯ z¯       1  ∗  1 ∗ 1 1  ∗ ∗ 5 z¯ + 24 q , Bqr ¯ 04 + 24 q , Br02 r03 + 24 q , Br03 r02 + 8 q , Cqr        ¯ 02 r02 6 1 1 1 ∗ ∗ ∗ 5 + 12 q , Cq¯ qr ¯ 03 + 12 q , Dq¯ q¯ qr ¯ 02 + 120 q , Eq¯ q¯ q¯ q¯ q¯ z¯ + o |z| + o(|z| |r|). The Poincar´e normal form of (10) (for K in a neighborhood of K0 ) can be written in the form dz dt

=

(α + iω) z + z

M  j=1

c j (K) (z¯z) j ,

where α + iω = α(K) + iω(K), z is a complex variable, M ≥ 1 and c j (K) are complex-valued coefficients (see [6] pp. 28 and [11] pp. 1952). The quantity c1 (K0 ) is given by  i (36) c1 (K0 ) = 2ω0 g20 g11 − 2|g11 |2 − 13 |g02 |2 + g221 . (see [6] pp. 47; see also [11] pp. 1952 and [29] pp. 98). 13

Definition 1. The first Lyapunov coefficient is defined by  l1 (K0 ) =

Re(c1 (K0 )) ω0

=

1 Re 2ω20

ig20 g11 + ω0 g21 .

At the bifurcation point K0 where α(K0 ) = 0 and l1 (K0 ) = 0, the following definition gives a rather compact expression for l2 (K0 ) (see for example [29] pp. 310; see also [6] pp. 49). Definition 2. The second Lyapunov coefficient is defined by l2 (K0 ) =

Re(c2 (K0 )) ω0

 =

1 12ω40

  ω30 Reg32 + ω20 Im g20 g¯ 31 − g11 (4g31 + 3¯g22 ) − 13 g02 (g40 + g¯ 13 ) − g30 g12

+

     ω0 Re g20 g¯ 11 (3g12 − g¯ 30 ) + g02 (¯g12 − 13 g30 ) + 13 g¯ 02 g03 + g11 g¯ 02 ( 53 g¯ 30 + 3g12 )

+

1 ¯ 03 3 g02 g

+

!  Im(g20 g11 ) 3Re(g20 g11 ) − 2|g02 |2 .

     − 4g11 g30 + 3Im(g20 g11 )Im(g21 ) + Im g11 g¯ 02 g¯ 220 − 3¯g20 g11 − 4g211

Definition 3. i) A Hopf point (K0 , u0 , v0 , w0 ) is an equilibrium point of system (4) where L(K0 ) has a pair of purely imaginary eigenvalues Γ(K0 ) = ±iω0 , continuously differentiable in a, such that α(K0 ) = 0, ω(K0 ) = ω0 > 0, and admits no other critical eigenvalues i.e., located on the imaginary axis. i) A Hopf point is called transversal if the parameter dependent complex eigenvalues cross the imaginary axis with non-zero derivative. A Hopf point of codimension 1 is a Hopf point where l1  0. ii) A Hopf point of codimension 2 is a Hopf point where l1 vanishes. It is called transversal if α = 0 and l1 = 0 have transversal intersections, where α = α(K) is the real part of the critical eigenvalues. Remark 6. In a neighborhood of a codimension 1 transversal Hopf point, H1 point for succinctness, with l1  0 the dynamic behavior of the system (4), reduced to the family of parameter dependent continuations of the center manifold, is orbitally topologically equivalent to the following complex normal form z˙ = (η + i) z + l1 z|z|2 , where z ∈ C, α and l1 are real functions having derivatives of arbitrary high order, which are continuations of 0, ω0 and the first Lyapunov coefficient at the H1 point. Remark 7. In a neighborhood of a codimension 2 transversal Hopf point, H2 point for simplicity, with l2  0 the dynamic behavior of the system (4), reduced to the family of parameter dependent continuations of the center manifold, is orbitally topologically equivalent to the following complex normal form z˙ = (η1 + i) z + η2 z|z|2 + l2 z|z|4 , where z ∈ C, and η1 and η2 are unfolding parameters. We summarize the previous discussion in the following theorem. 14

Theorem 3. Assume that the Hopf hypothesis is satisfied. Then system (4) undergoes a Hopf bifurcation (see [6]). Furthermore 1 l1 (K0 ) < 0 (resp. > 0). If in addition all other i) The Hopf Bifurcation is supercritical (resp. subcritical) if α (K 0) eigenvalues of L(a0 ) have negative real parts, then the bifurcating periodic solutions are stable (resp. unstable) if l1 (K0 ) < 0 ( resp. l1 (K0 ) > 0).   ii) Suppose that for K = K0 , one has η1 (K0 ) = 0, l1 (K0 ) = 0 and the map a → η1 (K), l1 (K) is regular at K = K0 . Then, by the introduction of a complex variable, applying smooth invertible coordinate transformations that depend smoothly on the parameters, and performing smooth parameter and time changes, the system (4) reduced to the family of parameter-dependent continuations of the center manifold, is orbitally topologically equivalent to

z˙ = (η1 + i) z + η2 z|z|2 + l2 z|z|4 , where z ∈ C, η1 and η2 are unfolding parameters. The bifurcation diagrams for l2 < 0 can be found in [29], p. 313 (see also [10]). As an application of the method outlined in this section, in the next section we show that under certain hypothesis the system (1)-(2) undergoes a codimension 1 Hopf bifurcation. 3. Stability and bifurcation analysis for the system (1)-(2) In this section we discuss the its homogeneous equilibrium points (or homogeneous solutions) for the system (1)-(2) and investigate their asymptotic behavior. Next, we consider the related codimension 1 homogeneous and non-homogeneous Hopf bifurcation with the spatial domain Ω = (0, π),  ∈ R+ . It is important to highlight that in [8] the authors studied the corresponding Hopf bifurcation for the ODE version of system (1)-(2). However, to the best of our knowledge the direction of Hopf bifurcation and stability of the bifurcating periodic solutions were not studied in the current literature. Finally, we include a general algorithm to determine the direction and stability of the spatially non-homogeneous bifurcating periodic solutions for the system (1)-(2). 3.1. Equilibria and Stability Analysis In the following we discuss the homogeneous solutions of the system (1)-(2), that is the solutions of the system γ1 [1 + K − u − Kw]u, γ2 [u − v]v, γ3 [v − w]w,

= 0 = 0 = 0

(37)

Aside from the obvious homogeneous solutions (u0 , v0 , w0 ) = (0, 0, 0); (u1 , v1 , w1 ) = (1 + K, 0, 0) and (u2 , v2 , w2 ) = (1 + K, 1 + K, 0), the system (37) has the biologically interesting homogeneous solution (u3 , v3 , w3 ) = (1, 1, 1). To study the qualitative behavior of these homogeneous solutions, observe that the linearization of system (1)-(2) at a point (u∗ , v∗ , w∗ ) is given by ⎧ ∂u(x, t) ⎪ ⎪ = d1 Δu − γ1 (Kw∗ + 2u∗ − K − 1)u − γ1 Ku∗ w ⎪ ⎪ ∂t ⎪ ⎪ ⎨ ∂v(x, t) (38) = d2 Δv + γ2 v∗ u + γ2 (u∗ − 2v∗ )v ⎪ ⎪ ∂t ⎪ ⎪ ⎪ ⎪ ∂w(x, t) = d Δw + γ w∗ v + γ (v∗ − 2w∗ )w ⎩ 2 3 3 ∂t

with boundary and initial conditions given by (2). Let φ = (u, v, w), D = diag(d1 , d2 , d3 ) and J be the matrix ⎛ ∗ ∗ ⎜⎜⎜ −γ1 (Kw + 2u − K − 1) ⎜⎜⎜ ⎜⎜ γ2 v∗ J(K, u∗ , v∗ , w∗ ) = ⎜⎜⎜⎜⎜ ⎜⎜⎜ ⎜⎝ 0 15

0

−γ1 Ku∗

γ2 (u∗ − 2v∗ )

0

γ3 w∗

γ2 (v∗ − 2w∗ )

0

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ . ⎟⎟⎟ ⎟⎟⎠

Then (38) can be written as

∂φ = DΔφ + Jφ ∂t so, the solution of (38), (2) with initial condition φ(·, 0) = φ0 is given by φ(x, t) =

where φ0 , ψm  =

 Ω

∞ 

e(J−μm D)t φ0 , ψm ψm (x),

n=0

φ0 (x)ψm (x) dx.

It follows from the linearization principle that a homogeneous solution of (1)-(2) is locally asymptotically stable if the eigenvalues of all matrix J − μm D have negative real part; if there exists a m such that J − μm D has an eigenvalue with positive real part then the solution is unstable. Remark 8. We observe that if (u∗ , v∗ , w∗ ) is an unstable equilibrium for the flow of the corresponding ODE system of system (1)-(2) then (u∗ , v∗ , w∗ ) is also an unstable homogeneous equilibrium for the flow of (1)-(2), since the subspace of the function independent of x is invariant for the flow of (1)-(2). Proposition 6. The trivial homogeneous equilibria (u0 , v0 , w0 ) = (0, 0, 0), (u1 , v1 , w1 ) = (1 + K, 0, 0) and (u2 , v2 , w2 ) = (1 + K, 1 + K, 0) of system (1)-(2) are unstable for any strictly positive parameters γ1 , γ2 , γ3 , d1 , d2 , d2 and K. Proof. The instability of the equilibria (u0 , v0 , v0 ), (u1 , v1 , w1 ) and (u2 , v2 , w2 ) follow from [8] taking μ0 = 0 . In the following we will study the local stability of E = (u3 , v3 , w3 ) = (1, 1, 1). The linearization evaluated at E is given by ⎞ ⎛ 0 −γ1 K ⎟⎟⎟ ⎜⎜⎜ −(γ1 + μm d1 ) ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ . ⎜ −(γ2 + μm d2 ) 0 γ2 Lm (K) = ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎠ ⎜⎝ −(γ3 + μm d3 ) 0 γ3 (γ2 + γ3 )(γ1 + γ3 )(γ1 + γ2 ) , then the homogeneous equilibrium E = (1, 1, 1) of system γ1 γ2 γ3 (1)-(2) is locally asymptotically stable. Otherwise, it is unstable.

Proposition 7. If K < K0 =

Proof. The characteristic polynomial associated to Lm (K) is given by    p(λ) = λ3 + d1 μm + d2 μm + d3 μm + γ1 + γ2 + γ3 λ2 + d1 d2 μ2m + d1 d3 μ2m + d2 d3 μ2m + d1 γ2 μm + d1 γ3 μm + d2 γ1 μm  +d2 γ3 μm + d3 γ1 μm + d3 γ2 μm + γ1 γ2 + γ1 γ3 + γ2 γ3 λ + d1 d2 d3 μ3m + d1 d2 γ3 μ2m + d1 d3 γ2 μ2m + d2 d3 γ1 μ2m +γ1 γ2 γ3 K + d1 γ2 γ3 μm + d2 γ1 γ3 μm + d3 γ1 γ2 μm + γ1 γ2 γ3 . By the Routh-Hurwitz criteria the polynomial (39) is stable if   d1 μm + d2 μm + d3 μm + γ1 + γ2 + γ3 d1 d2 μ2m + d1 d3 μ2m + d2 d3 μ2m + d1 γ2 μm + d1 γ3 μm + d2 γ1 μm   +d2 γ3 μm + d3 γ1 μm + d3 γ2 μm + γ1 γ2 + γ1 γ3 + γ2 γ3 − + d1 d2 d3 μ3m + d1 d2 γ3 μ2m + d1 d3 γ2 μ2m + d2 d3 γ1 μ2m

(40)

 +γ1 γ2 γ3 K + d1 γ2 γ3 μm + d2 γ1 γ3 μm + d3 γ1 γ2 μm + γ1 γ2 γ3 > 0, (d2 μm + d3 μm + γ2 + γ3 )(d1 μm + d3 μm + γ1 + γ3 )(d1 μm + d2 μm + γ1 + γ2 ) . Therefore, if γ1 γ2 γ3 the characteristic polynomial (39) is stable for m = 0, that is, if K < K0 , then it will be stable for all m = 1, 2, .... Hence, from Theorem 1.7 given in [29] p. 36, the homogeneous solution E of (1)-(2) is locally asymptotically stable. On the other hand, if K < K0 , a straightforward calculation shows that for m = 0 p(λ) has a positive eigenvalue and therefore E is unstable. which is equivalent to K <

16

(39)

We obtained in Propositions 7 the local stability for the homogeneous equilibrium E of the system (1)-(2). Accordingly, in the next subsection we concentrate on the investigation of codimension 1 Hopf bifurcation in the system (1)-(2).

3.2. Spatially homogeneous and non-homogeneous Hopf bifurcation for the constant coexistence steady state E In this subsection we confine ourselves to the study of spatially homogeneous and non-homogeneous codimension 1 Hopf bifurcation in the system (1)-(2). Furthermore, we determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. To this purpose, we consider K as the bifurcation parameter to show (γ2 + γ3 )(γ1 + γ3 )(γ1 + γ2 ) that at K = K0 = the system (1)-(2) undergoes a codimension 1 Hopf bifurcation and γ1 γ2 γ3 the bifurcating periodic orbits are either spatially homogeneous, which coincides with the periodic solution of the corresponding ODE system, or spatially non-homogeneous. Furthermore, we show that there exist infinity Hopf bifurcation points where spatially non-homogeneous periodic orbits bifurcate. To follow our discussion according to   the framework of Section 2, we translate the critical point E = 1, 1, 1 of system (1)-(2) to the origin by the change of coordinates uˆ vˆ wˆ

= u − 1, = v − 1, = w − 1,

and still let u, v and w denote uˆ , vˆ and wˆ respectively. Thus, the system (1)-(2) take the form ut − d1 u xx = −γ1 (u + Kw)(u + 1), vt − d2 v xx = γ2 (u − v)(v + 1), wt − d3 w xx = γ3 (v − w)(w + 1),

x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0, x ∈ (0, π) , t > 0.

(41)

satisfying the conditions u x (0, t) = v x (0, t) = w x (0, t) = 0 u x (π, t) = v x (π, t) = w x (π, t) = 0 u (x, 0) = u0 (x) − 1, v (x, 0) = v0 (x) − 1, w (x, 0) = w0 (x) − 1

t > 0, t > 0, x ∈ (0, π) .

(42)

Setting f = f (K, u, v, w) = −γ1 (u + Kw)(u + 1), g = g(K, u, v, w) = γ2 (u − v)(v + 1) and h = h(K, u, v, w) = γ3 (v − w)(w + 1), as in Section 2, we consider the linearization of system (41)-(42) around the origin E0 = (0, 0, 0) ⎛ ∂2 ⎜⎜⎜ d1 ∂x 2 + A(K) ⎜⎜⎜ L (K) := ⎜⎜⎜ D(K) ⎝ G(K)

B(K) ∂ d2 ∂x 2 + E(K) H(K) 2

C(K) F(K) ∂2 d3 ∂x 2 + I(K)

⎛ ⎞ 2 ⎜⎜⎜ A(K) − d1m2 ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎟⎟⎟ and Lm (K) := ⎜⎜⎜ D(K) ⎜⎝ ⎠ G(K)

B(K) 2 E(K) − d2m2 H(K)

where C(K) = fw (K, 0, 0, 0) = −γ1 K B(K) = fv (K, 0, 0, 0) = 0, A(K) = fu (K, 0, 0, 0) = −γ1 , E(K) = gv (K, 0, 0, 0) = −γ2 , F(K) = gw (K, 0, 0, 0) = 0 D(K) = gu (K, 0, 0, 0) = γ2 , H(K) = hv (K, 0, 0, 0) = γ3 , G(K) = hu (K, 0, 0, 0) = 0, I(K) = hw (K, 0, 0, 0) = −γ3 . The characteristic equation of Lm (K) is given by λ3 + Mm (K)λ2 + Nm (K)λ + Pm (K) = 0, m = 0, 1, 2, · · · , 17

C(K) F(K) 2 I(K) − d3m2

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ (43) ⎟⎠

where Mm (K)

Nm (K) Pm (K)

(d1 + d2 + d3 )m2 2   (γ2 + γ3 )d1 + (γ1 + γ3 )d2 + (γ1 + γ2 )d3 m2 (d1 d2 + d1 d3 + d2 d3 )m4 = γ1 γ2 + γ1 γ3 + γ2 γ3 + + 2 4 = γ1 + γ2 + γ3 +

= (1 + K)γ1 γ2 γ3 +

(γ2 γ3 d1 + γ1 γ3 d2 + γ1 γ2 d3 )m2 (d1 d2 γ3 + d1 d3 γ2 + d2 d3 γ1 )m4 d1 d2 d3 m6 + + 2 4 6

(44)

and δm (K), λm (K) = αm (K) ± iωm (K) are eigenvalues of Lm , with δm (K) ∈ R and  δm (K)

=

" 36Mm Nm − 108Pm − 8Mm3 + 12 12Mm3 Pm − 3Mm2 Nm2 − 54Mm Nm Pm + 12Nm3 + 81P2m 2Nm −

−

"

1 3

2Mm2 3

36Mm Nm − 108Pm − 8Mm3 + 12 12Mm3 Pm − 3Mm2 Nm2 − 54Mm Nm Pm + 12Nm3 + 81P2m

αm (K)

=

 " 1 − 12 36Mm Nm − 108Pm − 8Mm3 + 12 12Mm3 Pm − 3Mm2 Nm2 − 54Mm Nm Pm + 12Nm3 + 81P2m Nm −

+

Mm2 3

" 36Mm Nm − 108Pm − 8Mm3 + 12 12Mm3 Pm − 3Mm2 Nm2 − 54Mm Nm Pm + 12Nm3 + 81P2m

ωm (K)

=

1 3

1 3

 

√

3 1 2 6

−

1 3

−

" 36Mm Nm − 108Pm − 8Mm3 + 12 12Mm3 Pm − 3Mm2 Nm2 − 54Mm Nm Pm + 12Nm3 + 81P2m 2Nm −

− 36Mm Nm − 108Pm −

8Mm3

+ 12

"

12Mm3 Pm

2Mm2 3

−

3Mm2 Nm2

− 54Mm Nm Pm +

12Nm3

+

1 3

1 Mm 3

1 Mm 3

(45)

1 3

 1 − Mm . 3

81P2m

We shall identify the Hopf bifurcation values K that satisfies the Hopf hypothesis that is, for which there exist a m ∈ N ∪ {0} such that (46) αm (Km ) = 0, ωm (Km ) > 0, and α j (Km )  0, for j  m; and for the unique pair of complex eigenvalues near the imaginary axis λ(K) = α(K) ± iω(K), we shall verify that α (Km )  0. Note that if K < K0 , we conclude from Proposition 7 that E (respectively E0 ) is locally asymptotically stable for the system (41)-(42). This suggest that any potential bifurcation point K must be in the interval [K0 , ∞). (γ2 + γ3 )(γ1 + γ3 )(γ1 + γ2 ) , is the smallest bifurcation value at which E0 The next proposition shows that K0 = γ1 γ2 γ3 (respectively E) is a spatially homogeneous transversal Hopf point for system (41)-(42) (respectively for system (1)-(2)). Furthermore, the next proposition shows that the system (1)-(2) has infinity points where undergoes nonhomogeneous Hopf bifurcation. Most of the calculations in Proposition 8 are too long and tedious therefore, its is conducted in MAPLE 13 to obtain the required conclusions. Proposition 8. At K = Km the system (41)-(42) satisfies the Hopf hypothesis. Consequently, E0 is a transversal Hopf point for the system (41)-(42) which undergoes a Hopf bifurcation at K = Km . for m = 0, 1, 2, 3... and K ∈ [K0 , ∞). 18

Proof. From (44) we obtain δ0 (K0 ) = −(γ1 + γ2 + γ3 ) < 0, λ0 (K0 ) = α0 (K0 ) ± iω0 (K0 ) where α0 (K0 ) = 0 and √ ω0 (K0 ) = γ1 γ2 + γ1 γ3 + γ2 γ3 > 0. Therefore, L0 (K0 ) has a negative eigenvalue and a pair of purely imaginary eigenvalues. Furthermore, we observe that for each m = 1, 2, 3, ... , αm (K) = 0 if and only if  (d + d )(γ + γ ) + (γ + γ )(d + d )(γ + γ ) 2 3 1 3 2 3 1 3 1 2 (γ2 + γ3 )(γ1 + γ3 )(γ1 + γ2 ) Km = + 2 γ1 γ2 γ3 γ1 γ2 γ3    (γ2 + γ3 )(γ1 + γ3 )(d1 + d2 ) 2 (d2 + d3 )(d1 + d3 )(γ1 + γ2 ) + + m + γ1 γ2 γ3 2 γ1 γ2 γ3 4  +

 (d2 + d3 )(γ1 + γ3 ) + (γ2 + γ3 )(d1 + d3 ) (d1 + d2 )  γ1 γ2 γ3

m4 +

4

(47)

(d2 + d3 )(d1 + d3 )(d1 + d2 ) 6 m . γ1 γ2 γ3 6

Hence, αm (K0 )  0 since K0 < Km , m = 1, 2, 3, .... Consequently, all other eigenvalues of Lm (K0 ) has negative real part for m = 1, 2, 3, .... On the other hand, for K = Km we have αm (Km ) = 0,

ωm (Km ) = and

√4

 γ1 γ2 +4 γ1 γ3 +4 γ2 γ3 +2 m2 d1 γ2 +2 m2 d1 γ3 +2 m2 d2 γ1 +2 m2 d2 γ3 +2 m2 d3 γ1 +2 m2 d3 γ2 +m4 d1 d2 +m4 d1 d3 +m4 d2 d3 2

>0

2 γ1 + 2 γ2 + 2 γ3 + m2 d1 + m2 d2 + m2 d3 < 0. 2   Furthermore, αm K j  0 if j  m. Consequently, all other eigenvalues of Lm (K j ) has no null real part for j  m. d Finally, with the help of MAPLE 13 we can check that dK αm (Km ) = 0 if, and only if γ1 = 0 or γ2 = 0 or d γ3 = 0. Hence, dK αm (Km )  0 for all m = 0, 1, 2, 3... where Km ∈ [K0 , ∞), since all parameters in model (1)-(2) are positive.

δm (Km ) = −

Remark 9. The result in Proposition 8 for the case m = 0 corresponds to the Hopf bifurcation of spatially homogeneous periodic solution which is known from previous studies [8]. Apparently K = K0 is also the unique value of K for the Hopf bifurcation of spatially homogeneous periodic solution for any  > 0. In the following we concentrate on the study of the spatially non-homogeneous codimension 1 Hopf bifurcation in system (41)-(42) for m ≥ 1 that is, we shall identify the values of the bifurcation parameter K which satisfies the Hopf hypothesis. To this purpose, as observed previously, any potential bifurcation point K must be in the interval [K0 , ∞). For d1 , d2 , d3 , γ1 , γ2 , γ3 and  fixed, we define Km , m = 0, 1, 2, 3, ..., as the roots of the equation α(K) = 0 where Km is given by (47). Obviously, these roots satisfies K0 < K1 < K2 < · · · < Km < · · ·. Remark 10. For the bifurcation values K = K j with j = 1, 2, . . . , m, . . ., we have from (47) that ⎧ ⎪ α (K ) = ⎪ ⎪ ⎨ m j α ⎪ m (K j )  ⎪ ⎪ ⎩ α (K )  m j

0 if m = j 0 if m < j 0 if m > j,

since K j < Km if j < m or vice-versa. We collect the above results and apply Theorem 3 to obtain the following Theorem: Theorem 4. Assume that d1 , d2 , d3 , γ1 , γ2 , γ3 and  are positive constants. Then there exist infinity points K j (m), j = 0, 1, 2, . . . , m, . . . satisfying K0 < K1 < K2 < · · · < Km < · · · such that the system (41)-(42) undergoes a bifurcation at these points. 19

i) The bifurcating periodic solutions from K = K0 are spatially homogeneous, which coincides with the periodic solution of the corresponding ODE system. ii) The bifurcating periodic solutions from K = K j , j = 1, 2, . . . , m, . . . are spatially non-homogeneous. In the following, we consider the bifurcation direction and stability of the bifurcating spatially homogeneous periodic solutions. Theorem 5. Consider the one-parameter family of differential equations (41)-(42). At (γ2 + γ3 )(γ1 + γ3 )(γ1 + γ2 ) the system (41)-(42) undergoes a codimension 1 Hopf bifurcation around the K = K0 = γ1 γ2 γ3 l1 (K0 )  < 0 (respectively, l1 (K0 )  > 0), the one-parameter family of differential equations (41)equilibrium E0 . If d dK

α0 (K0 )

d dK

α0 (K0 )

(42) undergoes a supercritical (respectively, subcritical) codimension 1 Hopf bifurcation at K = K0 . Furthermore, the bifurcating spatially homogeneous periodic solutions from K = K0 are locally asymptotically stable (respectively, unstable) if l1 (K0 ) < 0 (respectively, l1 (K0 ) > 0), which coincides with the periodic solution of the corresponding ODE system. Proof. According with Theorem 3, to determine the stability and bifurcation direction of the bifurcating periodic solution, we shall calculate Re(c1 (K0 )) where c1 (K0 ) is given by (36). Following the framework of Section 2, note that   ⎛ ⎞ 2 3 3 ⎜⎜⎜ −γ1 0 − γ1γ+γ + γ1γ+γ + γ2γ+γ + 2 γ1 ⎟⎟⎟ 3 2 1 ⎜ ⎟⎟⎟⎟ . L0 (K0 ) = ⎜⎜⎜⎜⎜ γ2 −γ2 0 ⎟⎟⎠ ⎝ −γ3 0 γ3 √ The eigenvalues associated with L0 (K0 ) are ±iω(K0 ) and δ0 (K0 ) = −(γ1 +γ2 +γ3 ) where ω0 (K0 ) = γ1 γ2 + γ1 γ3 + γ2 γ3 > 0; the eigenvectors of L0 (K0 ) associated to the complex and real eigenvalues are given, respectively, by  ⎞  ⎞ ⎛ (γ2 +γ3 )  ⎛ (γ2 +γ3 )  √ √ ⎜⎜⎜ − γ γ γ1 − i γ1 γ2 + γ1 γ3 + γ2 γ3 ⎟⎟⎟ ⎛⎜⎜ a0 ⎞⎟⎟ ⎜⎜⎜ − γ γ γ1 + i γ1 γ2 + γ1 γ3 + γ2 γ3 ⎟⎟⎟ ⎛⎜⎜ a¯ 0 ⎞⎟⎟ 2 3 2 3    ⎟ ⎟ ⎜⎜⎜ ⎟ ⎟⎟⎟ ⎜⎜⎜ ⎜ ⎜ √ √ ⎟⎟⎟⎟ = ⎜⎜⎜⎜ b0 ⎟⎟⎟⎟ , q¯ = ⎜⎜⎜⎜ ⎟⎟⎟ = ⎜⎜⎜ b¯ 0 ⎟⎟⎟⎟⎟ q = ⎜⎜⎜ 1 + γ13 i γ1 γ2 + γ1 γ3 + γ2 γ3 1 − γ13 i γ1 γ2 + γ1 γ3 + γ2 γ3 ⎠⎟ ⎠ ⎜⎝ ⎟⎟⎠ ⎝⎜ ⎟⎠ ⎝ ⎜⎜⎝ c0 c¯ 0 1 1 and

 (γ

1 +γ2 )(γ1 +γ3 )

γ2 γ 3

,−

T (γ1 +γ3 ) γ3 , 1 .

On the other hand we have ⎛ ⎜⎜⎜ ⎜ L0∗ (K0 ) = ⎜⎜⎜⎜⎜  ⎝ − γ1 +γ2 + γ3

−γ1 0 γ1 +γ3 γ2 +

γ2 +γ3 γ1

 + 2 γ1

γ2 −γ2 0

0 γ3 −γ3

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ , ⎠

whose eigenvalues are the same as the eigenvalues of L0 (K0 ) and the eigenvector associated to the complex and real eigenvalues given, respectively, by ⎛ γ2 γ3 (γ3 +i √γ1 γ2 +γ1 γ3 +γ2 γ3 ) ⎜⎜⎜ − (γ1 +γ√2 )(γ1 +γ3 )(γ2 +γ3 ) ⎜⎜ ∗ q = ⎜⎜⎜⎜⎜ γ3 (γ2 −i γ1 γ2 +γ1 γ3 +γ2 γ3 ) (γ1 +γ2 )(γ2 +γ3 ) ⎜⎝ 1

⎛ γ2 γ3 (γ3 −i √γ1 γ2 +γ1 γ3 +γ2 γ3 ) ⎞ ⎛ ⎜⎜⎜ − ⎟⎟⎟ ⎜ a∗ ⎞⎟ (γ1 +γ√2 )(γ1 +γ3 )(γ2 +γ3 ) ⎜ ⎟⎟⎟ ⎜⎜⎜ 0 ⎟⎟⎟ ⎟⎟⎟ = ⎜⎜⎜ b∗ ⎟⎟⎟ , q¯ ∗ = ⎜⎜⎜⎜⎜ γ3 (γ2 +i γ1 γ2 +γ1 γ3 +γ2 γ3 ) ⎜⎜⎝ ⎟⎟⎠ ⎜⎝ ∗0 ⎟⎠ (γ1 +γ2 )(γ2 +γ3 ) c0 1

20

⎞ ⎛ ⎛ ⎟⎟⎟ ⎜ a¯ ∗ ⎞⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ 0 ⎟⎟⎟ ⎟⎟⎟ = ⎜⎜⎜ b¯ ∗ ⎟⎟⎟ and ⎜⎜⎜⎜⎜ ⎟⎟⎠ ⎜⎝ ∗0 ⎟⎠ ⎜⎝ c¯ 0

γ 2 γ3 (γ1 +γ3 )(γ2 +γ3 ) −γ3 (γ1 +γ3 )

1

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ . ⎠

Furthermore, fuu (K0 , 0, 0)

=

−2γ1 , fuv (K0 , 0, 0) = 0, fvv (K0 , 0, 0) = 0, fuw (K0 , 0, 0) = −K0 γ1 , fvw (K0 , 0, 0) = 0, fww (K0 , 0, 0) = 0,

fuuu (K0 , 0, 0)

=

0, fvvv (K0 , 0, 0) = 0, fwww (K0 , 0, 0) = 0, fuuv (K0 , 0, 0) = 0, fuuw (K0 , 0, 0) = 0, fuvv (K0 , 0, 0) = 0,

fuvw (K0 , 0, 0)

=

0, fuww (K0 , 0, 0) = 0, fvvw (K0 , 0, 0) = 0, fvww (K0 , 0, 0) = 0,

guu (K0 , 0, 0)

=

0, guv (K0 , 0, 0) = γ2 , gvv (K0 , 0, 0) = −2γ2 , guw (K0 , 0, 0) = 0, gvw (K0 , 0, 0) = 0, gww (K0 , 0, 0) = 0,

guuu (K0 , 0, 0)

=

0, gvvv (K0 , 0, 0) = 0, gwww (K0 , 0, 0) = 0, guuv (K0 , 0, 0) = 0, guuw (K0 , 0, 0) = 0, guvv (K0 , 0, 0) = 0,

guvw (K0 , 0, 0)

=

0, guww (K0 , 0, 0) = 0, gvvw (K0 , 0, 0) = 0, gvww (K0 , 0, 0) = 0,

huu (K0 , 0, 0)

=

0, huv (k0, 0, 0) = 0, hvv (K0 , 0, 0) = 0, huw (K0 , 0, 0) = 0, hvw (K0 , 0, 0) = γ3 , hww (K0 , 0, 0) = −2γ3 ,

huuu (K0 , 0, 0)

=

0, hvvv (K0 , 0, 0) = 0, hwww (K0 , 0, 0) = 0 huuv (K0 , 0, 0) = 0, huuw (K0 , 0, 0) = 0, huvv (K0 , 0, 0) = 0,

huvw (K0 , 0, 0)

=

0, huww (K0 , 0, 0) = 0, hvvw (K0 , 0, 0) = 0, hvww (K0 , 0, 0) = 0.

Thus, we conclude that ⎞ ⎛ ⎛ ⎞ ⎜⎜⎜ −2γ1 x1 y1 − K0 γ1 x3 y1 − K0 γ1 x1 y3 ⎟⎟⎟ ⎜⎜⎜ 0 ⎟⎟⎟ ⎟⎟⎟ ⎜ ⎜ ⎟ B (X, Y) = ⎜⎜⎜⎜ γ2 x2 y1 + γ2 x1 y2 − 2γ2 x2 y2 ⎟⎟⎠ and C (X, Y, Z) = ⎜⎜⎜⎜⎝ 0 ⎟⎟⎟⎟⎠ . ⎝ γ3 x3 y2 + γ3 x2 y3 − 2γ3 x3 y3 0   Consequently, we obtain with all derivatives evaluated at (K0 , 0, 0, 0) B(q, q)

=

B(q, q) ¯ =

B(q, ¯ q) ¯ =

⎛ ⎜⎜⎜ ⎜⎜⎜ ⎜⎜⎝ ⎛ ⎜⎜⎜ ⎜⎜⎜ ⎜⎜⎜ ⎜⎝ ⎛ ⎜⎜⎜ ⎜⎜⎜ ⎜⎜⎜ ⎜⎝

⎞ −2K0 γ1 (a0 + K0 ) ⎟⎟⎟ ⎟⎟⎟ 2γ2 b0 (a0 − b0 ) ⎟⎟⎠ 2γ3 (b0 − 1) ¯ 0 − a¯e0 − K0 a0) γ1 (−2a 0e  γ2 b0 e¯ 0 + a0 b¯ 0 − 2b0 b¯ 0   γ3 b¯ 0 + b0 − 2 ⎞ −2γ1 a¯0 (¯a0 + K)  ⎟⎟⎟⎟ 2γ2 b¯ 0 a¯ 0 − b¯ 0 ⎟⎟⎟⎟ ⎟⎟⎠   2γ3 b¯ 0 − 1

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ ⎟⎠

⎛ ⎜⎜⎜ ⎜ = ⎜⎜⎜⎜ ⎝ ⎛ ⎜⎜⎜ ⎜ = ⎜⎜⎜⎜ ⎝ =

δ0 ε0 ζ0 ρ0 ξ0 ϑ0

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

⎞ ⎛ ⎜⎜⎜ δ¯ 0 ⎟⎟⎟ ⎟ ⎜⎜⎜ ⎜⎜⎝ ε¯ 0 ⎟⎟⎟⎟⎠ . ¯ζ0

A straightforward calculation gives   = −2a0 a∗0 γ1 (a0 + K0 ) + 2γ2 b0 b∗0 (a0 − b0 ) + 2γ3 (b0 − 1) q∗ , Bqq       q∗ , Bqq¯ = γ1 a∗0 (−2a0 a¯ 0 − K0 a¯ 0 − K0 a0 ) + γ2 b∗0 b0 a¯ 0 + a0 b¯ 0 − 2b0 b¯ 0 + γ3 b¯ 0 + b0 − 2   = −2a0 a¯ ∗0 γ1 (a0 + K0 ) + 2γ2 b0 b¯ ∗0 (a0 − b0 ) + 2γ3 (b0 − 1) q¯ ∗ , Bqq       q¯ ∗ , Bqq¯ = γ1 a¯ ∗0 (−2a0 a¯ 0 − K0 a¯ 0 − K0 a0 ) + γ2 b¯ ∗0 b0 a¯ 0 + a0 b¯ 0 − 2b0 b¯ 0 + γ3 b¯ 0 + b0 − 2       = −2γ1 a¯ 0 a∗0 (¯a0 + K0 ) + 2γ2 b¯ 0 b∗0 a¯ 0 − b¯ 0 + 2γ3 b¯ 0 − 1 q∗ , Bq¯ q¯       = −2γ1 a¯ 0 a¯ ∗0 (¯a0 + K0 ) + 2γ2 b¯ 0 b¯ ∗0 a¯ 0 − b¯ 0 + 2γ3 b¯ 0 − 1 q¯ ∗ , Bq¯q¯   q∗ , Cqqq¯ = 0. Thus, we calculate H20 and H11 as follow     ⎛ ⎜⎜⎜ 2γ1 (a0 a0 + K0 ) −1 + a0 a∗0 + a¯ 0 a¯ ∗0 − 2γ2 b0 (a0 − b0 ) a0 b∗0 + a¯ 0 b¯ ∗0 − 2γ3 (b0 − 1) (a0 + a¯ 0 )       ⎜⎜ H20 = ⎜⎜⎜⎜⎜ 2γ1 a0 (a0 + K0 ) b0 a∗0 + a¯ ∗0 b¯ 0 + 2γ2 b0 (a0 − b0 ) 1 − b0 b∗0 − b¯ ∗0 b¯ 0 − 2γ3 (b0 − 1) b0 + b¯ 0     ⎜⎝ −2γ3 (b0 − 1) + 2a0 γ1 (a0 + K0 ) a∗0 + a¯ ∗0 − 2γ2 (a0 − b0 ) b0 b∗0 + b0 b¯ ∗0 21

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

       ⎛ ⎜⎜⎜ γ1 (−2a0 a¯ 0 − K0 a¯ 0 − K0 a0 ) 1 − a0 a∗0 − a¯ 0 a¯ ∗0 − γ2 b0 a¯ 0 + (a0 − 2b0 ) b¯ 0 a0 b∗0 + a¯ 0 b¯ ∗0 − γ3 b¯ 0 + b0 − 2 (a0 + a¯ 0 )         ⎜⎜⎜ H11 = ⎜⎜⎜⎜ −γ1 (−2a0 a¯ 0 − K0 a¯ 0 − K0 a0 ) a∗0 b0 + a¯ ∗0 b¯ 0 + γ2 b0 a¯ 0 + (a0 − 2b0 ) b¯ 0 1 − b0 b∗0 − b¯ ∗0 b¯ 0 − γ3 b¯ 0 + b0 − 2 b0 + b¯ 0         ⎜⎝ γ1 (2a0 a¯ 0 + K0 a¯ 0 + K0 a0 ) a∗0 + a¯ ∗0 − γ2 b0 a¯ 0 + (a0 − 2b0 ) b¯ 0 b∗0 + b¯ ∗0 − γ3 b¯ 0 + b0 − 2

Since 2iω0 is not an eigenvalue of L0 (K0 ), the operators

where

⎛ ⎜⎜⎜ −γ1 ⎜ L0 (K0 ) = ⎜⎜⎜⎜ γ2 ⎝ 0

0 −γ2 γ3

L0 (K0 ) 2

− iω0 I and

r20

=

−(L0 (K0 ) − 2iω0 I)−1 H20 ,

r11

=

−(L0 (K0 ))−1 H11 ,

r02

=

−(L0 (K0 ) + 2iω0 I)−1 H02 .

−K0 γ1 0 −γ3

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ,

L0−1 (K0 )

1 = (1 + K0 ) γ1 γ2 γ3

⎛ ⎜⎜⎜ −γ1 − 2iω0 ⎜ γ2 (L0 (K0 ) − 2iω0 I) = ⎜⎜⎜⎜ ⎝ 0

0 −γ2 − 2iω0 γ3

⎛ ⎜⎜⎜ − (γ3 + 2iω0 ) (γ2 + 2iω0 ) −γ2 (γ3 + 2iω0 ) (L(K0 ) − 2iω0 I)−1 = (γ +2iω )(γ +2iω )(1γ +2iω )+K γ γ γ ⎜⎜⎜⎜⎝ 3 0 2 0 1 0 0 1 2 3 −γ2 γ3

L0 (K0 ) 2

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ . ⎟⎟⎠

+ iω0 I are invertibles and

⎛ ⎜⎜⎜ −γ2 γ3 ⎜⎜⎜ ⎜⎜⎝ −γ2 γ3 −γ2 γ3

K0 γ1 γ3 −γ1 γ3 −γ1 γ3

−K0 γ1 0 −γ3 − 2iω0

K0 γ1 γ2 K0 γ1 γ2 −γ1 γ2

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ,

γ3 (γ1 + 2iω0 ) − (γ3 + 2iω0 ) (γ1 + 2iω0 ) −γ3 (γ1 + 2iω0 )

K0 γ1 (γ2 + 2iω0 ) (γ2 + 2iω0 ) (γ1 + 2iω0 ) − (γ2 + 2iω0 ) (γ1 + 2iω0 )

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎠

and ⎛ ⎜⎜⎜ (γ3 + 2iω0 ) (γ2 + 2iω0 ) γ2 (γ3 + 2iω0 ) (2iω0 I − L(K0 ))−1 = (γ +2iω )(γ +2iω )(1γ +2iω )+K γ γ γ ⎜⎜⎜⎜⎝ 3 0 2 0 1 0 0 1 2 3 γ2 γ3

−γ3 (γ1 + 2iω0 ) (γ3 + 2iω0 ) (γ1 + 2iω0 ) γ3 (γ1 + 2iω0 )

−K0 γ1 (γ2 + 2iω0 ) − (γ2 + 2iω0 ) (γ1 + 2iω0 ) (γ2 + 2iω0 ) (γ1 + 2iω0 )

⎞ ⎟⎟⎟ ⎟⎟⎟ . ⎟⎠

Finally, with the above expressions we can calculate   c1 (K0 ) = ωi0 g20 g11 − |g11 |2 − 13 |g02 |2 + g221 = 2ωi 0 q∗ , Bqq  · q∗ , Bqq¯  + q∗ , Br11 q  + 12 q∗ , Br20 q¯ . if

Hence, according with Theorem 3 and Proposition 8 the Hopf bifurcation is supercritical (respectively, subcritical) l1 (K0 )  < 0 (respectively, l1 (K0 )  > 0). Furthermore, the bifurcating spatially homogeneous periodic solutions

d dK

α0 (K0 )

d dK

α0 (K0 )

from K = K0 are locally asymptotically stable (respectively unstable) if l1 (K0 ) < 0 (respectively, l1 (K0 ) > 0). The following example highlight the analytical results obtained in Theorem 5 for the system (1)-(2). To this purpose, we will follow the framework of section 2. Example 1. We choose the following values for the parameters in system (1): γ1 = 2, γ2 = 3, γ3 = 4, m = 0,  = 1, d1 = 0.5, d2 = 0.6 and d3 = 0.7. Thus, we obtain:    γ1 +γ3 γ2 +γ3 d 2 Re α + + + 2 = 8.75; where K is given by (47); (K ) = −0.05608, where αm (K) is K0 = γ1γ+γ m 0 0 γ γ dK 3 2 1 given by (45); ⎞ ⎞ ⎛ ⎞ ⎛ ⎛ ⎞ ⎛ ⎜⎜⎜ −1.1667 + 2.9744i ⎟⎟⎟ ⎜⎜⎜ a0 ⎟⎟⎟ ⎜⎜⎜ −1.1667 − 2.9744i ⎟⎟⎟ ⎜⎜⎜ a¯ 0 ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ ¯ ⎟⎟⎟ 1 + 1.2748i 1 − 1.2748i q = ⎜⎜ q¯ = ⎜⎜ ⎟⎟⎠ = ⎜⎜⎝ b0 ⎟⎟⎠ , ⎟⎟⎠ = ⎜⎜⎝ b0 ⎟⎟⎠ , ⎝ ⎝ c 1 1 ⎛ ⎛ ⎞ ⎛ 0 ⎞ ⎞ ⎛ c0 ⎞ ⎜⎜⎜ −403.20 − 513.98i ⎟⎟⎟ ⎜⎜⎜ a∗0 ⎟⎟⎟ ⎜⎜⎜ −403.20 + 513.98i ⎟⎟⎟ ⎜⎜⎜ a¯ ∗0 ⎟⎟⎟ ⎜ ⎜ ⎟ ⎟ ⎟ ⎜ ⎟ ⎜ q∗ = ⎜⎜⎜⎜ 16.80 − 28.55i ⎟⎟⎟⎟ = ⎜⎜⎜⎜ b∗0 ⎟⎟⎟⎟ , q¯ ∗ = ⎜⎜⎜⎜ 16.80 + 28.55i ⎟⎟⎟⎟ = ⎜⎜⎜⎜ b¯ ∗0 ⎟⎟⎟⎟ . ⎝ ⎝ ⎠ ⎝ ∗ ⎠ ⎠ ⎝ ∗ ⎠ c0 c¯ 0 1 1 22

⎛ ⎜⎜⎜ 70.7778 − 76.3437i ⎜ = ⎜⎜⎜⎜ −26.0000 − 6.3738i ⎝ 10.1980i

Furthermore, we get Bqq   q∗ , Bqq   q¯ ∗ , Bqq   q∗ , Bq¯ q¯   q∗ , Cqqq¯

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ,

Bqq¯

= 104 (1.0447 + 6.6321i) , = 103 (−68.396 − 4.9511i) , = 103 (−68.396 + 4.9511i) , = 0.

⎛ ⎜⎜⎜ 1. 523 × 10−3 ⎜ 0 = ⎜⎜⎜⎜ ⎝ 0

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ,

  q∗ , Bqq¯ = 10−12 (−2.8649 + 3.6521i) ,   q¯ ∗ , Bqq¯ = 10−12 (−2.8649 − 3.6521i) ,   q¯ ∗ , Bq¯ q¯ = 103 (−68.396 − 4.9511i) ,

To obtain r20 and r11, note that = =

H20 H11

106 (0.0236 − 0.1101i, 1.7958 − 1.3426i, 0.0574 − 0.0631i)T , 10−9 (−0.0284, 0.1506 + 0.0278i, 0.0057)T .

Thus, we get r20 r11

= =

1021 (−5.6931 − 0.0877i, −1.5023 + 2.4657i, 0.6251 + 1.6688i)T , 10−10 (−0.4780 + 0.0830i, 0.0241 − 0.0095i, 0.0384 − 0.0095i)T ,

and hence we obtain Br20 q¯ Br11 q

= =

1023 (1.7165 − 0.7039i, −0.0562 + 0.4737i, −0.0851 − 0.0319i)T , 10−9 (0.6424 + 0.4948i, −0.1357 − 0.2137i, −0.0057)T , 

q∗ , Br20 q¯  ∗ q , Br11 q

1 2



= 1025 (−1.7241 + 5.8619i) , = 10−7 (−5.0952 + 1.2325i) .

With the above expression we get from (36)        c1 (k0 ) = ωi0 q∗ , Bqq q∗ , Bqq¯ + q∗ , Br11 q + 12 q∗ , Br20 q¯ = 1025 (−1.7241 + 5.8619i) . Hence, according with Theorem 5 the Hopf bifurcation is subcritical and the bifurcating periodic solution is locally asymptotically stable since l1 (K0 )  > 0 and l1 (K0 ) < 0. d dK

α0 (K0 )

The bifurcation direction and stability of the spatially non-homogeneous bifurcating periodic solutions are determined in the following theorem. Theorem 6. The Hopf bifurcation for the systems (41)-(42) at K = K j , j = 1, 2, ..., m, ...., is supercritial (respectively, l (K ) l (K ) subcritical) if 1 j  < 0 (respectively, if 1 j  > 0). Furthermore, the bifurcating spatially non-homogeneous d dK

α j (K j )

d dK

α0 (K j )

periodic solutions are stable (respectively, unstable) if l1 (K j ) < 0 (respectively, l1 (K j ) > 0). We will omit the proof of Theorem 6 since the procedure described in Theorem 5 can be adapted in connection with the determination of the first Lyapunov coefficient for any m ∈ N − {0}, that is, for the non-homogeneous Hopf bifurcation, with minor modifications. instead of this, in the next subsection we present a numerical example highlighting the occurrence of non-homogeneous Hopf bifurcation in the system (1)-(2). Furthermore, we include a general algorithm conducted in MATLAB to calculate the first Lyapunov coefficient for any m ∈ N − {0}. 3.3. Numerical example for the non-homogeneous hopf bifurcation and the general algorithm The following example highlight the analytical results of Theorem 6 for the system (1)-(2). To this purpose, we will follow the framework of section 2. 23

Example 2. We choose the following values for the parameters in system (1): γ1 = 2, γ2 = 3, γ3 = 4, m = 1,  = 1, d1 = 0.5, d2 = 0.6 and d3 = 0.7. Thus, we obtain: K1 = 15.189, where Km is given by (47); ⎛ ⎞ ⎜⎜⎜ −2.5 0 −30.378 ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ −3.6 0 L1 (K1 ) = ⎜⎜ 3 ⎟⎟⎠ ; ⎝ 0 4 −4.7 ±iω0 = ±i6.1376;

δ1 (K1 ) = −10.8;

⎛ ⎜⎜⎜ 0.9034 ⎜ q = cos x ⎜⎜⎜⎜ 0.1927 − 0.3285i ⎝ −0.0743 − 0.1825i

d dK



  Re α1 (K1 ) = −0.0388, where αm (K) is given by (45);

⎞ ⎞ ⎛ ⎛ ⎟⎟⎟ ⎜⎜⎜ ⎜⎜⎜ a1 ⎟⎟⎟ 0.9034 ⎟⎟⎟ ⎟⎟⎟ ⎜ ⎜⎜⎜ ⎟⎟⎠ = cos x ⎜⎜⎝ b1 ⎟⎟⎠ ; q¯ = cos x ⎜⎜⎜⎜⎝ 0.1927 + 0.3285i c1 −0.0743 + 0.1825i

⎞ ⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ a¯ 1 ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜ ¯ ⎟⎟⎟ ⎟⎟⎠ = cos x ⎜⎜⎝ b1 ⎟⎟⎠ . c¯ 1

On the other hand we have ⎛ ⎜⎜⎜ −2.5 ⎜ L1∗ (K) = ⎜⎜⎜⎜ 0 ⎝ −2K

3 0 −3.6 4 0 −4.7

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

The operator L1∗ (K1 ) has the same eigenvalues as the eigenvalues of L1 (K1 ), and eigenvectors given by ⎛ ⎛ ⎞ ⎞ ⎛ ∗ ⎞ ⎛ ⎜⎜⎜ −0.1317 − 0.1719i ⎟⎟⎟ ⎜⎜⎜ −0.1317 + 0.1719i ⎟⎟⎟ ⎜⎜⎜ a1 ⎟⎟⎟ ⎜⎜⎜ ⎜ ⎜ ⎟ ⎟ ⎟ ⎜ ⎜ q∗ = cos x ⎜⎜⎜⎜ 0.2420 − 0.4127i ⎟⎟⎟⎟ = cos x ⎜⎜⎜⎜ b∗1 ⎟⎟⎟⎟ ; q¯ ∗ = cos x ⎜⎜⎜⎜ 0.2420 + 0.4127i ⎟⎟⎟⎟ = cos x ⎜⎜⎜⎜ ⎝ ⎝ ⎠ ⎠ ⎝ ∗ ⎠ ⎝ c1 0.8510 0.8510 Furthermore, we get ⎛ ⎜⎜⎜ 0.8161 − 10.0177i 2 ⎜ Bqq = cos x ⎜⎜⎜⎜ 1.4693 + 1.0210i ⎝ −0.3720 + 0.3031i

a¯ ∗1 b¯ ∗1 c¯ ∗1

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

⎛ ⎛ ⎛ ⎞ ⎞ ⎞ ⎞ ⎜⎜⎜ δ1 ⎟⎟⎟ ⎜⎜⎜ 0.8161 ⎟⎟⎟ ⎜⎜⎜ ρ1 ⎟⎟⎟ ⎟⎟⎟ ⎜ ⎜ ⎟ ⎟ ⎟⎟⎟ ⎟ 2 ⎜ 2 2 ⎟⎟⎠ = cos x ⎜⎜⎜⎜⎝ ε1 ⎟⎟⎟⎟⎠ ; Bqq¯ = cos x ⎜⎜⎜⎜⎝ 0.1741 ⎟⎟⎟⎟⎠ = cos x ⎜⎜⎜⎜⎝ ξ1 ⎟⎟⎟⎟⎠ ; ζ1 ϑ1 0.0544

and ⎛ ⎛ ⎞ ⎞ ⎜⎜⎜ 0.8161 ⎟⎟⎟ ⎜⎜⎜ ρ1 ⎟⎟⎟ ⎜ ⎜ ⎟ ⎟ Bqq¯ = cos2 x ⎜⎜⎜⎜ 0.1741 ⎟⎟⎟⎟ = cos2 x ⎜⎜⎜⎜ ξ1 ⎟⎟⎟⎟ . ⎝ ⎝ ⎠ ⎠ ϑ1 0.0544

To obtain r20 and r11 , note that ⎛ ⎛ ⎞ ⎜⎜⎜ −4 ⎜⎜⎜ 4 + 12.2752i 0 −30.378 ⎟⎟⎟ ⎜⎜⎜ ⎜⎜⎜ ⎟⎟⎟ 0 −3 L2 (K1 ) = ⎜⎜ 3 −5.4 − L (K ) = ; 2iω ⎜⎜⎝ 0 2 1 ⎟ ⎟⎠ ⎝ 0 4 −6.8 0

(2iω0 − L2 (K1 )) ⎛ ⎜⎜⎜ −2 ⎜ L0 (K1 ) = ⎜⎜⎜⎜ 3 ⎝ 0

0 −3 4

−1

0 5.4 + 12.2752i −4

⎛ ⎜⎜⎜ 0.0225 − 0.0873i 0.0538 − 0.0242i ⎜ = ⎜⎜⎜⎜ −0.0158 − 0.0125i 0.0299 − 0.0815i ⎝ −0.0053 + 0.0024i −0.0169 − 0.0185i

−30.378 0 −4

0.1470 + 0.1324i 0.0403 − 0.0182i 0.0341 − 0.0768i

⎛ ⎞ ⎜⎜⎜ 2 + 12.2752i ⎟⎟⎟ 0 ⎜ ⎟⎟⎟ −3 3 + 12.2752i ⎟⎟⎠ ; 2iω0 − L0 (K1 ) = ⎜⎜⎜⎜⎝ 0 −4 24

30.378 0 6.8 + 12.2752i ⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ;

30.378 0 4 + 12.2752i

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ;

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ ;

(2iω0 − L0 (K1 ))−1

⎛ ⎜⎜⎜ 0.0008 − 0.0901i 0.0359 − 0.0567i ⎜ = ⎜⎜⎜⎜ −0.0207 − 0.0053i 0.0077 − 0.0883i ⎝ −0.0035 + 0.0056i −0.0253 − 0.0108i

0.2010 + 0.0675i 0.0269 − 0.0425i 0.0140 − 0.0857i

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

Consequently, we obtain ⎞ ⎛  ⎜⎜⎜⎜ δ1 ⎟⎟⎟⎟ ⎜ cos (2x) (2iω0 − L2 (K1 )) + (2iω0 − L0 (K1 )) ⎜⎜⎜ ε1 ⎟⎟⎟⎟ ⎠ ⎝ ζ1 ⎞ ⎛ ⎛ ⎞ ⎜⎜⎜ 0.8437 − 0.2713i ⎟⎟⎟ ⎜⎜⎜ 0.9324 − 0.0300i ⎟⎟⎟ ⎜ ⎜ ⎟ ⎟ 1 (2x) ⎜⎜⎜⎜ 0.0523 − 0.3246i ⎟⎟⎟⎟ + 12 ⎜⎜⎜⎜ −0.0659 − 0.3420i ⎟⎟⎟⎟ ; 2 cos ⎝ ⎝ ⎠ ⎠ −0.1383 + 0.0067i −0.1080 − 0.0430i

=

r20

1 2

=

r11



= =

−1

−1

⎞ ⎛  ⎜⎜⎜⎜ ρ1 ⎟⎟⎟⎟ − 12 cos (2x) (L2 (K1 ))−1 + (L0 (K1 ))−1 ⎜⎜⎜⎜ ξ1 ⎟⎟⎟⎟ ⎠ ⎝ ϑ1⎞ ⎛ ⎛ ⎞ ⎜⎜⎜ 0.0420 ⎟⎟⎟ ⎜⎜⎜ 0.0037 ⎟⎟⎟ ⎜ ⎜ ⎟ ⎟ − 12 cos (2x) ⎜⎜⎜⎜ −0.0302 ⎟⎟⎟⎟ − 12 ⎜⎜⎜⎜ −0.0160 ⎟⎟⎟⎟ . ⎝ ⎝ ⎠ ⎠ −0.0296 −0.0274 

Using the formula given in (21) we get ⎞ ⎛ ⎟⎟⎟ ⎜⎜⎜ cos(2x) cos x (0.750 1i − 3. 663 9) + cos x (−2. 763 7 + 1. 402 3i) ⎟⎟⎟ ⎜⎜ ⎜ (0.634 x 19 + 0.690 01i) Br20 q¯ = ⎟ ⎜⎜⎜ cos(2x) cos x (−0.782 12i − 0.392 90) − cos ⎝ cos(2x) cos x 0.211 79 − 8. 453 5 × 10−2 i + cos x (0.206 24 − 0.178 84i) ⎟⎟⎠ ⎞ ⎛ ⎛ ⎞ ⎟⎟⎟ ⎜⎜⎜ −2. 763 7 + 1. 402 3i ⎟⎟⎟ ⎜⎜⎜ −3. 663 9 + 0.750 1i ⎟⎟ ⎜⎜ ⎜⎜ ⎟ ⎟ ⎜ ⎜ = cos(2x) cos x ⎜⎜ −0.392 90 − 0.782 12i ⎟⎟ + cos x ⎜⎜ −0.634 19 − 0.690 01i ⎟⎟⎟⎟ ; ⎝ ⎝ ⎠ −2 ⎠ 0.206 24 − 0.178 84i 0.211 79 − 8. 453 5 × 10 i

Br11 q

=

=

  ⎛ ⎜⎜⎜ cos(2x) cos x 1. 030 5 × 10−2 i − 0.372 82 − cos x (0.378 1 + 0.116 41i)   ⎜⎜⎜ ⎜⎜⎜ cos(2x) cos x(2. 236 8 × 10−2 − 3. 157 0 × 10−2 i) + cos x 3. 206 8 × 10−4 − 3. 649 8 × 10−2 i ⎜⎜⎝     cos(2x) cos x 9. 017 0 × 10−3 i + 1. 419 1 × 10−2 + cos x 1. 784 6 × 10−2 + 3. 690 3 × 10−3 i ⎞ ⎛ ⎛ −2 ⎜⎜⎜ −0.378 1 − 0.116 41i ⎜⎜⎜⎜ −0.372 82 + 1. 030 5 × 10 i ⎟⎟⎟⎟ ⎜ −2 ⎜ ⎜ cos(2x) cos x ⎜⎜ 2. 236 8 × 10 − 3. 157 0 × 10−2 i ⎟⎟⎟⎟ + cos x ⎜⎜⎜⎜ 3. 206 8 × 10−4 − 3. 649 8 × 10−2 i ⎠ ⎝ ⎝ 1. 419 1 × 10−2 + 9. 017 0 × 10−3 i 1. 784 6 × 10−2 + 3. 690 3 × 10−3 i

Finally, noticing that 1 ∗ 2 q , Br20 q¯ 

⎛ ⎞ ⎞ ⎛ ⎜⎜⎜ −0.1317 + 0.1719i ⎟⎟⎟ ⎜⎜⎜ ⎟⎟⎟ π −3. 663 9 + 0.750 1i ⎜⎜⎜ ⎟⎟⎟ ⎟⎟⎟ ⎜⎜⎜ 2 0.2420 + 0.4127i −0.392 90 − 0.782 12i = · ⎜⎜⎝ ⎟⎟ cos(2x) cos x ⎟⎟⎠ ⎜⎜⎝ −2 ⎠ 0 0.8510 ⎛ ⎞ ⎛0.211 79 − 8. 453 5 × 10 ⎞ i ⎜⎜⎜ −0.1317 + 0.1719i ⎟⎟⎟ ⎜⎜⎜ −2. 763 7 + 1. 402 3i ⎟⎟⎟ π ⎜ ⎟ ⎜ ⎟ + 12 ⎜⎜⎜⎜ 0.2420 + 0.4127i ⎟⎟⎟⎟ · ⎜⎜⎜⎜ −0.634 19 − 0.690 01i ⎟⎟⎟⎟ cos2 x ⎝ ⎠ ⎝ ⎠ 0.8510 0.206 24 − 0.178 84i 0 π π   = 12 cos(2x) cos2 x (0.761 53 + 1. 152 0i) + 12 cos2 x (0.429 73 + 1. 240 7i) 1 2

=

0 1 (0.761 53 8π

0

+ 1. 152 0i) + 14 π (0.429 73 + 1. 240 7i) ,

25

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠ .

and

q∗ , Br11 q 

⎞ ⎛ ⎞ ⎛ ⎜⎜⎜ −0.1317 − 0.1719i ⎟⎟⎟ ⎜⎜⎜ −0.372 82 + 1. 030 5 × 10−2 i ⎟⎟⎟ π ⎟ ⎜⎜⎜ ⎟⎟⎟ ⎜⎜⎜ = ⎜⎜ 0.2420 − 0.4127i ⎟⎟ · ⎜⎜ 2. 236 8 × 10−2 − 3. 157 0 × 10−2 i ⎟⎟⎟⎟ cos(2x) cos2 x ⎝ ⎠ ⎝ −2 −3 ⎠ 0 0.8510 ⎛ ⎞ ⎛1. 419 1 × 10 + 9. 017 0 × 10 i ⎞ ⎟⎟⎟ π ⎜⎜⎜ −0.1317 − 0.1719i ⎟⎟⎟ ⎜⎜⎜ −0.378 1 − 0.116 41i ⎟ ⎜ ⎟ ⎜ + ⎜⎜⎜⎜ 0.2420 − 0.4127i ⎟⎟⎟⎟ · ⎜⎜⎜⎜ 3. 206 8 × 10−4 − 3. 649 8 × 10−2 i ⎟⎟⎟⎟ cos2 x ⎝ ⎠ ⎝ −2 −3 ⎠ 0 1. 784 6 × 10 + 3. 690 3 × 10 i 0.8510  π   π  −2 = 5. 533 2 × 10 − 5. 353 3 × 10−2 i cos(2x) cos2 xdx + 2. 998 7 × 10−2 − 7. 450 2 × 10−2 i cos2 x 0 0    = π4 5. 533 2 × 10−2 − 5. 353 3 × 10−2 i + π2 2. 998 7 × 10−2 − 7. 450 2 × 10−2 i ,

we get from (36) Re(c1 (K ∗ ))

1 ∗ ∗ = Re q  , Br11 q  + 2 Re q  , Br20 q¯   π = 4 5. 533 2 × 10−2 + π2 2. 998 7 × 10−2 + 18 π (0.761 53) + 14 π (0.429 73) = 0.231 45π > 0.

Hence, according with Theorem 5 the Hopf bifurcation is supercritical and the bifurcating periodic solutions are unstable since l1 (K1 )  < 0 and l1 (K1 ) > 0. d dK

α1 (K1 )

In the following, we include the MATLAB code used to simulate calculate the first lyapunov coefficient for the non-homogeneous Hopf bifurcation. The same code is appropriately modified for calculating this coefficient for any m ∈ N − {0}. # Declaring the Parameters composing the system (1) g1=input(’Enter with a value for gamma 1: ’); g2=input(’Enter with a value for gamma 2: ’); g3=input(’Enter with a value for gamma 3: ’); d1=input(’Enter with a value for d1: ’); d2=input(’Enter with a value for d2: ’); d3=input(’Enter with a value for d3: ’); l=input(’Enter with a value for : ’); m=input(’Enter with a value for m: ’); # Finding the bifurcation values K1=((g1+g2)*(g1+g3)*(g2+g3))/(g1*g2*g3) K2=((((d2+d3)*(g1+g3)+(d1+d3)*(g2+g3))*(g1+g2))/((g1*g2*g3)*lˆ2) +((g2+g3)*(g1+g3)*(d1+d2))/((g1*g2*g3)*lˆ2))*(mˆ2) K3=(((d2+d3)*(d1+d3)*(g1+g2))/(g1*g2*g3*lˆ4)+((d2+d3)*(g1+g3) +(g2+g3)*(d1+d3))*(d1+d2)/(g1*g2*g3*lˆ4))*(mˆ4) K4=(((d2+d3)*(d1+d3)*(d1+d2))/(g1*g2*g3*lˆ6))*(mˆ6) km=K1+K2+K3+K4 # Determining the coefficients for the characteristic polynomial Mm=(g1+g2+g3)+((d1+d2+d3)*(mˆ2)/(lˆ2)) Nm=((g1*g2)+(g1*g3)+(g2*g3))+((g2+g3)*d1+(g1+g3)*d2+(g1+g2)*d3)*(mˆ2)/(lˆ2)+ ((d2*d3)+(d1*d3)+(d1*d2))*(mˆ4)/(lˆ4) Pm=(1+km)*g1*g2*g3+(g2*g3*d1+g1*g3*d2+g1*g2*d3)*(mˆ2)/(lˆ2)+ (d2*d3*g1+d1*d3*g2+d1*d2*g3)*(mˆ4)/(lˆ4)+d1*d2*d3*mˆ6/lˆ6 26

# Determining the matrix Lm (K0 ) and ω0 A=[-g1-(d1*mˆ2)/lˆ2 0 -(km)*g1;g2 -g2-(d2*mˆ2)/lˆ2 0;0 g3 -g3-(d3*mˆ2)/lˆ2] G=[g1 g2 g3] [vecpA valpA]=eig(A) omega=imag(valpA(1,1)) ∗ # Obtaining the eigenvalues and eigenvectors associated to Lm (K0 ) and Lm (K0 )

q0=[vecpA(1,2) vecpA(2,2) vecpA(3,2)] bq0=[vecpA(1,3) vecpA(2,3) vecpA(3,3)] N=A’ [vecpN valpN]=eig(N) sq0=[vecpN(1,2) vecpN(2,2) vecpN(3,2)] bsq0=[vecpN(1,3) vecpN(2,3) vecpN(3,3)] # Calculating the vectorsBqq and Bqq¯ Bqq=[-g1*(2*q0(1,1)*q0(1,1)+(km)*(q0(1,3)*q0(1,1)+q0(1,1)*q0(1,3))) g2*(q0(1,2)*q0(1,1) +q0(1,1)*q0(1,2)-2*q0(1,2)*q0(1,2)) g3*(q0(1,3)*q0(1,2)+q0(1,2)*q0(1,3)-2*q0(1,3)*q0(1,3))] Bqbq=[-g1*(2*q0(1,1)*conj(q0(1,1))+(km)*(q0(1,3)*conj(q0(1,1)) +q0(1,1)*conj(q0(1,3)))) g2*(q0(1,2)*conj(q0(1,1))+q0(1,1)*conj(q0(1,2)) 2-*q0(1,2)*conj(q0(1,2))) g3*(q0(1,3)*conj(q0(1,2))+q0(1,2)*conj(q0(1,3))-2*q0(1,3)*conj(q0(1,3)))] # Obtaining the operators L2 (K0 ), L0 (K0 ), (2iω0 I − L2 (K0 ))−1 , (2iω0 I − L0 (K0 ))−1 L2=[-g1-(d1*4/lˆ2) 0 -(km)*g1;g2 -g2-(d2*4/lˆ2) 0;0 g3 -g3-(d2*4/lˆ2)] IL2=inv(L2) L2om=[2*1i*omega+g1+(d1*4/lˆ2) 0 (km)*g1;-g2 2*1i*omega+g2+(d2*4/lˆ2) 0;0 -g3 2*1i*omega+g3+(d2*4/lˆ2)] IL2om=inv(L2om) L0=[-g1 0 -(km)*g1;g2 -g2 0;0 g3 -g3] IL0=inv(L0) L0om=[2*1i*omega+g1 0 (km)*g1;-g2 2*1i*omega+g2 0;0 -g3 2*1i*omega+g3] IL0om=inv(L0om) # Calculating the vectors r20 , r11 w201=IL2om*Bqq’ w202=IL0om*Bqq’ w111=IL2*Bqbq’ w112=IL0*Bqbq’ r20=(1/2)*(cos(2*x)*w201+w202) r11=-(1/2)*(cos(2*x)*w111+w112) # Obtaining the vectors Br20 q¯ , Br11 q Br20bq0=[-g1*(2*r20(1,1)*bq0(1,1)+(km)*(r20(3,1)*bq0(1,1)+r20(1,1)*bq0(1,3))) g2*(r20(2,1)*bq0(1,1) +r20(1,1)*bq0(1,2)-2*r20(2,1)*bq0(1,2)) g3*(r20(3,1)*bq0(1,2)+r20(2,1)*bq0(1,3)-2*r20(3,1)*bq0(1,3))] Br11q0=[-g1*(2*r11(1,1)*q0(1,1)+(km)*(r11(3,1)*q0(1,1)+r11(1,1)*q0(1,3))) g2*(r11(2,1)*q0(1,1)+ r11(1,1)*q0(1,2)-2*r11(2,1)*q0(1,2)) g3*(r11(3,1)*q0(1,2)+r11(2,1)*q0(1,3)-2*r11(3,1)*q0(1,3))]     # Calculating the vectors q∗ , Br20 q¯ , q∗ , Br11 q and the Lyapunov coefficient sqBr20bq0=(1/2)*sq0*Br20bq0’ sqBr11q0=sq0*Br11q0’ Real1=real(sqBr20bq0) Real2=real(sqBr11q0) 27

4. Simulation results In this section we present some simulations carried in MATLAB to verify the veracity of the analytical results obtained for system (1)-(2). If we consider system (1) with values of parameters: γ1 = 2, γ2 = 3, γ3 = 4, d1 = 0.6, d2 = 0.6 and d3 = 0.7 we get K1 < K = 15.6524, where K1 = 15.189. In this case, we obtain Figures 1 to 3 illustrating that the solutions with initial data close to E, that is, u0 (x) = v0 (x) = w0 (x) = 1.1, approximate a periodic function when t becomes large. u(x,t)

7 6 5 4 3 2 1 0 5 4.5 4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5 1

1 0.5

0.5 0

0

Time t

Distance x

Figure 1: Numerical simulations of the system (1)-(2): the non-homogeneous case for the function u(t,x).

v(x,t)

4 3.5 3 2.5 2 1.5 1 0.5 0 5 4.5 4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5 1

1 0.5

0.5 0

0

Time t

Distance x

Figure 2: Numerical simulations of the system (1)-(2): the non-homogeneous case for the function v(t,x).

5. Conclusions and remarks It is well know that more realist populational models are those describing the dynamics of populations that live in a 1, 2 or 3 dimensional space. In these cases the dynamics is described, for example, by a system of reaction-diffusion equations, that is, a nonlinear or rather quasilinear system of parabolic partial differential equations. In response to the need of considering and analyzing temporal and spatial variation in mathematical models describing the dynamics interaction among groups of individuals, in the Section 2 of this paper we analyzed the general reaction-diffusion (4)-(6) subject to Neumann boundary conditions on spatial domain Ω = (0; lπ), with l ∈ R+ . We display expressions for the first and second Lyapunov coefficients for the system (4)-(6), pushing forward the method found in [11]. The formula derived for the first and second Lyapunov coefficients allow us to verify the non-degeneracy conditions for the codimension 1 and 2 Hopf bifurcation. Furthermore, this allows to determine the number, types and positions of bifurcating small amplitude periodic orbits. The procedure to obtain the first and second Lyapunov coefficients can be adapted in connection for the determination of the third and fourth Lyapunov coefficients with few modifications. 28

w(x,t)

2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 5 4.5 4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5 1

1 0.5

0.5 0

0

Time t

Distance x

Figure 3: Numerical simulations of the system (1)-(2): the non-homogeneous case for the function w(t,x).

As an application of the Lyapunov coefficient formula, the predator-prey system (1)-(2) is analyzed. This system is a special subset of three-dimensional Lotka-Volterra dynamical system with diffusion which we can interpret biologically as follow. In the system (1)-(2) , the first equation shows that w is a predator and u is a prey. In the second equation , v is a predator and u is a prey. Finally, in the third equation, w is a predator and v is a prey. In Section 3 we showed that the system (1)-(2) undergoes a codimension 1 Hopf bifurcation; in other words, we carry out a Hopf bifurcation analysis for the reaction-diffusion model (1)-(2) and show that the bifurcating periodic orbits are either spatially homogeneous, which coincides with the periodic solution of the corresponding ODE system, or spatially non-homogeneous. The previous bifurcation result means that the species coexist in a periodic manner due to the occurrence of Andronov-Hopf bifurcation. It is worth observing that in previous studies the ODE version of system (1)-(2) was also analyzed and similar bifurcation results were obtained in [8] (see also [26, 27]). It shows that diffusion is a process that generally helps to maintain the qualitative behavior of a dynamical system. Consequently, the conclusions obtained for the ODE version of system (1)-(2) related to the coexistence among the species sustain for the spatial version. Finally, we would like to stress that the cascade of Hopf and steady state bifurcations shown in this paper represents the rich self-organized spatiotemporal dynamics of the diffusive predator-prey system (1)-(2). Furthermore, although this work ultimately focuses the specific three dimensional system of partial differential equations given by (1)-(2), the method of analysis and calculations explained in Section 3 can be adapted to the study of other systems with three or more phase variables and depending on more parameters. Acknowledgement The authors wish to thank the anonymous reviewers for their constructive suggestions, which have led to the improvement of our earlier presentation. Jocirei D. Ferreira acknowledges the funding support from CNPq-Brazil under grants No. 445356/2014-6 and 200628/2014-3/PDE, and from INCTMat/CNPq. Aida P. G. Nieva and Wilmer M. Yepez acknowledge support from the Departamento de Matem´aticas, Universidad del Cauca, Popay´an, Colombia. [1] A. A. Andronov, E. A. Leontovich et al., Theory of Bifurcations of Dynamic Systems on a Plane. Halsted Press, J. Wiley & Sons, New York, (1973). [2] A. Bocs´o, M. Farkas, Political and economic rationality leads to velcro bifurcation. Applied Mathematics and Computation, 140, pp. 381-389, (2003). [3] A. Gasull and J. Torregrosa, A new approach to the computation of the Lyapunov Constants. Comp. and Appl. Math., 20 , pp. 149-177 (2001). [4] A.M. Turing, The chemical basis of morphogenesis. Phil. Trans. Royal Soc. London, Series B, Biological Sciences, Vol. 237, No. 641, pp. 37-72, (1952). [5] A. Okubo, Diffusion and Ecological Problems: Mathematical Models. Springer-Verlag, Berlin, (1980). [6] B.D Hassard, N.D Kazarinoff, Y.H Wan Theory and Application of Hopf Bifurcation. Cambridge Univ. Press, Cambridge, (1981). [7] D. Tilman and P. Kareiva, Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions. Monographs in Population Biology, 30, Princeton Univ. press, NJ, USA, (1997). [8] E. Hadidi, Bifurcation of Limit Cycle for Three-Dimensional Lotka-Volterra Dynamical System. Applied Mathematical Sciences, 7, no. 139, 6909-6916, (2013).

29

[9] E. S´aez, E. Stange, and I. Sz´ant´o, Simultaneous Zip Bifurcation and Limit Cycles in Three Dimensional Competition Models, SIAM J. Applied Dynamical Systems, 05, 1-11, (2006). [10] F. Takens, Unfoldings of certain singularities of vectorfields: Generalized Hopf bifurcations, J. Diff. Equat., 14, 476 - 493, (1973). [11] F. Yi, J. Wei and J. shi. Bifurcatyion and Spatitemporal Patterns in a Homogeneous Diffusive predator-prey system. Journal of Differential Equations, Vol. 246, pp. 1944-1977, (2009). [12] J. D. Ferreira and L. A. F. de Oliveira, Hopf and zip bifurcation in a specific (n+1)-dimensional competitive system. Matem´aticas: Ense˜nanza Universitaria, Vol. 15, pp. 35-41, (2008). [13] J. D. Ferreira and L. A. F. de Oliveira, Zip Bifurcation in a Competitive System with Diffusion . Differential Equations and Dynamical Systems: An International Journal for Theory Applications & Computer Simulations, Vol. 17, pp 37-53, (2009). [14] J. D. Murray, Mathematical Biology I: An Introduction. Springer, USA, Third Edition, 2002. [15] J. D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, USA, Third Edition, 2003. [16] J. Y. Wakano, K. Ikeda, T. Miki and M. Mimura, Effective dispersal rate is a function of habitat size and corridor shape: mechanistic formulation of a two-patch compartment model for spatially continuous systems. Oikos, 120, 1712-1720, (2011). [17] J. Sotomayor, L.F. Mello and D. de Caravallo Braga, Lyaponuv Coefficients for Degenerate Hopf Bifurcations. arXiv:0709.3949v1 [math.DS], Sep. 25 (2007), home page: http://arxiv.org/abs/0709.3949. [18] J. Smoller, Shock waves and Reaction-Diffusion Equations. Springer-Verlag, New York, Second Edition, (1994). [19] M. Cavani and M. Farkas, Bifurcations in a Predator-prey Model with Memory and Diffusion II: Turing Bifurcation. Acta Math. Hungar., 63, 375-393, (1994). [20] M. Cristofol and L. Roques, Biological Invasions: Deriving the Regions at Risk from Partial Measurements. Institut National de la Recherche Agronomique and Math´ematiques et Informatique Appliqu´es, Research Report No. 32. http://www.biosp.org. [21] M. Farkas, Zip Bifurcation in a Competition Model. Nonlinear Analysis, Theory, Methods & Applications, 8, 1295-1309, (1984). [22] M. Farkas, Periodic Motions. Springer-Verlag, New York, (1994). [23] M. Farkas and J. D. Ferreira, Zip Bifurcation in a Reaction-Diffusion System. Differential Equations and Dynamical Systems-An International Journal for Theory, Applications and Computer Simulations, 15, 169-183, (2007). [24] M. Farkas, E. S´aez, and I. Sz´ant´o, Velcro bifurcation in competition models with generalized Holling functional response. Miskolc Mathematical Notes, 6, no. 2, 185–195, (2005). [25] P. A. Hancock and H. C. J. Godfray, Modelling the spread of Wolbachia in spatially heterogeneous environments. J. R. Soc. Interface, doi:10.1098/rsif.2012.0253, (2012). [26] S. M. Shahruz and D.A. Kalkin, Limit cycle behavior in three or higher- dimensional nonlinear systems : the Lotka-Volterra example. Journal of sound and vibration, 246, no. 2, 379-386, (2001). [27] S. M. Shahruz, Generation of self-pulsation in passively Q-switched lasers. Physica D, 142, 291-305, (2000). [28] W. Ko and K. Ryu, A qualitiative study on general Gause-type predator-prey models with constant diffusion rates. J. Math. Anal. and Appl., 344, 217-230, (2008). [29] Y. A. Kuznetsov, Elements of applied Bifurcation theory. Springer-Verlag, (1998).

Appendix 1 The complex numbers and the components of vectors involved in system (19) given in Remark 2 can be computed by the formulae:

Gi j

=

¯ 10,i G

=

¯ 01,i G

=

¯ 110,i G

=

¯ 210,i G

=

¯ 120,i G

=

¯ 011,i G

=

¯ 1110,i G

=

#

# ∂ j+i

q∗ , F(K0 , zq + z¯q¯ + r)## ∂z j ∂¯zi # ## z=0 ∂2 ∗ , F(K , zq + z ¯ q ¯ + r)

q #z=0,r=0 0 ∂ri ∂z #

#

q , F(K0 , zq + z¯q¯ + r)## ##z=0,r=0 ∂3

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z2 ##z=0,r=0 ∂3 ∗ ¯ + r)## ∂ri ∂z∂¯z q , F(K0 , zq + z¯q ##z=0,r=0 ∂3 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z∂z z=0,r=0 # # ∂3 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂¯z ##z=0,r=0 ∂3 ∗ ¯q¯ + r)## 2 q , F(K0 , zq + z ∂2 ∂ri ∂¯z

∗

∂ri ∂¯z

z=0,r=0

30

,

2≤i+ j≤5

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

¯ 1120,i G

=

¯ 1210,i G

=

¯ 1220,i G

=

¯ 2110,i G

=

¯ 2120,i G

=

¯ 2210,i G

=

¯ 2220,i G

=

¯ 11110,i G

=

¯ 11120,i G

=

¯ 11210,i G

=

¯ 11220,i G

=

¯ 12110,i G

=

¯ 12120,i G

=

¯ 12210,i G

=

¯ 12220,i G

=

¯ 21110,i G

=

¯ 21120,i G

=

¯ 21210,i G

=

¯ 21220,i G

=

¯ 22110,i G

=

¯ 22120,i G

=

¯ 22210,i G

=

¯ 22220,i G

=

¯ 111110,i G

=

¯ 111120,i G

=

¯ 111210,i G

=

¯ 112110,i G

=

¯ 121110,i G

=

¯ 211110,i G

=

¯ 111220,i G

=

¯ 112210,i G

=

¯ 121210,i G

=

¯ 211210,i G

=

#

# ∂3

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z2 z=0,r=0 ## ∂4 ∗ ¯ + r)## ∂ri ∂z∂¯z∂z q , F(K0 , zq + z¯q # z=0,r=0

# ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z2 ∂z ##z=0,r=0 ∂4 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂z ∂¯z ##z=0,r=0 ∂4 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z∂z∂¯z # ## z=0,r=0 ∂4 ∗ , F(K , zq + z ¯ q ¯ + r)

q # 0 ∂ri ∂z∂¯z2 ## z=0,r=0 ∂4 ∗ #

q , F(K0 , zq + z¯q¯ + r)# ∂ri ∂¯z3 ##z=0,r=0 ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 4 ∂ri ∂z z=0,r=0 ## 5 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 3 ∂ri ∂¯z∂z z=0,r=0 ## ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂z∂¯z∂z2 ## z=0,r=0 5 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z2 ∂z2 ##z=0,r=0 5 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂z ∂¯z∂z ##z=0,r=0 5 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z∂z∂¯z∂z # ## z=0,r=0 ∂5 ∗ , F(K , zq + z ¯ q ¯ + r)

q # 0 ∂ri ∂z∂¯z2 ∂z ## z=0,r=0 ∂5

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z3 ∂z ##z=0,r=0 5 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z3 ∂¯z z=0,r=0 ## ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂¯z∂z ∂¯z ##z=0,r=0 ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂z∂¯z∂z∂¯z # ## z=0,r=0 ∂5 ∗ , F(K , zq + z ¯ q ¯ + r)

q # 0 ∂ri ∂¯z2 ∂z∂¯z ## z=0,r=0 ∂5 ∗ #

q , F(K0 , zq + z¯q¯ + r)# ∂ri ∂z2 ∂¯z2 ##z=0,r=0 ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂¯z∂z∂¯z ## z=0,r=0 5 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z∂¯z3 ## z=0,r=0 ∂5 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z4 ##z=0,r=0 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z5 z=0,r=0 ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 4 ∂ri ∂¯z∂z z=0,r=0 ## 6 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 3 ∂ri ∂z∂¯z∂z ##z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂z2 ∂¯z∂z2 # # z=0,r=0 ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z3 ∂¯z∂z ## z=0,r=0 ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z4 ∂¯z ##z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 3 ∂ri ∂¯z ∂z z=0,r=0 ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 2 ∂ri ∂z∂¯z ∂z z=0,r=0 ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂z∂¯z∂z∂¯z∂z ## z=0,r=0 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z∂¯z∂z2 ∂¯z z=0,r=0 4

31

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

¯ 112120,i G

=

¯ 121120,i G

=

¯ 211120,i G

=

¯ 122110,i G

=

¯ 221110,i G

=

¯ 212110,i G

=

¯ 112220,i G

=

¯ 121220,i G

=

¯ 122120,i G

=

¯ 12221,i G

=

¯ 211220,i G

=

¯ 212120,i G

=

¯ 212210,i G

=

¯ 221120,i G

=

¯ 222110,i G

=

¯ 221210,i G

=

¯ 122220,i G

=

¯ 212220,i G

=

¯ 221220,i G

=

¯ 222120,i G

=

¯ 222210,i G

=

¯ 222220,i G

=

#

# ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z∂z∂¯z∂z2 ##z=0,r=0 ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z∂z2 ∂¯z∂z ## z=0,r=0 ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z∂z3 ∂¯z ##z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 2 ∂ri ∂z ∂¯z ∂z z=0,r=0 # # ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 3 2 ∂ri ∂z ∂¯z z=0,r=0 ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂z∂z∂¯z∂z∂¯z ## z=0,r=0 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z3 ∂z2 z=0,r=0 ## 6 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂¯z ∂z∂¯z∂z ##z=0,r=0 6 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂¯z∂z∂¯z ∂z # ## z=0,r=0 ∂6 ∗ , F(K , zq + z ¯ q ¯ + r)

q #z=0,r=0 0 ∂ri ∂z∂¯z3 ∂z ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z2 ∂z2 ∂¯z z=0,r=0 ## 6 ∂ ∗ ¯ + r)## ∂ri ∂¯z∂z∂¯z∂z∂¯z q , F(K0 , zq + z¯q ## z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 2 ∂ri ∂z∂¯z ∂z∂¯z # ## z=0,r=0 ∂6 ∗ , F(K , zq + z ¯ q ¯ + r)

q # 0 ∂ri ∂¯z∂z2 ∂¯z2 ## z=0,r=0 ∂6

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z2 ∂¯z3 z=0,r=0 ## 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂z∂¯z∂z∂¯z2 ## z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 4 ∂ri ∂¯z ∂z z=0,r=0 ## ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 3 ∂ri ∂¯z ∂z∂¯z ##z=0,r=0 ∂6 ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 ∂ri ∂¯z2 ∂z∂¯z2 ## z=0,r=0 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z∂z∂¯z3 ## z=0,r=0 6 ∂ ∗ ## , F(K , zq + z ¯ q ¯ + r)

q 0 4 ∂ri ∂z∂¯z ## z=0,r=0 6 ∂

q∗ , F(K0 , zq + z¯q¯ + r)## ∂ri ∂¯z5 z=0,r=0

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

,

i = 1, 2, 3

and Hi j

=

∂i+ j F(K0 , zq ∂zi ∂¯z j

+ z¯q) ¯ z=0 − Gi j q − G¯ i j q¯

, 2 ≤ i + j ≤ 5.

Appendix 2 In this appendix, we present the expressions for Bqq , Bqq¯ , Bq¯ q¯ , Cqqq , Cqqq¯ , Cqq¯ q¯ , Cq¯ q¯ q¯ , Dqqqq , Dqqqq¯ , Dqqq¯ q¯ , Dqq¯ q¯ q¯ , Dq¯ q¯ q¯ q¯ , Eqqqqq , Eqqqqq¯ , Eqqqq¯ q¯ , Eqqq¯ q¯ q¯ , Eqq¯ q¯ q¯ q¯ and Eq¯ q¯ q¯ q¯ q¯ , which are presented in Proposition 4. We have that:

32

Bqq

Bqq¯

⎛ $ mx % ⎜⎜⎜ δ1m ⎜⎜⎜ = cos2 ⎜ ε1m l ⎜⎝ ζ 1m

⎛ $ mx % ⎜⎜⎜ ρ1m ⎜⎜⎜ = cos2 ⎜ ξ1m l ⎜⎝ ϑ 1m

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

where

⎧ ⎪ δ1m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ε1m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ζ1m ⎪ ⎪ ⎪ ⎩ ⎧ ⎪ ρ1m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ξ1m ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ϑ1m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

= = =

=

=

=

fuu a2m + fvv b2m + fww c2m +2 fuv am bm + 2 fuw am cm + 2 fvw bm cm , guu a2m + gvv b2m + gww c2m +2guv am bm + 2guw am cm + 2gvw bm cm , huu a2m + hvv b2m + hww c2m +2huv am bm + 2hguw am cm + 2hvw bm cm ;   fuu |am |2 + fuv am b¯ m + a¯ m bm + fvv |bm |2 + fuw (a  m c¯ m + a¯ m cm ) + fww |cm |2 + fvw bm c¯ m + b¯ m cm ,   guu (K0 , 0) |am |2 + guv am b¯ m + a¯ m bm +gvv |bm |2 + guw (am c¯ m + a¯ m cm )  +gww |cm |2 + gvw (K0 , 0) bm c¯ m + b¯ m cm ,   huu |am |2 + huv am b¯ m + a¯ m bm +hvv |bm |2 + huw (a  m c¯ m + a¯ m cm ) 2 +hww |cm | + hvw bm c¯ m + b¯ m cm ;

(48)

¯ qq ; Bq¯ q¯ = B

Cqqq

⎛ $ mx % ⎜⎜⎜ δ2m ⎜⎜⎜ = cos3 ⎜ ε2m l ⎜⎝ ζ 2m

⎧ ⎪ δ2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎞ ⎪ ⎪ ⎪ ⎟⎟⎟ ⎪ ⎪ ⎨ ε2m ⎟⎟⎟ ⎟⎠⎟ where ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ζ2m ⎪ ⎪ ⎪ ⎩

Cqqq¯ ⎧ ⎪ ⎪ ρ2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ξ2m ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ϑ2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

=

fuuu a3m + fvvv b3m + fwww c3m + 3 fuuv a2m bm + 3 fuuw a2m cm + 3 fuvv b2m am +6 fuvw am bm cm + 3 fuww am c2m + 3 fvvw b2m cm + 3 fvww bm c2m , = guuu a3m + gvvv b3m + gwww c3m + 3guuv a2m bm + 3guuw a2m cm + 3guvv b2m am +6guvw am bm cm + 3guww am c2m + 3gvvw b2m cm + 3gvww bm c2m , = huuu a3m + hvvv b3m + hwww c3m + 3huuv a2m bm + 3huuw a2m cm + 3huvv b2m am +6huvw am bm cm + 3huww am c2m + 3hvvw b2m cm + 3hvww bm c2m ; ⎛ $ mx % ⎜⎜⎜ ρ2m ⎜⎜⎜ = cos3 ⎜ ξ2m l ⎜⎝ ϑ 2m

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

(49)

where

  fuuu am |am |2 + fvvv bm |bm |2 + fwww cm |cm |2 + fuuv 2bm am a¯ m + a2m b¯ m       + fuuw a2m c¯ m + 2am a¯ m cm + fuvv a¯ m b2m + 2am |bm |2 + fuvw 2¯am bm cm + 2am b¯ m cm + 2am bm c¯ m       + fuww c2m a¯ m + 2am |cm |2 + fvvw 2cm |bm |2 + b2m c¯ m + fvww c2m b¯ m + 2 |cm |2 bm   = guuu am |am |2 + gvvv bm |bm |2 + gwww cm |cm |2 + guuv 2bm am a¯ m + a2m b¯ m       +guuw a2m c¯ m + 2am a¯ m cm + guvv a¯ m b2m + 2am |bm |2 + guvw 2¯am bm cm + 2am b¯ m cm + 2am bm c¯ m       +guww c2m a¯ m + 2am |cm |2 + gvvw 2cm |bm |2 + b2m c¯ m + gvww c2m b¯ m + 2 |cm |2 bm   = huuu am |am |2 + hvvv bm |bm |2 + hwww cm |cm |2 + huuv 2bm am a¯ m + a2m b¯ m       +huuw a2m c¯ m + 2am a¯ m cm + huvv a¯ m b2m + 2am |bm |2 + huvw 2¯am bm cm + 2am b¯ m cm + 2am bm c¯ m       +huww c2m a¯ m + 2am |cm |2 + hvvw 2cm |bm |2 + b2m c¯ m + hvww c2m b¯ m + 2 |cm |2 bm ; =

Cqq¯ q¯

⎛ $ mx % ⎜⎜⎜ σ2m ⎜⎜⎜ = cos ⎜ τ2m l ⎜⎝  2m 3

33

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

(50)

⎧ ⎪ ⎪ σ2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ τ2m ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

Cq¯ q¯ q¯

  fuuu |am |2 a¯ m + fvvv |bm |2 b¯ m + fwww |cm |2 c¯ m + fuuv bm a¯ 2m + 2 |am |2 b¯ m       + fuuw 2 |am |2 c¯ m + cm a¯ 2m + fuvv 2 |bm |2 a¯ m + am b¯ 2m + fuvw 2¯am b¯ m cm + 2¯am bm c¯ m + 2am b¯ m c¯ m       + fuww 2 |cm |2 a¯ m + am c¯ 2m + fvvw cm b¯ 2m + 2 |bm |2 c¯ m + fvww 2 |cm |2 b¯ m + bm c¯ 2m   = guuu |am |2 a¯ m + gvvv |bm |2 b¯ m + gwww |cm |2 c¯ m + guuv bm a¯ 2m + 2 |am |2 b¯ m       +guuw 2 |am |2 c¯ m + cm a¯ 2m + guvv 2 |bm |2 a¯ m + am b¯ 2m + guvw 2¯am b¯ m cm + 2¯am bm c¯ m + 2am b¯ m c¯ m       +guww 2 |cm |2 a¯ m + am c¯ 2m + gvvw cm b¯ 2m + 2 |bm |2 c¯ m + gvww 2 |cm |2 b¯ m + bm c¯ 2m   = huuu |am |2 a¯ m + hvvv |bm |2 b¯ m + hwww |cm |2 c¯ m + huuv bm a¯ 2m + 2 |am |2 b¯ m       +huuw 2 |am |2 c¯ m + cm a¯ 2m + huvv 2 |bm |2 a¯ m + am b¯ 2m + huvw 2¯am b¯ m cm + 2¯am bm c¯ m + 2am b¯ m c¯ m       +huww 2 |cm |2 a¯ m + am c¯ 2m + hvvw cm b¯ 2m + 2 |bm |2 c¯ m + hvww 2 |cm |2 b¯ m + bm c¯ 2m ; =

⎛ $ mx % ⎜⎜⎜ γ2m ⎜⎜⎜ = cos3 ⎜ κ2m l ⎜⎝ ς 2m

⎧ ⎪ γ2m ⎪ ⎪ ⎪ ⎪ ⎪ ⎞ ⎪ ⎪ ⎪ ⎟⎟⎟ ⎪ ⎪ ⎨ κ2m ⎟⎟⎟ where ⎪ ⎪ ⎟⎟⎠ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ς2m ⎪ ⎪ ⎪ ⎩

Dqqqq ⎧ ⎪ δ3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ε ⎪ ⎪ ⎨ 3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ζ3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

(51)

fuuu a¯ 3m + fvvv b¯ 3m + fwww c¯ 3m + 3 fuuv a¯ 2m b¯ m + 3 fuuw a¯ 2m c¯ m + 3 fuvv a¯ m b¯ 2m +6 fuvw c¯ m b¯ m a¯ m + 3 fuww a¯ m c¯ 2m + 3 fvvw b¯ 2m c¯ m + 3 fvww b¯ m c¯ 2m = guuu a¯ 3m + gvvv b¯ 3m + gwww c¯ 3m + 3guuv a¯ 2m b¯ m + 3guuw a¯ 2m c¯ m + 3guvv a¯ m b¯ 2m +6guvw c¯ m b¯ m a¯ m + 3guww a¯ m c¯ 2m + 3gvvw b¯ 2m c¯ m + 3gvww b¯ m c¯ 2m = huuu a¯ 3m + hvvv b¯ 3m + hwww c¯ 3m + 3huuv a¯ 2m b¯ m + 3huuw a¯ 2m c¯ m + 3huvv a¯ m b¯ 2m +6huvw c¯ m b¯ m a¯ m + 3huww a¯ m c¯ 2m + 3hvvw b¯ 2m c¯ m + 3hvww b¯ m c¯ 2m ; =

⎛ $ mx % ⎜⎜⎜ δ3m ⎜⎜⎜ = D(q, q, q, q) = cos ⎜ ε3m l ⎜⎝ ζ 3m 4

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

(52)

where

=

fuuuu a4m + fvvvv b4m + fwwww c4m + 4a3m bm fuuuv + 4a3m cm fuuuw + 6a2m b2m fuuvv +12a2m bm cm fuuvw + 6a2m c2m fuuww + 4am b3m fuvvv + 4am c3m fuwww + 12am b2m cm fuvvw +12am bm c2m fuvww + 4b3m cm fvvvw + 6b2m c2m fvvww + 4bm c3m fvwww = guuuu a4m + gvvvv b4m + gwwww c4m + 4a3m bm guuuv + 4a3m cm guuuw + 6a2m b2m guuvv +12a2m bm cm guuvw + 6a2m c2m guuww + 4am b3m guvvv + 4am c3m guwww + 12am b2m cm guvvw +12am bm c2m guvww + 4b3m cm gvvvw + 6b2m c2m gvvww + 4bm c3m gvwww = huuuu a4m + hvvvv b4m + hwwww c4m + 4a3m bm huuuv + 4a3m cm huuuw + 6a2m b2m huuvv +12a2m bm cm huuvw + 6a2m c2m huuww + 4am b3m huvvv + 4am c3m huwww + 12am b2m cm huvvw +12am bm c2m huvww + 4b3m cm hvvvw + 6b2m c2m hvvww + 4bm c3m hvwww ;

Dqqqq¯

⎛ $ mx % ⎜⎜⎜ ρ3m ⎜⎜⎜ = D(q, q, q, q) ¯ = cos ⎜ ξ3m l ⎜⎝ ϑ 3m 4

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

ρ3m

=

fuuuu a2m |am |2 + fvvvv b2m |bm |2 + fwwww c2m |cm |2       + b3m a¯ m + 3am bm |bm |2 fuvvv + 3a2m |bm |2 + 3b2m |am |2 fuuvv + a3m b¯ m + 3am bm |am |2 fuuuv       + c3m a¯ m + 3am cm |cm |2 fuwww + 3c2m |am |2 + 3a2m |cm |2 fuuww + a3m c¯ m + 3am cm |am |2 fuuuw       + c3m b¯ m + 3bm cm |cm |2 fvwww + 3c2m |bm |2 + 3b2m |cm |2 fvvww + b3m c¯ m + 3bm cm |bm |2 fvvvw   + 3am c2m b¯ m + 3bm c2m a¯ m + 6am bm |cm |2 fuvww   + 6am cm |bm |2 + 3b2m cm a¯ m + 3am b2m c¯ m fuvvw   + 3a2m bm c¯ m + 3a2m cm b¯ m + 6 |am |2 bm cm fuuvw

ξ3m

=

guuuu a2m |am |2 + gvvvv b2m |bm |2 + gwwww c2m |cm |2       + b3m a¯ m + 3am bm |bm |2 guvvv + 3a2m |bm |2 + 3b2m |am |2 guuvv + a3m b¯ m + 3am bm |am |2 guuuv 34

(53)

      + c3m a¯ m + 3am cm |cm |2 guwww + 3c2m |am |2 + 3a2m |cm |2 guuww + a3m c¯ m + 3am cm |am |2 guuuw       + c3m b¯ m + 3bm cm |cm |2 gvwww + 3c2m |bm |2 + 3b2m |cm |2 gvvww + b3m c¯ m + 3bm cm |bm |2 gvvvw   + 3am c2m b¯ m + 3bm c2m a¯ m + 6am bm |cm |2 guvww   + 6am cm |bm |2 + 3b2m cm a¯ m + 3am b2m c¯ m guvvw   + 3a2m bm c¯ m + 3a2m cm b¯ m + 6 |am |2 bm cm guuvw

ϑ3m

= huuuu a2m |am |2 + hvvvv b2m |bm |2 + hwwww c2m |cm |2    + b3m a¯ m + 3am bm |bm |2 huvvv + 3a2m |bm |2 + 3b2m |am |2 huuvv + a3m b¯ m + 3am bm |am |2 huuuv       + c3m a¯ m + 3am cm |cm |2 huwww + 3c2m |am |2 + 3a2m |cm |2 huuww + a3m c¯ m + 3am cm |am |2 huuuw       + c3m b¯ m + 3bm cm |cm |2 hvwww + 3c2m |bm |2 + 3b2m |cm |2 hvvww + b3m c¯ m + 3bm cm |bm |2 hvvvw   + 3am c2m b¯ m + 3bm c2m a¯ m + 6am bm |cm |2 huvww   + 6am cm |bm |2 + 3b2m cm a¯ m + 3am b2m c¯ m huvvw   + 3a2m bm c¯ m + 3a2m cm b¯ m + 6 |am |2 bm cm huuvw

Dqqqq¯

⎛ $ mx % ⎜⎜⎜ σ3m ⎜⎜⎜ = D(q, q, q, ¯ q) ¯ = cos4 ⎜ τ3m l ⎜⎝  3m

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

σ3m

=

fuuuu |am |4 + fvvvv |bm |4 + fwwww |cm |4     + 2bm a¯ m |bm |2 + 2am b¯ m |bm |2 fuvvv + 2cm a¯ m |cm |2 + 2am c¯ m |cm |2 fuwww     + 2bm a¯ m |am |2 + 2am b¯ m |am |2 fuuuv + 4 |am |2 |bm |2 + b2m a¯ 2m + a2m b¯ 2m fuuvv     + 4 |cm |2 |am |2 + c2m a¯ 2m + a2m c¯ 2m fuuww + 4cm c¯ m |bm |2 + b2m c¯ 2m + c2m b¯ 2m fvvww     + 2 |cm |2 cm b¯ m + 2bm c¯ m |cm |2 fvwww + 2 |bm |2 cm b¯ m + 2bm c¯ m |bm |2 fvvvw   + 2 |am |2 cm a¯ m + 2am c¯ m |am |2 fuuuw   + 2a2m b¯ m c¯ m + 4bm c¯ m |am |2 + 4cm b¯ m |am |2 + 2cm bm a¯ 2m fuuvw   + 2b2m a¯ m c¯ m + 4am c¯ m |bm |2 + 4cm a¯ m |bm |2 + 2am cm b¯ 2m fuvvw   + 2c2m a¯ m b¯ m + 4am b¯ m |cm |2 + 4bm a¯ m |cm |2 + 2am bm c¯ 2m fuvww

τ3m

=

guuuu |am |4 + gvvvv |bm |4 + gwwww |cm |4     + 2bm a¯ m |bm |2 + 2am b¯ m |bm |2 guvvv + 2cm a¯ m |cm |2 + 2am c¯ m |cm |2 guwww     + 2bm a¯ m |am |2 + 2am b¯ m |am |2 guuuv + 4 |am |2 |bm |2 + b2m a¯ 2m + a2m b¯ 2m guuvv     + 4 |cm |2 |am |2 + c2m a¯ 2m + a2m c¯ 2m guuww + 4cm c¯ m |bm |2 + b2m c¯ 2m + c2m b¯ 2m gvvww     + 2 |cm |2 cm b¯ m + 2bm c¯ m |cm |2 gvwww + 2 |bm |2 cm b¯ m + 2bm c¯ m |bm |2 gvvvw   + 2 |am |2 cm a¯ m + 2am c¯ m |am |2 guuuw   + 2a2m b¯ m c¯ m + 4bm c¯ m |am |2 + 4cm b¯ m |am |2 + 2cm bm a¯ 2m guuvw   + 2b2m a¯ m c¯ m + 4am c¯ m |bm |2 + 4cm a¯ m |bm |2 + 2am cm b¯ 2m guvvw   + 2c2m a¯ m b¯ m + 4am b¯ m |cm |2 + 4bm a¯ m |cm |2 + 2am bm c¯ 2m guvww 35

(54)

3m

=

4 4 4 huuuu  |cm |    |am | + hvvvv |bm | + hwwww 2 2 + 2bm a¯ m |bm | + 2am b¯ m |bm | huvvv + 2cm a¯ m |cm |2 + 2am c¯ m |cm |2 huwww     + 2bm a¯ m |am |2 + 2am b¯ m |am |2 huuuv + 4 |am |2 |bm |2 + b2m a¯ 2m + a2m b¯ 2m huuvv     + 4 |cm |2 |am |2 + c2m a¯ 2m + a2m c¯ 2m huuww + 4cm c¯ m |bm |2 + b2m c¯ 2m + c2m b¯ 2m hvvww     + 2 |cm |2 cm b¯ m + 2bm c¯ m |cm |2 hvwww + 2 |bm |2 cm b¯ m + 2bm c¯ m |bm |2 hvvvw   + 2 |am |2 cm a¯ m + 2am c¯ m |am |2 huuuw   + 2a2m b¯ m c¯ m + 4bm c¯ m |am |2 + 4cm b¯ m |am |2 + 2cm bm a¯ 2m huuvw   + 2b2m a¯ m c¯ m + 4am c¯ m |bm |2 + 4cm a¯ m |bm |2 + 2am cm b¯ 2m huvvw   + 2c2m a¯ m b¯ m + 4am b¯ m |cm |2 + 4bm a¯ m |cm |2 + 2am bm c¯ 2m huvww ;

Dqq¯ q¯ q¯

γ3m

=

κ3m

=

⎛ $ mx % ⎜⎜⎜ γ3m ⎜⎜⎜ = D(q, q, ¯ q, ¯ q) ¯ = cos ⎜ κ3m l ⎜⎝ ς 3m 4

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

  fuuuu |am |2 a¯ 2m + 2 |bm |2 a¯ m b¯ m + am b¯ 3m + a¯ m b¯ m |bm |2 fuvvv   fvvvv |bm |2 b¯ 2m + 2 |cm |2 a¯ m c¯ m + am c¯ 3m + a¯ m c¯ m |cm |2 fuwww   fwwww |cm |2 c¯ 2m + 2 |am |2 a¯ m b¯ m + bm a¯ 3m + a¯ m b¯ m |am |2 fuuuv     + cm a¯ 3m + 3 |am |2 a¯ m c¯ m fuuuw + bm c¯ 3m + 3 |cm |2 b¯ m c¯ m fvwww     + cm b¯ 3m + 3 |bm |2 b¯ m c¯ m fvvvw + 3 |am |2 b¯ 2m + 3 |bm |2 a¯ 2m fuuvv     + 3 |bm |2 c¯ 2m + 3 |cm |2 b¯ 2m fvvww + 3 |cm |2 a¯ 2m + 3 |am |2 c¯ 2m fuuww   + 3cm b¯ m a¯ 2m + 3bm c¯ m a¯ 2m + 6 |am |2 b¯ m c¯ m fuuvw   + 3am b¯ m c¯ 2m + 3bm a¯ m c¯ 2m + 6¯am b¯ m |cm |2 fuvww   + 3am c¯ m b¯ 2m + 3cm a¯ m b¯ 2m + 6¯am c¯ m |bm |2 fuvvw

  guuuu |am |2 a¯ 2m + 2 |bm |2 a¯ m b¯ m + am b¯ 3m + a¯ m b¯ m |bm |2 guvvv   gvvvv |bm |2 b¯ 2m + 2 |cm |2 a¯ m c¯ m + am c¯ 3m + a¯ m c¯ m |cm |2 guwww   gwwww |cm |2 c¯ 2m + 2 |am |2 a¯ m b¯ m + bm a¯ 3m + a¯ m b¯ m |am |2 guuuv     + cm a¯ 3m + 3 |am |2 a¯ m c¯ m guuuw + bm c¯ 3m + 3 |cm |2 b¯ m c¯ m gvwww     + cm b¯ 3m + 3 |bm |2 b¯ m c¯ m gvvvw + 3 |am |2 b¯ 2m + 3 |bm |2 a¯ 2m guuvv     + 3 |bm |2 c¯ 2m + 3 |cm |2 b¯ 2m gvvww + 3 |cm |2 a¯ 2m + 3 |am |2 c¯ 2m guuww   + 3cm b¯ m a¯ 2m + 3bm c¯ m a¯ 2m + 6 |am |2 b¯ m c¯ m guuvw   + 3am b¯ m c¯ 2m + 3bm a¯ m c¯ 2m + 6¯am b¯ m |cm |2 guvww   + 3am c¯ m b¯ 2m + 3cm a¯ m b¯ 2m + 6¯am c¯ m |bm |2 guvvw

36

(55)

ς3m

=

  huuuu |am |2 a¯ 2m + 2 |bm |2 a¯ m b¯ m + am b¯ 3m + a¯ m b¯ m |bm |2 huvvv   hvvvv |bm |2 b¯ 2m + 2 |cm |2 a¯ m c¯ m + am c¯ 3m + a¯ m c¯ m |cm |2 huwww   hwwww |cm |2 c¯ 2m + 2 |am |2 a¯ m b¯ m + bm a¯ 3m + a¯ m b¯ m |am |2 huuuv     + cm a¯ 3m + 3 |am |2 a¯ m c¯ m huuuw + bm c¯ 3m + 3 |cm |2 b¯ m c¯ m hvwww     + cm b¯ 3m + 3 |bm |2 b¯ m c¯ m hvvvw + 3 |am |2 b¯ 2m + 3 |bm |2 a¯ 2m huuvv     + 3 |bm |2 c¯ 2m + 3 |cm |2 b¯ 2m hvvww + 3 |cm |2 a¯ 2m + 3 |am |2 c¯ 2m huuww   + 3cm b¯ m a¯ 2m + 3bm c¯ m a¯ 2m + 6 |am |2 b¯ m c¯ m huuvw   + 3am b¯ m c¯ 2m + 3bm a¯ m c¯ 2m + 6¯am b¯ m |cm |2 huvww   + 3am c¯ m b¯ 2m + 3cm a¯ m b¯ 2m + 6¯am c¯ m |bm |2 huvvw ; ⎛ $ mx % ⎜⎜⎜ ψ3m ⎜⎜⎜ D(q, ¯ q, ¯ q, ¯ q) ¯ = cos4 ⎜ φ3m l ⎜⎝ θ 3m

⎧ ⎪ ψ3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ φ ⎪ ⎪ ⎨ 3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ θ3m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

= = =

=

=

=

where

fuuuu a¯ 4m + fvvvv b¯ 4m + fwwww c¯ 4m + 4¯a3m b¯ m fuuuv + 4¯a3m cm fuuuw + 6¯a2m b2m fuuvv +12¯a2m b¯ m c¯ m fuuvw + 6¯a2m c¯ 2m fuuww + 4¯am b¯ 3m fuvvv + 4¯am c¯ 3m fuwww + 12¯am b¯ 2m c¯ m fuvvw +12¯am b¯ m c¯ 2m fuvww + 4b¯ 3m c¯ m fvvvw + 6b¯ 2m c¯ 2m fvvww + 4b¯ m c¯ 3m fvwww guuuu a¯ 4m + gvvvv b¯ 4m + gwwww c¯ 4m + 4¯a3m b¯ m guuuv + 4¯a3m cm guuuw + 6¯a2m b2m guuvv +12¯a2m b¯ m c¯ m guuvw + 6¯a2m c¯ 2m guuww + 4¯am b¯ 3m guvvv + 4¯am c¯ 3m guwww + 12¯am b¯ 2m c¯ m guvvw +12¯am b¯ m c¯ 2m guvww + 4b¯ 3m c¯ m gvvvw + 6b¯ 2m c¯ 2m gvvww + 4b¯ m c¯ 3m gvwww huuuu a¯ 4m + hvvvv b¯ 4m + hwwww c¯ 4m + 4¯a3m b¯ m huuuv + 4¯a3m cm huuuw + 6¯a2m b2m huuvv +12¯a2m b¯ m c¯ m huuvw + 6¯a2m c¯ 2m huuww + 4¯am b¯ 3m huvvv + 4¯am c¯ 3m huwww + 12¯am b¯ 2m c¯ m huvvw +12¯am b¯ m c¯ 2m huvww + 4b¯ 3m c¯ m hvvvw + 6b¯ 2m c¯ 2m hvvww + 4b¯ m c¯ 3m hvwww ; ⎛ $ mx % ⎜⎜⎜ δ4m ⎜⎜⎜ E(q, q, q, q, q) = cos5 ⎜ ε4m l ⎜⎝ θ 4m

⎧ ⎪ δ4m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ε4m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ζ4m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

37

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

(57)

where

5 5 5 fuuuuu   am + fvvvvv bm + fwwwww cm 4 4 +5 am bm fuuuuv + am cm fuuuuw + am b4m fuvvvv + am c4m fuwwww + bm c4m fvwwww + b4m cm fvvvvw   +10 a3m b2m fuuuvv + a3m c2m fuuuww + a2m b3m fuuvvv + a2m c3m fuuwww + b3m c2m fvvvww + b2m c3m fvvwww   +20 am b3m cm fuvvvw + am bm c3m fuvvww + a3m bm cm fuuuvw   +30 am b2m c2m fuvvww + a2m b2m cm fuuvvw + a2m bm c2m fuuvww 5 5 5 guuuuu   am + gvvvvv bm + gwwwww cm 4 4 +5 am bm guuuuv + am cm guuuuw + am b4m guvvvv + am c4m guwwww + bm c4m gvwwww + b4m cm gvvvvw   +10 a3m b2m guuuvv + a3m c2m guuuww + a2m b3m guuvvv + a2m c3m guuwww + b3m c2m gvvvww + b2m c3m gvvwww   +20 am b3m cm guvvvw + am bm c3m guvvww + a3m bm cm guuuvw   +30 am b2m c2m guvvww + a2m b2m cm guuvvw + a2m bm c2m guuvww 5 5 5 huuuuu   am + hvvvvv bm + hwwwww cm 4 4 +5 am bm huuuuv + am cm huuuuw + am b4m huvvvv + am c4m huwwww + bm c4m hvwwww + b4m cm hvvvvw   +10 a3m b2m huuuvv + a3m c2m huuuww + a2m b3m huuvvv + a2m c3m huuwww + b3m c2m hvvvww + b2m c3m hvvwww   +20 am b3m cm huvvvw + am bm c3m huvvww + a3m bm cm huuuvw   +30 am b2m c2m huvvww + a2m b2m cm huuvvw + a2m bm c2m huuvww

⎛ $ mx % ⎜⎜⎜ ρ4m ⎜⎜⎜ E(q, q, q, q, q) ¯ = cos5 ⎜ ξ4m l ⎜⎝ ϑ 4m

(56)

where

(58)

ρ4m

=

fuuuuu a3m |am |2 + fvvvvv b3m |bm |2 + fwwwww c3m |cm |2       a4m b¯ m + 4a2m |am |2 bm fuuuuv + a4m c¯ m + 4a2m |am |2 cm fuuuuw + a¯ m b4m + 4am b2m |bm |2 fuvvvv       + a¯ m c4m + 4am |cm |2 c2m fuwwww + b¯ m c4m + 4bm |cm |2 c2m fvwwww + b4m c¯ m + 4cm b2m |bm |2 fvvvvw     + 6 |am |2 b2m am + 4a3m |bm |2 fuuuvv + 6am |am |2 c2m + 4a3m |cm |2 fuuuww     + 6a2m |bm |2 bm + 4 |am |2 b3m fuuvvv + 6a2m |cm |2 cm + 4 |am |2 b3m fuuwww     + 6 |bm | bm c2m + 4b3m |cm |2 fvvvww + 6b2m |cm |2 cm + 4 |bm |2 c3m fvvwww     + 12am bm cm |bm |2 + 4am b3m c¯ m + 4b3m cm a¯ m fuvvvw + 12am bm cm |cm |2 + 4am c3m b¯ m + 4am c3m a¯ m fuvwww     + 12am bm cm |am |2 + 4a3m cm b¯ m + 4a3m bm c¯ m fuuuvw + 12am b2m |cm |2 + 12am c2m |bm |2 + 6b2m c2m a¯ m fuvvww     + 12a2m cm |bm |2 + 12b2m cm |am |2 + 6a2m b2m c¯ m fuuvvw + 12a2m bm |cm |2 + 12bm c2m |am |2 + 6a2m c2m b¯ m fuuvww

ξ4m

=

fuuuuu a3m |am |2 + fvvvvv b3m |bm |2 + fwwwww c3m |cm |2       a4m b¯ m + 4a2m |am |2 bm fuuuuv + a4m c¯ m + 4a2m |am |2 cm fuuuuw + a¯ m b4m + 4am b2m |bm |2 fuvvvv       + a¯ m c4m + 4am |cm |2 c2m fuwwww + b¯ m c4m + 4bm |cm |2 c2m fvwwww + b4m c¯ m + 4cm b2m |bm |2 fvvvvw     + 6 |am |2 b2m am + 4a3m |bm |2 fuuuvv + 6am |am |2 c2m + 4a3m |cm |2 fuuuww     + 6a2m |bm |2 bm + 4 |am |2 b3m fuuvvv + 6a2m |cm |2 cm + 4 |am |2 b3m fuuwww     + 6 |bm | bm c2m + 4b3m |cm |2 fvvvww + 6b2m |cm |2 cm + 4 |bm |2 c3m fvvwww     + 12am bm cm |bm |2 + 4am b3m c¯ m + 4b3m cm a¯ m fuvvvw + 12am bm cm |cm |2 + 4am c3m b¯ m + 4am c3m a¯ m fuvwww     + 12am bm cm |am |2 + 4a3m cm b¯ m + 4a3m bm c¯ m fuuuvw + 12am b2m |cm |2 + 12am c2m |bm |2 + 6b2m c2m a¯ m fuvvww     + 12a2m cm |bm |2 + 12b2m cm |am |2 + 6a2m b2m c¯ m fuuvvw + 12a2m bm |cm |2 + 12bm c2m |am |2 + 6a2m c2m b¯ m fuuvww

ϑ4m

=

2 2 2 3 3 3 fuuuuu am |am | + fvvvvvbm |bm | + fwwwww cm |cm |    a4m b¯ m + 4a2m |am |2 bm fuuuuv + a4m c¯ m + 4a2m |am |2 cm fuuuuw + a¯ m b4m + 4am b2m |bm |2 fuvvvv       + a¯ m c4m + 4am |cm |2 c2m fuwwww + b¯ m c4m + 4bm |cm |2 c2m fvwwww + b4m c¯ m + 4cm b2m |bm |2 fvvvvw     + 6 |am |2 b2m am + 4a3m |bm |2 fuuuvv + 6am |am |2 c2m + 4a3m |cm |2 fuuuww     + 6a2m |bm |2 bm + 4 |am |2 b3m fuuvvv + 6a2m |cm |2 cm + 4 |am |2 b3m fuuwww     + 6 |bm | bm c2m + 4b3m |cm |2 fvvvww + 6b2m |cm |2 cm + 4 |bm |2 c3m fvvwww     + 12am bm cm |bm |2 + 4am b3m c¯ m + 4b3m cm a¯ m fuvvvw + 12am bm cm |cm |2 + 4am c3m b¯ m + 4am c3m a¯ m fuvwww     + 12am bm cm |am |2 + 4a3m cm b¯ m + 4a3m bm c¯ m fuuuvw + 12am b2m |cm |2 + 12am c2m |bm |2 + 6b2m c2m a¯ m fuvvww     + 12a2m cm |bm |2 + 12b2m cm |am |2 + 6a2m b2m c¯ m fuuvvw + 12a2m bm |cm |2 + 12bm c2m |am |2 + 6a2m c2m b¯ m fuuvww (59)

⎛ $ mx % ⎜⎜⎜ σ4m ⎜⎜⎜ E(q, q, q, q, ¯ q) ¯ = cos5 ⎜ τ4m l ⎜⎝  4m

38

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

σ4m

=

4 4 4 fuuuuu    am |am | + fvvvvv bm |bm| + fwwwww  cm |cm |  2 ¯ 4 2 2 + 2am |am | bm + 3 |am | bm fuuuuv + 2am |am |2 c¯ m + 3 |am |4 cm fuuuuw + 2b2m |bm |2 a¯ m + 3 |bm |4 am fuvvvv       + 2c2m |cm |2 a¯ m + 3 |cm |4 am fuwwww + 2c2m |cm |2 b¯ m + 3 |cm |4 bm fvwwww + 2b2m |bm |2 c¯ m + 3 |bm |4 cm fvvvvw     + a3m b¯ 2m + 6am |am |2 |bm |2 + 3¯am |am |2 b2m fuuuvv + a3m c¯ 2m + 6am |am |2 |cm |2 + 3¯am |am |2 c2m fuuuww     + b3m a¯ 2m + 6bm |am |2 |bm |2 + 3b¯ m |bm |2 a2m fuuvvv + c3m a¯ 2m + 6cm |am |2 |cm |2 + 3a2m |cm |2 c¯ m fuuwww     + b3m c¯ 2m + 6bm |bm |2 |cm |2 + 3c2m |bm |2 b¯ m fvvvww + c3m b¯ 2m + 6cm |bm |2 |cm |2 + 3b2m |cm |2 c¯ m fvvwww   + 6am bm |bm |2 c¯ m + 6am |bm |2 cm b¯ m + 6bm |bm |2 cm a¯ m + 2b3m a¯ m c¯ m fuvvvw   + 6am bm |cm |2 c¯ m + 6am |cm |2 cm b¯ m + 6bm |cm |2 cm a¯ m + 2c3m a¯ m b¯ m fuvwww   + 6am cm |am |2 b¯ m + 6am |am |2 bm c¯ m + 6 |am |2 bm cm a¯ m + 2a3m b¯ m c¯ m fuuuvw   + 12am |bm |2 |cm |2 + 6b2m |cm |2 a¯ m + 6c2m |bm |2 a¯ m + 3am c2m b¯ 2m + 3am b2m c¯ 2m fuvvww   + 12cm |am |2 |bm |2 + 6a2m |bm |2 c¯ m + 6b2m |am |2 c¯ m + 3cm a2m b¯ 2m + 3cm b2m a¯ 2m fuuvvw   + 12bm |am |2 |cm |2 + 6c2m |am |2 b¯ m + 6a2m |cm |2 b¯ m + 3bm a2m c¯ 2m + 3bm c2m a¯ 2m fuuvww

τ4m

=

4 4 4 guuuuu    am |am | + gvvvvv bm |bm| + gwwwww  cm |cm |  2 ¯ 4 2 2 2 + 2am |am | bm + 3 |am | bm guuuuv + 2am |am | c¯ m + 3 |am |4 cm guuuuw + 2b2m |bm |2 a¯ m + 3 |bm |4 am guvvvv       + 2c2m |cm |2 a¯ m + 3 |cm |4 am guwwww + 2c2m |cm |2 b¯ m + 3 |cm |4 bm gvwwww + 2b2m |bm |2 c¯ m + 3 |bm |4 cm gvvvvw     + a3m b¯ 2m + 6am |am |2 |bm |2 + 3¯am |am |2 b2m fuuuvv + a3m c¯ 2m + 6am |am |2 |cm |2 + 3¯am |am |2 c2m guuuww     + b3m a¯ 2m + 6bm |am |2 |bm |2 + 3b¯ m |bm |2 a2m fuuvvv + c3m a¯ 2m + 6cm |am |2 |cm |2 + 3a2m |cm |2 c¯ m guuwww     + b3m c¯ 2m + 6bm |bm |2 |cm |2 + 3c2m |bm |2 b¯ m fvvvww + c3m b¯ 2m + 6cm |bm |2 |cm |2 + 3b2m |cm |2 c¯ m gvvwww   + 6am bm |bm |2 c¯ m + 6am |bm |2 cm b¯ m + 6bm |bm |2 cm a¯ m + 2b3m a¯ m c¯ m guvvvw   + 6am bm |cm |2 c¯ m + 6am |cm |2 cm b¯ m + 6bm |cm |2 cm a¯ m + 2c3m a¯ m b¯ m guvwww   + 6am cm |am |2 b¯ m + 6am |am |2 bm c¯ m + 6 |am |2 bm cm a¯ m + 2a3m b¯ m c¯ m guuuvw   + 12am |bm |2 |cm |2 + 6b2m |cm |2 a¯ m + 6c2m |bm |2 a¯ m + 3am c2m b¯ 2m + 3am b2m c¯ 2m guvvww   + 12cm |am |2 |bm |2 + 6a2m |bm |2 c¯ m + 6b2m |am |2 c¯ m + 3cm a2m b¯ 2m + 3cm b2m a¯ 2m guuvvw   + 12bm |am |2 |cm |2 + 6c2m |am |2 b¯ m + 6a2m |cm |2 b¯ m + 3bm a2m c¯ 2m + 3bm c2m a¯ 2m guuvww

4m

=

4 4 4 huuuuu    am |am | + hvvvvv bm |bm| + hwwwww  cm |cm |  2 ¯ 4 2 2 2 + 2am |am | bm + 3 |am | bm huuuuv + 2am |am | c¯ m + 3 |am |4 cm huuuuw + 2b2m |bm |2 a¯ m + 3 |bm |4 am huvvvv       + 2c2m |cm |2 a¯ m + 3 |cm |4 am huwwww + 2c2m |cm |2 b¯ m + 3 |cm |4 bm hvwwww + 2b2m |bm |2 c¯ m + 3 |bm |4 cm hvvvvw     + a3m b¯ 2m + 6am |am |2 |bm |2 + 3¯am |am |2 b2m fuuuvv + a3m c¯ 2m + 6am |am |2 |cm |2 + 3¯am |am |2 c2m huuuww     + b3m a¯ 2m + 6bm |am |2 |bm |2 + 3b¯ m |bm |2 a2m fuuvvv + c3m a¯ 2m + 6cm |am |2 |cm |2 + 3a2m |cm |2 c¯ m huuwww     + b3m c¯ 2m + 6bm |bm |2 |cm |2 + 3c2m |bm |2 b¯ m fvvvww + c3m b¯ 2m + 6cm |bm |2 |cm |2 + 3b2m |cm |2 c¯ m hvvwww   + 6am bm |bm |2 c¯ m + 6am |bm |2 cm b¯ m + 6bm |bm |2 cm a¯ m + 2b3m a¯ m c¯ m huvvvw   + 6am bm |cm |2 c¯ m + 6am |cm |2 cm b¯ m + 6bm |cm |2 cm a¯ m + 2c3m a¯ m b¯ m huvwww   + 6am cm |am |2 b¯ m + 6am |am |2 bm c¯ m + 6 |am |2 bm cm a¯ m + 2a3m b¯ m c¯ m huuuvw   + 12am |bm |2 |cm |2 + 6b2m |cm |2 a¯ m + 6c2m |bm |2 a¯ m + 3am c2m b¯ 2m + 3am b2m c¯ 2m huvvww   + 12cm |am |2 |bm |2 + 6a2m |bm |2 c¯ m + 6b2m |am |2 c¯ m + 3cm a2m b¯ 2m + 3cm b2m a¯ 2m huuvvw   + 12bm |am |2 |cm |2 + 6c2m |am |2 b¯ m + 6a2m |cm |2 b¯ m + 3bm a2m c¯ 2m + 3bm c2m a¯ 2m huuvww (60)

⎛ $ mx % ⎜⎜⎜ γ4m ⎜⎜⎜ E(q, q, q, ¯ q, ¯ q) ¯ = cos5 ⎜ κ4m l ⎜⎝ ς 4m 39

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

γ4m

=

2 2 ¯ 2 2 2 2 fuuuuu   am |am | a¯ m + fvvvvv bm |bm | bm + fwwwww cm |cm | c¯ m 4 ¯ 2 2 4 + 3 |am | bm + 2 |am | a¯ m bm fuuuuv + 3 |am | c¯ m + 2 |am |2 a¯ 2m cm fuuuuw     + 3 |bm |4 a¯ m + 2 |bm |2 b¯ 2m am fuvvvv + 3 |cm |4 a¯ m + 2am |cm |2 c¯ 2m fuwwww     + 3 |cm |4 b¯ m + 2bm |cm |2 c¯ 2m fvwwww + 3 |bm |4 c¯ m + 2cm |bm |2 b¯ 2m fvvvvw   + 3am |am |2 b¯ 3m + 6 |am |2 |bm |2 a¯ m + b2m a¯ 3m fuuuvv   + 3am |am |2 c¯ 2m + 6 |am |2 |cm |2 a¯ m + c2m a¯ 3m fuuuww   + 3bm |bm |2 a¯ 2m + 6 |bm |2 |am |2 b¯ m + a2m b¯ 3m fuuvvv   + 3cm |cm |2 a¯ 2m + 6 |am |2 |cm |2 c¯ m + a2m c¯ 3m fuuwww   + 3bm |bm |2 c¯ 2m + 6 |bm |2 |cm |2 b¯ m + c2m b¯ 3m fvvvww   + 3cm |cm |2 b¯ 2m + 6 |bm |2 |cm |2 c¯ m + b2m c¯ 3m fvvwww   + 6am |bm |2 b¯ m c¯ m + 2am cm b¯ 3m + 6bm |bm |2 a¯ m c¯ m + 6 |bm |2 cm a¯ m b¯ m fuvvvw   + 6am |cm |2 b¯ m c¯ m + 2am bm c¯ 3m + 6cm |cm |2 a¯ m b¯ m + 6 |cm |2 bm a¯ m c¯ m fuvwww   + 6cm |am |2 b¯ m a¯ m + 2bm cm a¯ 3m + 6am |am |2 b¯ m c¯ m + 6 |am |2 bm a¯ m c¯ m fuuuvw   + 12 |bm |2 |cm |2 a¯ m + 6am |bm |2 c¯ 2m + 6am |cm |2 b¯ 2m + 3b2m a¯ m c¯ 2m + 3c2m a¯ m b¯ 2m fuvvww   + 12 |am |2 |bm |2 c¯ m + 6cm |am |2 b¯ 2m + 6cm |bm |2 a¯ 2m + 3a2m c¯ m b¯ 2m + 3b2m c¯ m a¯ 2m fuuvvw   + 12 |am |2 |cm |2 b¯ m + 6bm |am |2 c¯ 2m + 6bm |cm |2 a¯ 2m + 3a2m b¯ m c¯ 2m + 3c2m b¯ m a¯ 2m fuuvww

κ4m

=

2 2 ¯ 2 2 2 2 guuuuu   am |am | a¯ m + gvvvvv bm|bm | bm + gwwwww cm |cm | c¯ m 4 ¯ 2 2 4 + 3 |am | bm + 2 |am | a¯ m bm fuuuuv + 3 |am | c¯ m + 2 |am |2 a¯ 2m cm guuuuw     + 3 |bm |4 a¯ m + 2 |bm |2 b¯ 2m am fuvvvv + 3 |cm |4 a¯ m + 2am |cm |2 c¯ 2m guwwww     + 3 |cm |4 b¯ m + 2bm |cm |2 c¯ 2m fvwwww + 3 |bm |4 c¯ m + 2cm |bm |2 b¯ 2m gvvvvw   + 3am |am |2 b¯ 3m + 6 |am |2 |bm |2 a¯ m + b2m a¯ 3m guuuvv   + 3am |am |2 c¯ 2m + 6 |am |2 |cm |2 a¯ m + c2m a¯ 3m guuuww   + 3bm |bm |2 a¯ 2m + 6 |bm |2 |am |2 b¯ m + a2m b¯ 3m guuvvv   + 3cm |cm |2 a¯ 2m + 6 |am |2 |cm |2 c¯ m + a2m c¯ 3m guuwww   + 3bm |bm |2 c¯ 2m + 6 |bm |2 |cm |2 b¯ m + c2m b¯ 3m gvvvww   + 3cm |cm |2 b¯ 2m + 6 |bm |2 |cm |2 c¯ m + b2m c¯ 3m gvvwww   + 6am |bm |2 b¯ m c¯ m + 2am cm b¯ 3m + 6bm |bm |2 a¯ m c¯ m + 6 |bm |2 cm a¯ m b¯ m guvvvw   + 6am |cm |2 b¯ m c¯ m + 2am bm c¯ 3m + 6cm |cm |2 a¯ m b¯ m + 6 |cm |2 bm a¯ m c¯ m guvwww   + 6cm |am |2 b¯ m a¯ m + 2bm cm a¯ 3m + 6am |am |2 b¯ m c¯ m + 6 |am |2 bm a¯ m c¯ m guuuvw   + 12 |bm |2 |cm |2 a¯ m + 6am |bm |2 c¯ 2m + 6am |cm |2 b¯ 2m + 3b2m a¯ m c¯ 2m + 3c2m a¯ m b¯ 2m guvvww   + 12 |am |2 |bm |2 c¯ m + 6cm |am |2 b¯ 2m + 6cm |bm |2 a¯ 2m + 3a2m c¯ m b¯ 2m + 3b2m c¯ m a¯ 2m guuvvw   + 12 |am |2 |cm |2 b¯ m + 6bm |am |2 c¯ 2m + 6bm |cm |2 a¯ 2m + 3a2m b¯ m c¯ 2m + 3c2m b¯ m a¯ 2m guuvww

40

ς4m

=

2 2 ¯ 2 2 2 2 huuuuu   am |am | a¯ m + hvvvvv bm|bm | bm + hwwwww cm |cm | c¯ m 4 ¯ 2 2 4 + 3 |am | bm + 2 |am | a¯ m bm fuuuuv + 3 |am | c¯ m + 2 |am |2 a¯ 2m cm huuuuw     + 3 |bm |4 a¯ m + 2 |bm |2 b¯ 2m am fuvvvv + 3 |cm |4 a¯ m + 2am |cm |2 c¯ 2m huwwww     + 3 |cm |4 b¯ m + 2bm |cm |2 c¯ 2m fvwwww + 3 |bm |4 c¯ m + 2cm |bm |2 b¯ 2m hvvvvw   + 3am |am |2 b¯ 3m + 6 |am |2 |bm |2 a¯ m + b2m a¯ 3m huuuvv   + 3am |am |2 c¯ 2m + 6 |am |2 |cm |2 a¯ m + c2m a¯ 3m huuuww   + 3bm |bm |2 a¯ 2m + 6 |bm |2 |am |2 b¯ m + a2m b¯ 3m huuvvv   + 3cm |cm |2 a¯ 2m + 6 |am |2 |cm |2 c¯ m + a2m c¯ 3m huuwww   + 3bm |bm |2 c¯ 2m + 6 |bm |2 |cm |2 b¯ m + c2m b¯ 3m hvvvww   + 3cm |cm |2 b¯ 2m + 6 |bm |2 |cm |2 c¯ m + b2m c¯ 3m hvvwww   + 6am |bm |2 b¯ m c¯ m + 2am cm b¯ 3m + 6bm |bm |2 a¯ m c¯ m + 6 |bm |2 cm a¯ m b¯ m huvvvw   + 6am |cm |2 b¯ m c¯ m + 2am bm c¯ 3m + 6cm |cm |2 a¯ m b¯ m + 6 |cm |2 bm a¯ m c¯ m huvwww   + 6cm |am |2 b¯ m a¯ m + 2bm cm a¯ 3m + 6am |am |2 b¯ m c¯ m + 6 |am |2 bm a¯ m c¯ m huuuvw   + 12 |bm |2 |cm |2 a¯ m + 6am |bm |2 c¯ 2m + 6am |cm |2 b¯ 2m + 3b2m a¯ m c¯ 2m + 3c2m a¯ m b¯ 2m huvvww   + 12 |am |2 |bm |2 c¯ m + 6cm |am |2 b¯ 2m + 6cm |bm |2 a¯ 2m + 3a2m c¯ m b¯ 2m + 3b2m c¯ m a¯ 2m huuvvw   + 12 |am |2 |cm |2 b¯ m + 6bm |am |2 c¯ 2m + 6bm |cm |2 a¯ 2m + 3a2m b¯ m c¯ 2m + 3c2m b¯ m a¯ 2m huuvww

⎛ $ mx % ⎜⎜⎜ ψ4m ⎜⎜⎜ E(q, q, ¯ q, ¯ q, ¯ q) ¯ = cos ⎜ φ4m l ⎜⎝ θ 4m 5

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

(61)

where

ψ4m

=

2 3 2 ¯3 2 3 fuuuuu   |am | a¯ m + fvvvvv |bm | bm + fwwwww |cm | c¯ m 2 ¯ 2 2 4 + 4 |am | bm a¯ m + bm a¯ m fuuuuv + 4 |am | c¯ m a¯ 2m + cm a¯ 4m fuuuuw     + 4 |bm |2 a¯ m b¯ 2m + am b¯ 4m fuvvvv + 4 |cm |2 a¯ m c¯ 2m + am c¯ 4m fuwwww     + 4 |cm |2 b¯ m c¯ 2m + bm c¯ 4m fvwwww + 4 |bm |2 c¯ m b¯ 2m + cm b¯ 4m fvvvvw     + 6 |am |2 a¯ m b¯ 2m + 4 |bm |2 a¯ 3m fuuuvv + 6 |am |2 a¯ m c¯ 2m + 4 |cm |2 a¯ 3m fuuuww     + 6 |bm |2 b¯ m a¯ 2m + 4 |am |2 b¯ 3m fuuvvv + 6 |cm |2 c¯ m a¯ 2m + 4 |am |2 c¯ 3m fuuwww     + 6 |bm |2 b¯ m c¯ 2m + 4 |cm |2 b¯ 3m fvvvww + 6 |cm |2 c¯ m b¯ 2m + 4 |bm |2 c¯ 3m fvvwww     + 12 |bm |2 a¯ m b¯ m c¯ m + 4am c¯ m b¯ 3m + 4cm b¯ m c¯ 3m fuvvvw + 12 |cm |2 a¯ m b¯ m c¯ m + 4bm a¯ m c¯ 3m + 4am b¯ m c¯ 3m fuvwww     + 12 |am |2 a¯ m b¯ m c¯ m + 4bm c¯ m a¯ 3m + 4cm b¯ m a¯ 3m fuuuvw + 12 |bm |2 a¯ m c¯ 2m + 12 |cm |2 a¯ m b¯ 2m + 6am b¯ 2m c¯ 2m fuvvww     + 12 |am |2 c¯ m b¯ 2m + 12 |bm |2 c¯ m a¯ 2m + 6cm b¯ 2m a¯ 2m fuuvvw + 12 |am |2 b¯ m c¯ 2m + 12 |cm |2 b¯ m a¯ 2m + 6bm a¯ 2m c¯ 2m fuuvww

φ4m

=

2 3 2 ¯3 2 3 guuuuu   |am | a¯ m + gvvvvv |b m | bm + gwwwww |cm | c¯ m 2 ¯ 2 2 4 + 4 |am | bm a¯ m + bm a¯ m guuuuv + 4 |am | c¯ m a¯ 2m + cm a¯ 4m guuuuw     + 4 |bm |2 a¯ m b¯ 2m + am b¯ 4m guvvvv + 4 |cm |2 a¯ m c¯ 2m + am c¯ 4m guwwww     + 4 |cm |2 b¯ m c¯ 2m + bm c¯ 4m gvwwww + 4 |bm |2 c¯ m b¯ 2m + cm b¯ 4m gvvvvw     + 6 |am |2 a¯ m b¯ 2m + 4 |bm |2 a¯ 3m guuuvv + 6 |am |2 a¯ m c¯ 2m + 4 |cm |2 a¯ 3m guuuww     + 6 |bm |2 b¯ m a¯ 2m + 4 |am |2 b¯ 3m guuvvv + 6 |cm |2 c¯ m a¯ 2m + 4 |am |2 c¯ 3m guuwww     + 6 |bm |2 b¯ m c¯ 2m + 4 |cm |2 b¯ 3m gvvvww + 6 |cm |2 c¯ m b¯ 2m + 4 |bm |2 c¯ 3m gvvwww     + 12 |bm |2 a¯ m b¯ m c¯ m + 4am c¯ m b¯ 3m + 4cm b¯ m c¯ 3m guvvvw + 12 |cm |2 a¯ m b¯ m c¯ m + 4bm a¯ m c¯ 3m + 4am b¯ m c¯ 3m guvwww     + 12 |am |2 a¯ m b¯ m c¯ m + 4bm c¯ m a¯ 3m + 4cm b¯ m a¯ 3m guuuvw + 12 |bm |2 a¯ m c¯ 2m + 12 |cm |2 a¯ m b¯ 2m + 6am b¯ 2m c¯ 2m guvvww     + 12 |am |2 c¯ m b¯ 2m + 12 |bm |2 c¯ m a¯ 2m + 6cm b¯ 2m a¯ 2m guuvvw + 12 |am |2 b¯ m c¯ 2m + 12 |cm |2 b¯ m a¯ 2m + 6bm a¯ 2m c¯ 2m guuvww

41

θ4m

=

2 3 2 ¯3 2 3 huuuuu   |am | a¯ m + hvvvvv |b m | bm + hwwwww |cm | c¯ m 2 ¯ 2 2 4 + 4 |am | bm a¯ m + bm a¯ m huuuuv + 4 |am | c¯ m a¯ 2m + cm a¯ 4m huuuuw     + 4 |bm |2 a¯ m b¯ 2m + am b¯ 4m huvvvv + 4 |cm |2 a¯ m c¯ 2m + am c¯ 4m huwwww     + 4 |cm |2 b¯ m c¯ 2m + bm c¯ 4m hvwwww + 4 |bm |2 c¯ m b¯ 2m + cm b¯ 4m hvvvvw     + 6 |am |2 a¯ m b¯ 2m + 4 |bm |2 a¯ 3m huuuvv + 6 |am |2 a¯ m c¯ 2m + 4 |cm |2 a¯ 3m huuuww     + 6 |bm |2 b¯ m a¯ 2m + 4 |am |2 b¯ 3m huuvvv + 6 |cm |2 c¯ m a¯ 2m + 4 |am |2 c¯ 3m huuwww     + 6 |bm |2 b¯ m c¯ 2m + 4 |cm |2 b¯ 3m hvvvww + 6 |cm |2 c¯ m b¯ 2m + 4 |bm |2 c¯ 3m hvvwww     + 12 |bm |2 a¯ m b¯ m c¯ m + 4am c¯ m b¯ 3m + 4cm b¯ m c¯ 3m huvvvw + 12 |cm |2 a¯ m b¯ m c¯ m + 4bm a¯ m c¯ 3m + 4am b¯ m c¯ 3m huvwww     + 12 |am |2 a¯ m b¯ m c¯ m + 4bm c¯ m a¯ 3m + 4cm b¯ m a¯ 3m huuuvw + 12 |bm |2 a¯ m c¯ 2m + 12 |cm |2 a¯ m b¯ 2m + 6am b¯ 2m c¯ 2m huvvww     + 12 |am |2 c¯ m b¯ 2m + 12 |bm |2 c¯ m a¯ 2m + 6cm b¯ 2m a¯ 2m huuvvw + 12 |am |2 b¯ m c¯ 2m + 12 |cm |2 b¯ m a¯ 2m + 6bm a¯ 2m c¯ 2m huuvww (62)

⎛ $ mx % ⎜⎜⎜ η5m ⎜⎜⎜ E(q, ¯ q, ¯ q, ¯ q, ¯ q) ¯ = cos ⎜ ν5m l ⎜⎝ μ 5m 5

⎧ ⎪ η5m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ν5m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ μ5m ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

⎞ ⎟⎟⎟ ⎟⎟⎟ ⎟⎟⎠

where

5 5 ¯5 fuuuuu   a¯ m + fvvvvv bm + fwwwww c¯ m +5 a¯ 4m b¯ m fuuuuv + a¯ 4m c¯ m fuuuuw + b¯ 4m a¯ m fuvvvv + c¯ 4m a¯ m fuwwww + c¯ 4m b¯ m fvwwww + b¯ 4m c¯ m fvvvvw   +10 a¯ 3m b¯ 2m fuuuvv + a¯ 3m c¯ 2m fuuuww + b¯ 3m a¯ 2m fuuvvv + a¯ 2m c¯ 3m fuuwww + b¯ 3m c¯ 2m fvvvww + b¯ 2m c¯ 3m fvvwww   +20 a¯ m b¯ 3m c¯ m fuvvvw + a¯ m b¯ m c¯ 3m fuvwww + a¯ 3m b¯ m c¯ m fuuuvw   +30 a¯ m b¯ 2m c¯ 2m fuvvww + a¯ 2m b¯ 2m c¯ m fuuvvw + a¯ 2m b¯ m c¯ 2m fuuvww 5 5 ¯5 = guuuuu   a¯ m + gvvvvv bm + gwwwww c¯ m 4 ¯ 4 +5 a¯ m bm guuuuv + a¯ m c¯ m guuuuw + b¯ 4m a¯ m guvvvv + c¯ 4m a¯ m guwwww + c¯ 4m b¯ m gvwwww + b¯ 4m c¯ m gvvvvw   +10 a¯ 3m b¯ 2m guuuvv + a¯ 3m c¯ 2m guuuww + b¯ 3m a¯ 2m guuvvv + a¯ 2m c¯ 3m guuwww + b¯ 3m c¯ 2m gvvvww + b¯ 2m c¯ 3m gvvwww   +20 a¯ m b¯ 3m c¯ m guvvvw + a¯ m b¯ m c¯ 3m guvwww + a¯ 3m b¯ m c¯ m guuuvw   +30 a¯ m b¯ 2m c¯ 2m guvvww + a¯ 2m b¯ 2m c¯ m guuvvw + a¯ 2m b¯ m c¯ 2m guuvww 5 5 ¯5 = huuuuu   a¯ m + hvvvvv bm + hwwwww c¯ m +5 a¯ 4m b¯ m huuuuv + a¯ 4m c¯ m huuuuw + b¯ 4m a¯ m huvvvv + c¯ 4m a¯ m huwwww + c¯ 4m b¯ m hvwwww + b¯ 4m c¯ m hvvvvw   +10 a¯ 3m b¯ 2m huuuvv + a¯ 3m c¯ 2m huuuww + b¯ 3m a¯ 2m huuvvv + a¯ 2m c¯ 3m huuwww + b¯ 3m c¯ 2m hvvvww + b¯ 2m c¯ 3m hvvwww   +20 a¯ m b¯ 3m c¯ m huvvvw + a¯ m b¯ m c¯ 3m huvwww + a¯ 3m b¯ m c¯ m huuuvw   +30 a¯ m b¯ 2m c¯ 2m huvvww + a¯ 2m b¯ 2m c¯ m huuvvw + a¯ 2m b¯ m c¯ 2m huuvww

=

42

(63)