- Email: [email protected]
Degenerate Hopf bifurcation in two dimensions
,Wnnnlrneur Anolyr,.~, Theory, Printed m Great Britam.
Melhods
& Applicmunr,
Vol.
17.
No
3, pp
267.283.
1991 r
DEGENERATE
Budapest
HOPF BIFURCATION
University
of Technology,
Miiegyetem
0362-546X.91 $3,001 .OO 1991 Pergamon Press plc
IN TWO DIMENSIONS
rkp.3,
H-1521
Hungary
and R. E. KOOIJ Delft University
of Technology,
PO Box 356, 2600 AJ Delft, The Netherlands
(Received 13 March 1989; received for publication Key words and phrases: oscillations.
Degenerate
bifurcation,
Hopf
1 February 1991)
bifurcation,
computer
algebra,
nonlinear
INTRODUCTION
CONSIDER the differential
equation
in IR2 with one real parameter i = “0,
x, Y) (1)
j, = g(p,x,A
1’
wheref(p, 0, 0) = g(p, 0, 0) = 0, ,D E IR. Suppose at a critical parameter value, p = p,,, a Hopf bifurcation occurs. For the sake of simplicity let p, = 0 and let the coefficient matrix of the linear part of (1) at ,D = 0 be
[
0w 0#0. --o 0I;
The type of this bifurcation (super- or subcritical) is determined by the stability ,D = 0. If the first Poincare-Lyapunov constant g, is negative, then the origin and if g, > 0, then it is a repellor. The formula for g, is well known:
of the origin at is an attractor,
where the partial derivatives are to be taken at (p, x, y) = (0,0,O). If g, = 0, then a degeneracy occurs and the sign of the first nonvanishing PoincareLyapunov constant g, determines the stability of the origin. The greater the number k, the more drastically would increase the complexity of the formula determining g, as a function of the appropriate partial derivatives at (0, 0,O). We introduce here an algorithm for the determination of g, for any k (Section I). We use this algorithm in an algebraic computer program written in Macsyma 420 language and accomplished on a Symbolics 3620 computer. The computation times achieved are surprisingly short. The algorithm is rather simple and uses only elementary vector and matrix operations. However, the exact proof (Section 3) needs quite new mathematical apparatus (Section 2). This apparatus works with multilinear maps. Here we introduce these new tools only for the 267
268
V. KERTI!SI
R. E. KOOIJ
and
two-dimensional case where multilinear maps are represented by coordinates. These coordinate representations will be called arrays. Mathematical results for the n-dimensional case can be found in [l]. It seems the apparatus works in the infinite dimensional case as well. In the case of equation (I), the degenerate Hopf bifurcation is not typical and may not be of much interest. However, if (1) consists of more than one, say k, parameters then in the kdimensional parameter space there exists typically a point where the first k-l PoincareLyapunov constants vanish. This point, as an organizing center, determines the local dynamical behavior of the differential equation in its neighborhood. This behavior is typical and indispensable to revealing the global dynamical behavior. This is the reason that great efforts have been made to obtain an effective way for determining the type of the Hopf bifurcation [2-121. To find the organizing center one need the parametric handle of the algorithm used. For this purpose our algorithm has been incorporated into the algebraic programming system mentioned above. In Section 4 we show some simple examples with and without parameters. In Section 5, as an example of practical importance, a mathematical model with three parameters will be studied. It has already been discussed in several papers [ 13- 151; still its dynamical behavior is so complex that it has not yet succeeded in revealing the entire phase portrait. In this section the nondegenerate and degenerate Hopf bifurcation of this model will be thoroughly investigated. I. ALGORITHhl
Consider
FOR
DETI~RR1IN,4TiON
the two-dimensional
OF
differential
THE
POINC‘ARE-I
Y.4PUNOV
CONSTANTS
equation
k = f‘(s, y) (2) _ii = g(.u, y) 1 ’ where
f and
g
are
smooth
functions
and
f(0,0) = g(0,0) = f,(O, 0) = gv(O, 0) = 0,
f,.(O,O) = -g,(O,O)= 0 # 0. Let A be an n x 111 matrix with elements u!,. We shall use the notation A(i,j) = CI,,. Similarly, if cl E R” and c7 = (~1,, ~1~). . , [I,,), then let v(j) = L;,. T denotes the transpose. k 0n
1s the binomial
coefficient
defined
as
k! (k -
0 THEOREM
1
bk+, 9 L’h+I
ifO
(algorithm). Define the vectors and matrices c’?, e,,,,, 12A+3, Ezxtr, and dh _, ; k = 2, 3, . .; I = 1, 2, . . . ; k + I 2 4 as follows:
U2’= (l,O,
1); 1
ifi=
1
Ck(i) =
e, E R!“;
i0 I, is the
n)‘n’ . .
k
x
k identity
otherwise; matrix;
CZht2, A,,,,
269
Degenerate Hopf bifurcation in two dimensions
E,, is a (2k + 1) x (2k) matrix,
&(i,j)
=
1
ifi=j
0
otherwise;
1 CZk is a (2k) x (2k + 1) matrix,
G(i,j)
ifi=j ifj=2k+
=
landi=2/+
1
otherwise;
A,,, is a (k + I) x (k + 1) matrix; k +
a/+j-i- I
x
a;+j-i-l
1
_
ai+l-j $+1-J
g
j
al+j-ia’-.if
(0, 0);
bk, vk E IRk+‘; k-l bk
=
c /=2
h,k+l-/U/i
if k is odd
--/+,bk uk = i
’
if k is even;
-~&(C,A,,,&-lC,bl
where the inverses
always exist, and
d, E R2k+3; dk’ = The Poincare-Lyapunov
elk+3[z2k+3
constants
Note. The vectors dk are constants
-
A2k+2.,E2k+2(C2k+2A2k+2,,E2k+2)-1C2k+21.
of equations
(1) are gk = dTb2k+2.
not depending
d,r = $(3, 0, 1, 0, 3), 2. INDICES
on f and g. For example:
d,T = &(5, 0, 1, 0, 1, 0, 5). AND ARRAYS
and IK'= (1,2). Let lKnC1= IK"x IK', n E N. If i E IK" then let Let IN = {1,2,3,...) i = (i,,i, , . . . , i,), where ik E IK',1 5 k 5 n. We call i an index from IK"and i, is called the kth element if i. If i E IK"and j E IK'" then let (i, j) = (i,,i,,...,i,,j,,j,, ...,j,,)E lK"+"'. For convenience we use IK"as the empty set, and for n and m above it is allowed to equal zero. If n = 0 then (i, j) = (j, i) = j. If n = m = 0 then (i, j) is an “element” of IK'.Let ,$J= UnEN IK". 3 will be called the index set. Now we define some functions.
Definition
1. p: 3 x N -+ NJ; p(i, j) = m if exactly m elements
Definirion 2. 0’: 3 + IN; a(i) = p(i, 1) + p(i, 2).
of i have the value j.
270
V. KERT~SZ and R. E. Koou
Definition 3. y: 3 --t 2”; y(i) = {j E IK”(‘)1v k E IK’ p(i, k) = p(j, k)]. Let H be an arbitrary finite set. Then let N(H) denote the number of the elements
of H.
Definition 4. cp: 3 + N; q(i) = N(y(i)). The following relations are easy to prove:
Co(i) (o()
_ P(i, t) a(i)
+ 1 + 1 ’
where t E IK’.
Let R be the following equivalence relation on 3: i, j E 3, i = Rj if and only if y(i) = y(j). Let 3 j R be the set of equivalence classes defined by R. (We can define IK” 1R in the same way.) Obviously, if i = Rj then p(i, k) = p(j, k), k E IK’; a(i) = a(j); p(i) = r&j). Consequently, 0, cp and p(. , k) can be considered as functions defined on 3 1R. Obviously, iV( IK” 1R) = rn + 1. Definition 5. v: $1 R + iN; v(i) = p(i, 2) + 1. For any fixed k E IN, the restriction of v to IK” 1R is a bijection, this set. Let 6(k, .) be the inverse of this bijection.
namely
an enumeration
of
Definition 6. 79: N x N + JJ I R; 29(k, i) is the unique element of IKk I R for which the relation v(6(k, i)) = i holds, i 5 k + 1. We can derive an explicit formula for 6: 6(k, m) = (i, j), where the indices i, j E 3 I R are uniquely determined by the relations o(i) = p(i, 1) = k + 1 ~ m and cr(j) = p(j, 2) = m - 1. Functions of type A : IK” + R will be called arrays. The integer n is the type of array. The set of arrays of type n will be denoted by a”. Obviously, Q2 is the set of 2 x 2 matrices, and 0.’ is R2. Let 6!” be R. If an element of i E IK” is fixed then A E 6!” as a function of the remaining elements is obviously an element of Q”-‘. If i E IK” and i = (i, , i,, , i,,), then instead of A(i) if i E IK”] and or A((;, , iz, . . ., i,!)) we shall write simply /I(;, , i,, . . . . i,,), as well. Similarly, j E IK”1, where n, + n2 = n, then instead of A((i,j)) we shall write A(i,j). Note that if n = 1 or n = 2, then these notations correspond to vectors or matrices, respectively. Q” as a set of functions with real values is a linear space using the necessary usual way: (A + B)(i) = A(i) + B(i), etc. Definition 7. The array A E (I?“, n 2 2 is called (i) symmetrical, if A(i) = A(j), wherej E v(i), i E IK”; (ii) quasisymmetrical, if either n = 2 or n > 2 and fixing i, , A as a function is a symmetrical element of a”-‘; (iii) antisymmetrical, if xJt7(,) A(j) = 0, i E K”.
definitions
in the
of i,, i,, . . . , i,,
The set of symmetrical, quasisymmetrical and antisymmetrical elements of Q” will be denoted by C”, Q” and Z”, respectively. These are obviously linear subspaces of Q”. If S E [” then it can be considered as a function of type IK” / R - I?, as well. THEOREM
2. Q” = i” @ 2”.
Degenerate
Proof.
Hopf
bifurcation
Let A E Q” and S and 2 be defined S(i) = $)
c
as follows: Z(i)
A(
271
in two dimensions
= A(i)
i E IK".
- S(i);
JET(i)
Then A = S + Z and S is symmetrical, Z is antisymmetrical. We now show that this decomposition is unique. Suppose A = S, + Z,, where S, E [” and Z, E 2”. Then S - S, = Z, - Z. Let Sz = S - S, . S2 is symmetrical and antisymmetrical at the same time. Consequently C S(j) = &i)&(i) = 0. j E y(i); S,(i) = S,(j); jEy(r) Considering
that p(i) # 0, it follows
Let the symmetrical
component
n
S, = 0.
of A E Q.’ be denoted
by $&4).
8. Let A E a” and B E a.“‘, nm L 1. Then AB E @.n+mm2and
Definition
AB(i,j)
=
c
i E IKn-',
A(i, kVW,_i),
j
E
IK"-'.
k=l
Distributive
and associative
A(B provided
laws can be easily verified: (B + C)A
+ C) = AB + AC,
= BA + CA,
(AB)C
= A(BC),
in the last case B is not a vector.
Consider the special case where A E a”, n > 2, uk E R2, 1 5 k 5 p 5 n. The multiple product (. . .((Au,)u,)u, . . .)u,, will be denoted simply by Au, u2, . . ., up. If U, = u2 = ... = up = U, then we write Au”. Au” is a homogeneous polynomial of the elements of U. The degree of this polynomial is n. More preciously: Au” =
c
A(i)u~“‘r)U~(‘,‘)
;clK"
=
c
rp(i)S(A)(i)u~“.“u,P”.”
= $JA)u”;
uT = (u, ,
;tlK"IR
Zu”=Oforeveryu~/R~ifandonlyifZ~Z”.LetS~~”andu,~IR~,k= Then Su,, uk,, . . . , uk, = A%,4, where k, , k,, . . . , kp is any permutation
1,2,...,p
of 1,2, . . . ,p.
Let S E [” and let u be a differentiable
function
$ Su”(t)
of type R -+ I?‘. Then
= nSu”-‘(t&(t).
PROPOSITION 3. Let S E c”, B E ak; n, k > 1 and u E R2. Then
(Sun-‘)(Buk-‘) Proof.
(SU”-‘)(Buk-‘)
= Suu, . . . , u(Buk-‘)
= (SB)U”+~-~. = S(Buk-‘)u,
Now we define an Q” 4 R-type linear operator.
. . _, u = (SB)U”+~-~.
u2).
V.
272
Definition
R. E. KOOIJ
KERT~SZ and
9. Let YE O”, n 2 2. Then
v E I?“+’ and
Vi) = d6(n, i)E( V(W, Some basic properties
of this operator
are: ,,+ I Vu” = C V(i)u, p(run.iJ.1)u2/Hij(fl,l),2).
v;
giq=
0).
i=,
the map s(V)
+ v is a bijection; Vu” is identically zero if and only if v = 0, where V E (3”. Let now V E azk and uT = (u,, u,). VuZk = V(1, 1, . . , l)(uf + u:)~ for every u E I?* if and only if C,,v = 0. (For C,, see Section 1.) Definition (k + n -
10. Let 2 5 k be a given integer. Then for any determined by the equation
A E Q”, n 1 2, X,(A)
1) x (k + 1) matrix
where Wi E i”;
W,(I) =
PROPOSITION 4. Let A E (IV, n 2 2. Then
Proof.
Using the definition
1/do
if v(/) = j
0
otherwise.
for every
V E ck, k 2 2, VA = 3Kk(A)I/.
of W,, we have k+l v
=
c ,=I
w,w,,
which has the consequence h+I
VA(i) The equality
ktl
= 1 &j)WjA(i) ,= I
of the first and third expression
PROPOSITION 5. Let A, E
= c mk(A)(i,j)@j). ,= I is equivalent
n
to the statement.
Q” ‘, then
where p(s, , 1) = I + j - i, p(s, , 2) = i ~ j, p(s,, I)=/+,;-i-
l,p(s,,2)=i+
1 -j.
Proof. ntk(A,)(i,j)
= (W,A,)(i) _
= cp(ti(k + I ~ 1, i)S(W,A,)(ti(k (W/$)(p)
c
pt-p(ii(h+/-1.1)
_
cp(r,Ms,) v((r,,
1))
A,(l,s,)
=
c
'I,, ty(rY(k+l-1.1)
cP(r2)ds2) + (D((rL) 2)) A,@, sz),
+ I -
1, i))
[W,(r, l)A,(l,s)
+ y(r,
2M,(2,s)l
is a
Degenerate
Hopf
bifurcation
in two dimensions
273
where
(r,, 1) E ~(W,j)h (r2, 2)
dW,j)).
E
Note that k+l-j Ar,, 1) + 1 &r,) cp((r, , 1)) = a@,) + 1 = k Ar2, 2) + 1 _jdrd v)((r2, 2)) = a(r2) + 1 As,, 1) = &VI,
1. ’
k
s,>, 1) - Ar,,
1) = j + I - i;
As,, 1) = A(r2, s,>, 1) - Ar,, O(S,) = c&)
P(Sl)=
(p;;fl)) =(if j);
ds2)
(,~yJ
=
(2) can be rewritten
=
OF THE
+
Ak(l, i) =
n
J.
ALGORITHM
off
+ Al,22 + ... + A,d
and g with an appropriate
r E R\J:
+ o(lulr+‘),
pci, l)+P(i,2)
f a
k’ap(‘.1)xap”,2’y
a positive
definite u(u) = 1,and
whereVkEck;2
0);
polynomial
Ak
i)
=
F
1 (pci, l)+P(1,2) (0, 0). pci, ljx ap(i,2)y g
function
u: R2 + R of the form
v2u2+ v,u3 + ... + v,+,Ur+l, 0
[ 0
v,= the derivatives
t2,
Lyapunov
1
Compute
:
A, E Qk+‘, 1 5 k 5 r, 1
Now we define
(i
using the expansion
ic = A,u where 1.4’ = (x,y),
1) = j + I - i - 1;
= I;
3. PROOF
Equation
’
of u along the solutions
of (2’):
r+l + c (kVkA,
ti(U) = (2&4,)U2
k=3
where k-l Bk
=
c /=2
1
1’
II/IAk+l-,.
+ &)Uk
+ o(IU1r+2),
(2’)
274
V. KERT~SZand R. E. Koorr
Compute
312, (A ,). Obviously
-~ i k f t%,(A,)(i,j)
=
k +2 - i k \O
ifj=
i-
1
otherwise.
If k is odd then let (kV,A, + B,)uk = 0. Or, equivalently: knt,(A,)V’ + Bk = 0. Using the Gaussian elimination, we can easily see that Fm,(A,) is regular. Consequently: i;;, = -(l/k)X,(A,)-‘B,. If k is even then let (k&A, + B,)uk = c(x2 + y2)k’2, c E IF?.Or, equivalently: C,(kZVZ, (A ,)V, + Bk) = 0. Let the (k + 1)th element of Vk be zero. Then Ek ErV;. . Note that C, 3nk (A ,)E, is regular (use the Gaussian elimination). Now we can express _ v, : vk = -~El(Cx~,(A,)~~)-‘C,17,. Lastly,
note that k-l Bk
consider theorem
=
1 I=2
proposition 5 and use the notations 1 is completed. n
3. COMPUTATION
WVk+,-,I&;
u& = vk, bk = Bk, A, , = %,(A,).
BY
COMPUTER
The proof
of
ALGEBRA
We developed a computer program (on a Symbolics 3620 computer written in Macsyma 420 language) which carries out the algorithm explained in Section 1. The program computes the Lyapunov function and the Poincare-Lyapunov constants. We introduce here some computational examples. The simple examples 1 and 2 were solved by hand in order to show how the algorithm works. Examples 3 and 4 are more complex and were carried out with our computer program. A rather sophisticated example is shown in the next section.
Example 1. .t= -y+x2 j, The partial
derivatives
=x+x2.
in the origin are: f,,(O, 0) = g,(O,
0) = 2. Any other derivative
is zero.
Degenerate
A 3.1 =
s-i;
Hopf
bifurcation
_;;];
in two dimensions
-%t=i[i;
275
;;I;
f
I
uj = -fA,‘,
01
A 3,2
c
1
-2
6, =
I1 -41I 0
b4 = 20, + 3A,,,q
;
03 = 0;
I-
+ 50
$00
0
00’
-4
oot1
=
’
b, =
-4
;
g, = +(3
0
1
0
3)b, = -+.
0 o_Ii Check the result by computing
u4
10
-4 A 4.1
=
0
and derivating
0I the Lyapunov
2
00 00 -+E4(C4A4,,E4)-‘Cd
$
o-2
function.
-+E,(C,A,,,E,)-‘Cab,
=
0
E4 =
[
04
= 4
! 1 -3 -50
02
-101
042
053
;
216
+ I1
V. KEKTESZ and R. E. KOOIJ
X3
x3
-4)
VJ =
X2Y
(-1
-*
x’y
0 + 0)
XY2
x2y2
X_$
YS
Example
2,
0-
-1 6, = 2a,
v3 =
LIZ
1 -I_
6, = 3A,,2v3 :=
-2-
00
3 a, =
0
0
1
0
0
6
0
2
0
I 3 0
3
3
0
3 -2
0
2
0
0
0
0
1
-3 li
2
-?
= 1
0
0
0
0
1
YJ
= I!‘I 0
2 3
6
0
2
-3
o_ 0
-2
-2
_z 3
0 vq=
0
.
0 !I 0
Note that if bk = 0 then L’~= 0 and b,,, = kA,,, vk, k = 4,5, . . . . We can conclude vk = 0 if k 2 4 and g, = 0 for every k. Consequently, v = x2 i y2 + fx’ - +yv’ - 2xy2 is the first integral
of our equation.
Actually:
ti = 0.
Example 3. i=
-_v+x”
j = .Y g, = g, = g, = g, = 0;
231 g5 = m
7337 793 U = x2 + y2 - 512x1’_v - 1536xq’3
9273 - mx’y5
231 7623 1309 - ~ x5y7 ___ x3y9 - 512xy1’. 1280 512
Degenerate
Hopf
bifurcation
271
in two dimensions
Example 4 (see [2]). i = -y - ax2 + 2xy + by2 j = x + x2 + (5b - 3a)xy - y2 g, = +$(b - a)3(2b2 - ab + 1)
g, = g, = 0; 5. INVESTIGATION
Consider
the differential
OF A REACTION-KINETIC
MODEL
equation T = -/3T + BDaCe’, C,T>O
c = 1 - C - DaCer;
1 ’
(3)
where p, B and Da are positive constants (parameters). This equation represents the first-order exothermic chemical reaction A + B in a continuously stirred tank reactor and was investigated in great detail in [ 131 and [ 141. The same equation is used in [ 151 as a simplified model for a coal char particle combustion. Equation (3) possesses an extremely complex dynamical behavior which has not yet been entirely revealed. In order to be able to reveal the entire global dynamic behavior one needs the complete local investigation. Amongst others, local investigation contains the Hopf bifurcation study. Degenerate cases are not typical, still, studying the degenerate Hopf bifurcation, we can reveal typical bifurcation phenomena around the degenerate point as the organizing center. Let us consider a system with two real parameters:
k = f(Q, P, x, Y) i
=
1
(4)
g(& P, x, Y).
Suppose there is a function @: R + R such that if O( = @(/I) then a Hopf bifurcation occurs in (4). Assume moreover there exists a & such that if /I < /IO then the bifurcation is supercritical (g, < 0) and if /3 > /3(, then subcritical (g, > 0). Assume the periodic orbits appear in the supercritical case if 01 > Q(p) and in the subcritical case if (Y < 0. On this diagram the origin is the organizing center. Now we carry out a total Hopf bifurcation investigation on equation (3). First of all we introduce a new parameter instead of B: p = /i’/BDa. In this way we can determine that the location of the equilibra namely on ,u and Da. Let 8 and y the coordinates of an equilibrium: 0 = -/I0 0 = 1 -y
+ f ye” -
WDae”.
depends
only on two parameters,
V. KERT~Z and R. E. Koorr
278
At1
A/i = const.
t
UE: unstable ULC: unstable
SE: stable
equilibrium limit cycle
SLC: stable
Fig.
equilibrium limit
cycle
I.
Note that 0 < w < 1. We cannot express 0 and I,Uexplicitly as functions of ,Qand Da. Therefore we use ,8, 0 and IJ/ as new parameters, and we introduce x and y as new variables: x=T-H y=c-rp
By these, equation
(3) yields:
.k = P(H - 1)X +
‘fy +p&(e’
-
1 - x) + iy(e’
-
1))
(5) j, = (cc/ - 1)x - ty
+ (w - I)((e’ - I - x) + iy(e’
-
1)).
1
Degenerate
The coefficient
matrix
Hopf
of the linearized
bifurcation
279
in two dimensions
system is
A=
The necessary
conditions
of the Hopf bifurcation
1<8<1 w Substituting
are: det A > 0 and tr A = 0, that is p=
and
l (0 - L)cy’
the critical /I value in (5) we obtain 1 8 i = w” + (0 _ 1)w2Y
j = (tf - 1)x - :Y In order to obtain
the canonical
I9
1 - X) + +_v(e” - 1 - x)
(e" -
+ (0 -
lb
1
I
+ (w - 1) (e” - 1 - x) + t (e’ - 1) . I 1 form,
we transform
the variables:
As a result: if=ov+g (6)
)
i, = -_ou ii cog 1 where 1 + wQ2 g = w exp (1 _ vI)v2 (
’ +~2v~u) (1 - w)w
u-l-
+ (cou -
iU)(exp(:+_$G2z4
-
O
The next step is to determine w and w: g, = (w
2c02 + 1)2(l/o4
the Poincare-Lyapunov
+ 24U3w2 + l&02
- I+2
constants
1)
0 < 0.
g,,g, and g, as functions
of
- 2l/Y3 + 74!/2 - 5l+U+ 1)
S(I,V - 1)4y/6w2 g, = -(v/‘w2
+ 1)(186~8~‘o
1967y/‘o”
- 3051,~~o~ + 13951,~‘o~ - 145w6w8 - 1861,~~~~ + 1621,~~~~
+ 49561~1’0~ - 1244~‘~~
- 1098~~0~
- 6289~~~11~ + 12843~~~~
- 7768~~0~
+ 2808~~~~
+ 201661~1~0~ - 185041//‘w’
-
10819~5w2
+ 1860~ - 1201,~’ + 540~~ - 750~’
•t 345~‘~~
+ 3671y3w6 + 990~‘~~
+ 465~’
-
f 899w2w4 - 186tj~0” - 276~‘~~ + 8718~~~~
135~ + 15)/(1152(1~
- 2033~0’ -
1)8~‘206)
280
V. KERT~
g, = (‘//‘w’ + 1)“(1341751//‘2w’6 - 134175~9w’4
-
- 27960831t~/‘~w’~
- 89543353~“~’ + 8876630~swx + 66376480~‘~~ + 106174634~sw6 134175ww6 -
153360~“~’
- 474652~‘~~‘~
+ 6389651j/‘~w’~
- 3855099~‘~‘~ 6082518’~/“w’~
+ 429615O@w” 264523371y”w’
and R. E. KOOII
- 5378394~“~‘~
1219977~~~‘~
+ 46788711’/1’w”’ - 753500~~~‘~
-
+ 522950~~0~
185245120~8w”
-
-
+ 1032312’~1’w~ - 4620172~‘~~ + 7022916’~/~w’
+ 88800~’ 36075~
- 488400@
2775)/(442368(~
+ 932784~“~~
1375224~‘~~
18746700~‘~~
+ 9356650~~~~
- 2949788~‘~’
+ 976800~’
-
12385480’/“w6
- 237643272~~~~ + 1173485~~~~
-
+ 74127618’j/‘w4
- 268137446~~~~
+ 162416891y2w4 - 2246523~~~
+
+ 15359604~6ws
+ 283247706~‘~~
- 20262648~~~~
+
- 2963550~‘~‘~
- 3266418’//“w’
+ 13828311 l’,//‘w’ - 90400342~‘~~
- 180498774’+~‘w” + 276233900’/‘0
+ 36075~ 7308~’
- 3948321//“w’~
- 22318196~~~‘~
+ 88296~‘~~~
- 4255655~“~’
- 341149’$“w’4
+ 11365519~‘“w’2
+ 5026841~‘~‘~
+ 26360401y’“w4
- 66201688~‘~”
+ 1762334~“~‘~
+ 166747467~~~~
+ 1341750”
- 74544~~~~
- 10434862~~0~
+ 763343~~~~
-
197025~~
- 112505’//w2 +
- 999000’/1~ + 582750~”
-
197025ry’ +
- l)‘21y’xw”‘).
Note. Equation (6) contains two parameters, namely I,Yand w. In the general case of two parameters the typical situation is the following. On the parameter plane there is a curve which separates the domains where g, > 0 and g, < 0, respectively. Along this curve g, = 0. On the curve there is a point where g, = 0 and which separates the points where g, > 0 and g, < 0, respectively. We need g, only in this point. In case of equation (6) this typical situation does not occur as the point g, = 0 is outside of the possible domain of the parameters. We have shown the expression of g, above only in order to demonstrate the possibilities of our algebraic program. Continuing the computation we solve the equation g, = 0. As a result we obtain on the parameter plane (see Fig. 2) which have the equations y, : w = [l ~ ‘// ~ 2~’
+ (1 ~ ~)(81y
- 8~’ + l)“2/(2y/3)]“2
y2: w = [l - y/ - 2’~/~ - (1 - w)(8w3 - 8~’ + 1)“2/(2w3)]“2. The sign of g, outside
yl and y2 can be easily determined
(see the figure).
two curves
(7) (8)
Degenerate
Hopf
bifurcation
in two dimensions
Fig. 2.
b)
Da p1 = e”‘/H,* Da,=e”‘!(H,H, =3+
L-j 2
Fig. 3.
1)
281
282
V. KERT~SZ and R. E. KOOIJ
Now we determine the sign of g, along y, and yz . First we leave out the factors with positive sign from g, and then we substitute the expressions (7) and (8) for w. So we obtain g,, and gzz, respectively. g,, = (-64~~’ + 1282~”
+(8w3
- 81y2 + 1)“2(24v”
- 2496~~
+ 1736~~ - 3801,~~ -
- 35241,~” + 9240~” - 9001~1~- 98~’
- 228~”
- 14002~8
+ 111~ -
g,, = -(641c/12 + (8~~ - Q2
+ 8321,~’ - 1402’,~’ + 7081,~’
1704~~ + 112’~ - 18) + 72&//l’
+ 116601y’ - 3302~~
- 2840~’
+ 3008~~
18)/@
+ 1)“2(24w”
- 2281,~“’ + 832~~ -
1402~~ + 708~’
+ 1282ry’ - 24961,~~ + 17361~’ - 380~~ ~ 17Oy/’ + 112~ - 18) - 728~” + 3524~“’
- 9240~’
+ 900’//3 + 98I$ A further
computation
-
+ 14002’//’ - 1166Oy’
+ 3302~~
+ 2840~~
- 3008~”
112l/Y + 18)/‘J/h.
yields the result that g, < 0 on yl and g, > 0 on y2.
As a summary, consider Fig. 3(a) where the coordinate B of the equilibrium is plotted vs the parameters Da and ,u. The Hopf bifurcation diagram is shown on the @Da, ,u) surface. The projection of this diagram on to the Da, p plane can be seen in Fig. 3(b). At the critical value PC,,, = ((0 - l)cy)-’ of the third parameter p, a super- or subcritical Hopf bifurcation occurs according to the figure. On the curve g, = 0 the Hopf bifurcation is degenerate. The sign of g, determines the type of this degenerate Hopf bifurcation (see Fig. 2 as an illustration).
REFERENCES 1. KERT~SZ V., hlultilinear maps and bifurcation theory. To appear in Co//. math. Sot. Jdnos Solyur, Szeged (Hungary). 2. BAUTIN N. N., On the number of limit cycles which appear with the variation of coefficients from an equilibrium position of focus or center type, Mal. Sh. 30, (72) 181-196 (1952); Am. marh. Sot. Trunsl. No. 100 (1954). (In Russian.) 3. ANDRO~.OW A. A., LEONTOVI~~IE. A., GVRDON I. I. & MAIEK A. G., Qua/r/u/we Theory of Second Order D_vnarnicul S_vslems. Halsted Press, New York (1973). 4. MARSDEN J. E. & RICCRAC‘KEN bl., T/?eHopj’Bi~urcution und irs App/icotiom. Springer, New York (1976). 5. HASSARD B. & WAN Y. H.. Bifurcation formulae derived from center manifold theory, J. math. Ana/.vsis Applic. 63, 297-312 (1978). 6. SCHMIDT D. S., Hopf’s bifurcation theorem and the center theorem of Lyapunov with resonance cases, J. tnafh. .4na/.vsis Apphc. 63, 354-370 (I 978). 7. GIBBER F. 8; WILLAMOWSKI K. D., Lyapunov approach to multiple Hopf bifurcation, J. mar/z. Analysis Applrc. 71, 333-350 (1979). 8. NEGRINI P. & SALVADORI L.. Attractivity and Hopf bifurcation, Nonlinear Anu/_ysi;, 3, 87-99 (1979). 9. GOLLBITSKY hl. & LANWORD W. F., Clarsification and unfolding of degenerate Hopf bifurcation, J. d/yf. Eqns41, 357-415 (1981). 10. SHI SONGLIXG, A method of constructing cycles without contact around a weak focus, J. d/j__ Eyns 41, 301-312 (1981). 11. HASSERD B. D., KAZARINO~~ N. D. & \VAN Y-H., Theorj~andApplicaiion of HogfBifirrcarion. Cambridge Univertity Pres\, Cambridge (1981).
Degenerate
Hopf
bifurcation
in two dimensions
283
12. SHI SONGLING, On the structure of Poincav-Lyapunov constants for weak focus of Polynomial vector fields, J. diff. Eqns 52, 52-57 (1984). 13. UPPAL A., RAY W. H. & POORE A. B., On the dynamical behaviour of continuous stirred tank reactors, Chem. Engng Sri. 29, 967-985 (1974). 14. KEENER J. P., Infinite period bifurcation in simple chemical reactors, in Mode//kg qf Chemical Rem/ion Sysrems (Edited by K. H. EBERT, P. DEUFLHARD and W. JACER), pp. 126-137. Springer, Berlin (1981). 15. REMBNYIK., KERT~SZ V., Vii~bs L. & HORV~~THF., Dynamical behavior of coal char particle combustion. To appear in Cornbust. Flame.
Melhods
& Applicmunr,
Vol.
17.
No
3, pp
267.283.
1991 r
DEGENERATE
Budapest
HOPF BIFURCATION
University
of Technology,
Miiegyetem
0362-546X.91 $3,001 .OO 1991 Pergamon Press plc
IN TWO DIMENSIONS
rkp.3,
H-1521
Hungary
and R. E. KOOIJ Delft University
of Technology,
PO Box 356, 2600 AJ Delft, The Netherlands
(Received 13 March 1989; received for publication Key words and phrases: oscillations.
Degenerate
bifurcation,
Hopf
1 February 1991)
bifurcation,
computer
algebra,
nonlinear
INTRODUCTION
CONSIDER the differential
equation
in IR2 with one real parameter i = “0,
x, Y) (1)
j, = g(p,x,A
1’
wheref(p, 0, 0) = g(p, 0, 0) = 0, ,D E IR. Suppose at a critical parameter value, p = p,,, a Hopf bifurcation occurs. For the sake of simplicity let p, = 0 and let the coefficient matrix of the linear part of (1) at ,D = 0 be
[
0w 0#0. --o 0I;
The type of this bifurcation (super- or subcritical) is determined by the stability ,D = 0. If the first Poincare-Lyapunov constant g, is negative, then the origin and if g, > 0, then it is a repellor. The formula for g, is well known:
of the origin at is an attractor,
where the partial derivatives are to be taken at (p, x, y) = (0,0,O). If g, = 0, then a degeneracy occurs and the sign of the first nonvanishing PoincareLyapunov constant g, determines the stability of the origin. The greater the number k, the more drastically would increase the complexity of the formula determining g, as a function of the appropriate partial derivatives at (0, 0,O). We introduce here an algorithm for the determination of g, for any k (Section I). We use this algorithm in an algebraic computer program written in Macsyma 420 language and accomplished on a Symbolics 3620 computer. The computation times achieved are surprisingly short. The algorithm is rather simple and uses only elementary vector and matrix operations. However, the exact proof (Section 3) needs quite new mathematical apparatus (Section 2). This apparatus works with multilinear maps. Here we introduce these new tools only for the 267
268
V. KERTI!SI
R. E. KOOIJ
and
two-dimensional case where multilinear maps are represented by coordinates. These coordinate representations will be called arrays. Mathematical results for the n-dimensional case can be found in [l]. It seems the apparatus works in the infinite dimensional case as well. In the case of equation (I), the degenerate Hopf bifurcation is not typical and may not be of much interest. However, if (1) consists of more than one, say k, parameters then in the kdimensional parameter space there exists typically a point where the first k-l PoincareLyapunov constants vanish. This point, as an organizing center, determines the local dynamical behavior of the differential equation in its neighborhood. This behavior is typical and indispensable to revealing the global dynamical behavior. This is the reason that great efforts have been made to obtain an effective way for determining the type of the Hopf bifurcation [2-121. To find the organizing center one need the parametric handle of the algorithm used. For this purpose our algorithm has been incorporated into the algebraic programming system mentioned above. In Section 4 we show some simple examples with and without parameters. In Section 5, as an example of practical importance, a mathematical model with three parameters will be studied. It has already been discussed in several papers [ 13- 151; still its dynamical behavior is so complex that it has not yet succeeded in revealing the entire phase portrait. In this section the nondegenerate and degenerate Hopf bifurcation of this model will be thoroughly investigated. I. ALGORITHhl
Consider
FOR
DETI~RR1IN,4TiON
the two-dimensional
OF
differential
THE
POINC‘ARE-I
Y.4PUNOV
CONSTANTS
equation
k = f‘(s, y) (2) _ii = g(.u, y) 1 ’ where
f and
g
are
smooth
functions
and
f(0,0) = g(0,0) = f,(O, 0) = gv(O, 0) = 0,
f,.(O,O) = -g,(O,O)= 0 # 0. Let A be an n x 111 matrix with elements u!,. We shall use the notation A(i,j) = CI,,. Similarly, if cl E R” and c7 = (~1,, ~1~). . , [I,,), then let v(j) = L;,. T denotes the transpose. k 0n
1s the binomial
coefficient
defined
as
k! (k -
0 THEOREM
1
bk+, 9 L’h+I
ifO
(algorithm). Define the vectors and matrices c’?, e,,,,, 12A+3, Ezxtr, and dh _, ; k = 2, 3, . .; I = 1, 2, . . . ; k + I 2 4 as follows:
U2’= (l,O,
1); 1
ifi=
1
Ck(i) =
e, E R!“;
i0 I, is the
n)‘n’ . .
k
x
k identity
otherwise; matrix;
CZht2, A,,,,
269
Degenerate Hopf bifurcation in two dimensions
E,, is a (2k + 1) x (2k) matrix,
&(i,j)
=
1
ifi=j
0
otherwise;
1 CZk is a (2k) x (2k + 1) matrix,
G(i,j)
ifi=j ifj=2k+
=
landi=2/+
1
otherwise;
A,,, is a (k + I) x (k + 1) matrix; k +
a/+j-i- I
x
a;+j-i-l
1
_
ai+l-j $+1-J
g
j
al+j-ia’-.if
(0, 0);
bk, vk E IRk+‘; k-l bk
=
c /=2
h,k+l-/U/i
if k is odd
--/+,bk uk = i
’
if k is even;
-~&(C,A,,,&-lC,bl
where the inverses
always exist, and
d, E R2k+3; dk’ = The Poincare-Lyapunov
elk+3[z2k+3
constants
Note. The vectors dk are constants
-
A2k+2.,E2k+2(C2k+2A2k+2,,E2k+2)-1C2k+21.
of equations
(1) are gk = dTb2k+2.
not depending
d,r = $(3, 0, 1, 0, 3), 2. INDICES
on f and g. For example:
d,T = &(5, 0, 1, 0, 1, 0, 5). AND ARRAYS
and IK'= (1,2). Let lKnC1= IK"x IK', n E N. If i E IK" then let Let IN = {1,2,3,...) i = (i,,i, , . . . , i,), where ik E IK',1 5 k 5 n. We call i an index from IK"and i, is called the kth element if i. If i E IK"and j E IK'" then let (i, j) = (i,,i,,...,i,,j,,j,, ...,j,,)E lK"+"'. For convenience we use IK"as the empty set, and for n and m above it is allowed to equal zero. If n = 0 then (i, j) = (j, i) = j. If n = m = 0 then (i, j) is an “element” of IK'.Let ,$J= UnEN IK". 3 will be called the index set. Now we define some functions.
Definition
1. p: 3 x N -+ NJ; p(i, j) = m if exactly m elements
Definirion 2. 0’: 3 + IN; a(i) = p(i, 1) + p(i, 2).
of i have the value j.
270
V. KERT~SZ and R. E. Koou
Definition 3. y: 3 --t 2”; y(i) = {j E IK”(‘)1v k E IK’ p(i, k) = p(j, k)]. Let H be an arbitrary finite set. Then let N(H) denote the number of the elements
of H.
Definition 4. cp: 3 + N; q(i) = N(y(i)). The following relations are easy to prove:
Co(i) (o(
_ P(i, t) a(i)
+ 1 + 1 ’
where t E IK’.
Let R be the following equivalence relation on 3: i, j E 3, i = Rj if and only if y(i) = y(j). Let 3 j R be the set of equivalence classes defined by R. (We can define IK” 1R in the same way.) Obviously, if i = Rj then p(i, k) = p(j, k), k E IK’; a(i) = a(j); p(i) = r&j). Consequently, 0, cp and p(. , k) can be considered as functions defined on 3 1R. Obviously, iV( IK” 1R) = rn + 1. Definition 5. v: $1 R + iN; v(i) = p(i, 2) + 1. For any fixed k E IN, the restriction of v to IK” 1R is a bijection, this set. Let 6(k, .) be the inverse of this bijection.
namely
an enumeration
of
Definition 6. 79: N x N + JJ I R; 29(k, i) is the unique element of IKk I R for which the relation v(6(k, i)) = i holds, i 5 k + 1. We can derive an explicit formula for 6: 6(k, m) = (i, j), where the indices i, j E 3 I R are uniquely determined by the relations o(i) = p(i, 1) = k + 1 ~ m and cr(j) = p(j, 2) = m - 1. Functions of type A : IK” + R will be called arrays. The integer n is the type of array. The set of arrays of type n will be denoted by a”. Obviously, Q2 is the set of 2 x 2 matrices, and 0.’ is R2. Let 6!” be R. If an element of i E IK” is fixed then A E 6!” as a function of the remaining elements is obviously an element of Q”-‘. If i E IK” and i = (i, , i,, , i,,), then instead of A(i) if i E IK”] and or A((;, , iz, . . ., i,!)) we shall write simply /I(;, , i,, . . . . i,,), as well. Similarly, j E IK”1, where n, + n2 = n, then instead of A((i,j)) we shall write A(i,j). Note that if n = 1 or n = 2, then these notations correspond to vectors or matrices, respectively. Q” as a set of functions with real values is a linear space using the necessary usual way: (A + B)(i) = A(i) + B(i), etc. Definition 7. The array A E (I?“, n 2 2 is called (i) symmetrical, if A(i) = A(j), wherej E v(i), i E IK”; (ii) quasisymmetrical, if either n = 2 or n > 2 and fixing i, , A as a function is a symmetrical element of a”-‘; (iii) antisymmetrical, if xJt7(,) A(j) = 0, i E K”.
definitions
in the
of i,, i,, . . . , i,,
The set of symmetrical, quasisymmetrical and antisymmetrical elements of Q” will be denoted by C”, Q” and Z”, respectively. These are obviously linear subspaces of Q”. If S E [” then it can be considered as a function of type IK” / R - I?, as well. THEOREM
2. Q” = i” @ 2”.
Degenerate
Proof.
Hopf
bifurcation
Let A E Q” and S and 2 be defined S(i) = $)
c
as follows: Z(i)
A(
271
in two dimensions
= A(i)
i E IK".
- S(i);
JET(i)
Then A = S + Z and S is symmetrical, Z is antisymmetrical. We now show that this decomposition is unique. Suppose A = S, + Z,, where S, E [” and Z, E 2”. Then S - S, = Z, - Z. Let Sz = S - S, . S2 is symmetrical and antisymmetrical at the same time. Consequently C S(j) = &i)&(i) = 0. j E y(i); S,(i) = S,(j); jEy(r) Considering
that p(i) # 0, it follows
Let the symmetrical
component
n
S, = 0.
of A E Q.’ be denoted
by $&4).
8. Let A E a” and B E a.“‘, nm L 1. Then AB E @.n+mm2and
Definition
AB(i,j)
=
c
i E IKn-',
A(i, kVW,_i),
j
E
IK"-'.
k=l
Distributive
and associative
A(B provided
laws can be easily verified: (B + C)A
+ C) = AB + AC,
= BA + CA,
(AB)C
= A(BC),
in the last case B is not a vector.
Consider the special case where A E a”, n > 2, uk E R2, 1 5 k 5 p 5 n. The multiple product (. . .((Au,)u,)u, . . .)u,, will be denoted simply by Au, u2, . . ., up. If U, = u2 = ... = up = U, then we write Au”. Au” is a homogeneous polynomial of the elements of U. The degree of this polynomial is n. More preciously: Au” =
c
A(i)u~“‘r)U~(‘,‘)
;clK"
=
c
rp(i)S(A)(i)u~“.“u,P”.”
= $JA)u”;
uT = (u, ,
;tlK"IR
Zu”=Oforeveryu~/R~ifandonlyifZ~Z”.LetS~~”andu,~IR~,k= Then Su,, uk,, . . . , uk, = A%,4, where k, , k,, . . . , kp is any permutation
1,2,...,p
of 1,2, . . . ,p.
Let S E [” and let u be a differentiable
function
$ Su”(t)
of type R -+ I?‘. Then
= nSu”-‘(t&(t).
PROPOSITION 3. Let S E c”, B E ak; n, k > 1 and u E R2. Then
(Sun-‘)(Buk-‘) Proof.
(SU”-‘)(Buk-‘)
= Suu, . . . , u(Buk-‘)
= (SB)U”+~-~. = S(Buk-‘)u,
Now we define an Q” 4 R-type linear operator.
. . _, u = (SB)U”+~-~.
u2).
V.
272
Definition
R. E. KOOIJ
KERT~SZ and
9. Let YE O”, n 2 2. Then
v E I?“+’ and
Vi) = d6(n, i)E( V(W, Some basic properties
of this operator
are: ,,+ I Vu” = C V(i)u, p(run.iJ.1)u2/Hij(fl,l),2).
v;
giq=
0).
i=,
the map s(V)
+ v is a bijection; Vu” is identically zero if and only if v = 0, where V E (3”. Let now V E azk and uT = (u,, u,). VuZk = V(1, 1, . . , l)(uf + u:)~ for every u E I?* if and only if C,,v = 0. (For C,, see Section 1.) Definition (k + n -
10. Let 2 5 k be a given integer. Then for any determined by the equation
A E Q”, n 1 2, X,(A)
1) x (k + 1) matrix
where Wi E i”;
W,(I) =
PROPOSITION 4. Let A E (IV, n 2 2. Then
Proof.
Using the definition
1/do
if v(/) = j
0
otherwise.
for every
V E ck, k 2 2, VA = 3Kk(A)I/.
of W,, we have k+l v
=
c ,=I
w,w,,
which has the consequence h+I
VA(i) The equality
ktl
= 1 &j)WjA(i) ,= I
of the first and third expression
PROPOSITION 5. Let A, E
= c mk(A)(i,j)@j). ,= I is equivalent
n
to the statement.
Q” ‘, then
where p(s, , 1) = I + j - i, p(s, , 2) = i ~ j, p(s,, I)=/+,;-i-
l,p(s,,2)=i+
1 -j.
Proof. ntk(A,)(i,j)
= (W,A,)(i) _
= cp(ti(k + I ~ 1, i)S(W,A,)(ti(k (W/$)(p)
c
pt-p(ii(h+/-1.1)
_
cp(r,Ms,) v((r,,
1))
A,(l,s,)
=
c
'I,, ty(rY(k+l-1.1)
cP(r2)ds2) + (D((rL) 2)) A,@, sz),
+ I -
1, i))
[W,(r, l)A,(l,s)
+ y(r,
2M,(2,s)l
is a
Degenerate
Hopf
bifurcation
in two dimensions
273
where
(r,, 1) E ~(W,j)h (r2, 2)
dW,j)).
E
Note that k+l-j Ar,, 1) + 1 &r,) cp((r, , 1)) = a@,) + 1 = k Ar2, 2) + 1 _jdrd v)((r2, 2)) = a(r2) + 1 As,, 1) = &VI,
1. ’
k
s,>, 1) - Ar,,
1) = j + I - i;
As,, 1) = A(r2, s,>, 1) - Ar,, O(S,) = c&)
P(Sl)=
(p;;fl)) =(if j);
ds2)
(,~yJ
=
(2) can be rewritten
=
OF THE
+
Ak(l, i) =
n
J.
ALGORITHM
off
+ Al,22 + ... + A,d
and g with an appropriate
r E R\J:
+ o(lulr+‘),
pci, l)+P(i,2)
f a
k’ap(‘.1)xap”,2’y
a positive
definite u(u) = 1,and
whereVkEck;2
0);
polynomial
Ak
i)
=
F
1 (pci, l)+P(1,2) (0, 0). pci, ljx ap(i,2)y g
function
u: R2 + R of the form
v2u2+ v,u3 + ... + v,+,Ur+l, 0
[ 0
v,= the derivatives
t2,
Lyapunov
1
Compute
:
A, E Qk+‘, 1 5 k 5 r, 1
Now we define
(i
using the expansion
ic = A,u where 1.4’ = (x,y),
1) = j + I - i - 1;
= I;
3. PROOF
Equation
’
of u along the solutions
of (2’):
r+l + c (kVkA,
ti(U) = (2&4,)U2
k=3
where k-l Bk
=
c /=2
1
1’
II/IAk+l-,.
+ &)Uk
+ o(IU1r+2),
(2’)
274
V. KERT~SZand R. E. Koorr
Compute
312, (A ,). Obviously
-~ i k f t%,(A,)(i,j)
=
k +2 - i k \O
ifj=
i-
1
otherwise.
If k is odd then let (kV,A, + B,)uk = 0. Or, equivalently: knt,(A,)V’ + Bk = 0. Using the Gaussian elimination, we can easily see that Fm,(A,) is regular. Consequently: i;;, = -(l/k)X,(A,)-‘B,. If k is even then let (k&A, + B,)uk = c(x2 + y2)k’2, c E IF?.Or, equivalently: C,(kZVZ, (A ,)V, + Bk) = 0. Let the (k + 1)th element of Vk be zero. Then Ek ErV;. . Note that C, 3nk (A ,)E, is regular (use the Gaussian elimination). Now we can express _ v, : vk = -~El(Cx~,(A,)~~)-‘C,17,. Lastly,
note that k-l Bk
consider theorem
=
1 I=2
proposition 5 and use the notations 1 is completed. n
3. COMPUTATION
WVk+,-,I&;
u& = vk, bk = Bk, A, , = %,(A,).
BY
COMPUTER
The proof
of
ALGEBRA
We developed a computer program (on a Symbolics 3620 computer written in Macsyma 420 language) which carries out the algorithm explained in Section 1. The program computes the Lyapunov function and the Poincare-Lyapunov constants. We introduce here some computational examples. The simple examples 1 and 2 were solved by hand in order to show how the algorithm works. Examples 3 and 4 are more complex and were carried out with our computer program. A rather sophisticated example is shown in the next section.
Example 1. .t= -y+x2 j, The partial
derivatives
=x+x2.
in the origin are: f,,(O, 0) = g,(O,
0) = 2. Any other derivative
is zero.
Degenerate
A 3.1 =
s-i;
Hopf
bifurcation
_;;];
in two dimensions
-%t=i[i;
275
;;I;
f
I
uj = -fA,‘,
01
A 3,2
c
1
-2
6, =
I1 -41I 0
b4 = 20, + 3A,,,q
;
03 = 0;
I-
+ 50
$00
0
00’
-4
oot1
=
’
b, =
-4
;
g, = +(3
0
1
0
3)b, = -+.
0 o_Ii Check the result by computing
u4
10
-4 A 4.1
=
0
and derivating
0I the Lyapunov
2
00 00 -+E4(C4A4,,E4)-‘Cd
$
o-2
function.
-+E,(C,A,,,E,)-‘Cab,
=
0
E4 =
[
04
= 4
! 1 -3 -50
02
-101
042
053
;
216
+ I1
V. KEKTESZ and R. E. KOOIJ
X3
x3
-4)
VJ =
X2Y
(-1
-*
x’y
0 + 0)
XY2
x2y2
X_$
YS
Example
2,
0-
-1 6, = 2a,
v3 =
LIZ
1 -I_
6, = 3A,,2v3 :=
-2-
00
3 a, =
0
0
1
0
0
6
0
2
0
I 3 0
3
3
0
3 -2
0
2
0
0
0
0
1
-3 li
2
-?
= 1
0
0
0
0
1
YJ
= I!‘I 0
2 3
6
0
2
-3
o_ 0
-2
-2
_z 3
0 vq=
0
.
0 !I 0
Note that if bk = 0 then L’~= 0 and b,,, = kA,,, vk, k = 4,5, . . . . We can conclude vk = 0 if k 2 4 and g, = 0 for every k. Consequently, v = x2 i y2 + fx’ - +yv’ - 2xy2 is the first integral
of our equation.
Actually:
ti = 0.
Example 3. i=
-_v+x”
j = .Y g, = g, = g, = g, = 0;
231 g5 = m
7337 793 U = x2 + y2 - 512x1’_v - 1536xq’3
9273 - mx’y5
231 7623 1309 - ~ x5y7 ___ x3y9 - 512xy1’. 1280 512
Degenerate
Hopf
bifurcation
271
in two dimensions
Example 4 (see [2]). i = -y - ax2 + 2xy + by2 j = x + x2 + (5b - 3a)xy - y2 g, = +$(b - a)3(2b2 - ab + 1)
g, = g, = 0; 5. INVESTIGATION
Consider
the differential
OF A REACTION-KINETIC
MODEL
equation T = -/3T + BDaCe’, C,T>O
c = 1 - C - DaCer;
1 ’
(3)
where p, B and Da are positive constants (parameters). This equation represents the first-order exothermic chemical reaction A + B in a continuously stirred tank reactor and was investigated in great detail in [ 131 and [ 141. The same equation is used in [ 151 as a simplified model for a coal char particle combustion. Equation (3) possesses an extremely complex dynamical behavior which has not yet been entirely revealed. In order to be able to reveal the entire global dynamic behavior one needs the complete local investigation. Amongst others, local investigation contains the Hopf bifurcation study. Degenerate cases are not typical, still, studying the degenerate Hopf bifurcation, we can reveal typical bifurcation phenomena around the degenerate point as the organizing center. Let us consider a system with two real parameters:
k = f(Q, P, x, Y) i
=
1
(4)
g(& P, x, Y).
Suppose there is a function @: R + R such that if O( = @(/I) then a Hopf bifurcation occurs in (4). Assume moreover there exists a & such that if /I < /IO then the bifurcation is supercritical (g, < 0) and if /3 > /3(, then subcritical (g, > 0). Assume the periodic orbits appear in the supercritical case if 01 > Q(p) and in the subcritical case if (Y <
+ f ye” -
WDae”.
depends
only on two parameters,
V. KERT~Z and R. E. Koorr
278
At1
A/i = const.
t
UE: unstable ULC: unstable
SE: stable
equilibrium limit cycle
SLC: stable
Fig.
equilibrium limit
cycle
I.
Note that 0 < w < 1. We cannot express 0 and I,Uexplicitly as functions of ,Qand Da. Therefore we use ,8, 0 and IJ/ as new parameters, and we introduce x and y as new variables: x=T-H y=c-rp
By these, equation
(3) yields:
.k = P(H - 1)X +
‘fy +p&(e’
-
1 - x) + iy(e’
-
1))
(5) j, = (cc/ - 1)x - ty
+ (w - I)((e’ - I - x) + iy(e’
-
1)).
1
Degenerate
The coefficient
matrix
Hopf
of the linearized
bifurcation
279
in two dimensions
system is
A=
The necessary
conditions
of the Hopf bifurcation
1<8<1 w Substituting
are: det A > 0 and tr A = 0, that is p=
and
l (0 - L)cy’
the critical /I value in (5) we obtain 1 8 i = w” + (0 _ 1)w2Y
j = (tf - 1)x - :Y In order to obtain
the canonical
I9
1 - X) + +_v(e” - 1 - x)
(e" -
+ (0 -
lb
1
I
+ (w - 1) (e” - 1 - x) + t (e’ - 1) . I 1 form,
we transform
the variables:
As a result: if=ov+g (6)
)
i, = -_ou ii cog 1 where 1 + wQ2 g = w exp (1 _ vI)v2 (
’ +~2v~u) (1 - w)w
u-l-
+ (cou -
iU)(exp(:+_$G2z4
-
O
The next step is to determine w and w: g, = (w
2c02 + 1)2(l/o4
the Poincare-Lyapunov
+ 24U3w2 + l&02
- I+2
constants
1)
0 < 0.
g,,g, and g, as functions
of
- 2l/Y3 + 74!/2 - 5l+U+ 1)
S(I,V - 1)4y/6w2 g, = -(v/‘w2
+ 1)(186~8~‘o
1967y/‘o”
- 3051,~~o~ + 13951,~‘o~ - 145w6w8 - 1861,~~~~ + 1621,~~~~
+ 49561~1’0~ - 1244~‘~~
- 1098~~0~
- 6289~~~11~ + 12843~~~~
- 7768~~0~
+ 2808~~~~
+ 201661~1~0~ - 185041//‘w’
-
10819~5w2
+ 1860~ - 1201,~’ + 540~~ - 750~’
•t 345~‘~~
+ 3671y3w6 + 990~‘~~
+ 465~’
-
f 899w2w4 - 186tj~0” - 276~‘~~ + 8718~~~~
135~ + 15)/(1152(1~
- 2033~0’ -
1)8~‘206)
280
V. KERT~
g, = (‘//‘w’ + 1)“(1341751//‘2w’6 - 134175~9w’4
-
- 27960831t~/‘~w’~
- 89543353~“~’ + 8876630~swx + 66376480~‘~~ + 106174634~sw6 134175ww6 -
153360~“~’
- 474652~‘~~‘~
+ 6389651j/‘~w’~
- 3855099~‘~‘~ 6082518’~/“w’~
+ 429615O@w” 264523371y”w’
and R. E. KOOII
- 5378394~“~‘~
1219977~~~‘~
+ 46788711’/1’w”’ - 753500~~~‘~
-
+ 522950~~0~
185245120~8w”
-
-
+ 1032312’~1’w~ - 4620172~‘~~ + 7022916’~/~w’
+ 88800~’ 36075~
- 488400@
2775)/(442368(~
+ 932784~“~~
1375224~‘~~
18746700~‘~~
+ 9356650~~~~
- 2949788~‘~’
+ 976800~’
-
12385480’/“w6
- 237643272~~~~ + 1173485~~~~
-
+ 74127618’j/‘w4
- 268137446~~~~
+ 162416891y2w4 - 2246523~~~
+
+ 15359604~6ws
+ 283247706~‘~~
- 20262648~~~~
+
- 2963550~‘~‘~
- 3266418’//“w’
+ 13828311 l’,//‘w’ - 90400342~‘~~
- 180498774’+~‘w” + 276233900’/‘0
+ 36075~ 7308~’
- 3948321//“w’~
- 22318196~~~‘~
+ 88296~‘~~~
- 4255655~“~’
- 341149’$“w’4
+ 11365519~‘“w’2
+ 5026841~‘~‘~
+ 26360401y’“w4
- 66201688~‘~”
+ 1762334~“~‘~
+ 166747467~~~~
+ 1341750”
- 74544~~~~
- 10434862~~0~
+ 763343~~~~
-
197025~~
- 112505’//w2 +
- 999000’/1~ + 582750~”
-
197025ry’ +
- l)‘21y’xw”‘).
Note. Equation (6) contains two parameters, namely I,Yand w. In the general case of two parameters the typical situation is the following. On the parameter plane there is a curve which separates the domains where g, > 0 and g, < 0, respectively. Along this curve g, = 0. On the curve there is a point where g, = 0 and which separates the points where g, > 0 and g, < 0, respectively. We need g, only in this point. In case of equation (6) this typical situation does not occur as the point g, = 0 is outside of the possible domain of the parameters. We have shown the expression of g, above only in order to demonstrate the possibilities of our algebraic program. Continuing the computation we solve the equation g, = 0. As a result we obtain on the parameter plane (see Fig. 2) which have the equations y, : w = [l ~ ‘// ~ 2~’
+ (1 ~ ~)(81y
- 8~’ + l)“2/(2y/3)]“2
y2: w = [l - y/ - 2’~/~ - (1 - w)(8w3 - 8~’ + 1)“2/(2w3)]“2. The sign of g, outside
yl and y2 can be easily determined
(see the figure).
two curves
(7) (8)
Degenerate
Hopf
bifurcation
in two dimensions
Fig. 2.
b)
Da p1 = e”‘/H,* Da,=e”‘!(H,H, =3+
L-j 2
Fig. 3.
1)
281
282
V. KERT~SZ and R. E. KOOIJ
Now we determine the sign of g, along y, and yz . First we leave out the factors with positive sign from g, and then we substitute the expressions (7) and (8) for w. So we obtain g,, and gzz, respectively. g,, = (-64~~’ + 1282~”
+(8w3
- 81y2 + 1)“2(24v”
- 2496~~
+ 1736~~ - 3801,~~ -
- 35241,~” + 9240~” - 9001~1~- 98~’
- 228~”
- 14002~8
+ 111~ -
g,, = -(641c/12 + (8~~ - Q2
+ 8321,~’ - 1402’,~’ + 7081,~’
1704~~ + 112’~ - 18) + 72&//l’
+ 116601y’ - 3302~~
- 2840~’
+ 3008~~
18)/@
+ 1)“2(24w”
- 2281,~“’ + 832~~ -
1402~~ + 708~’
+ 1282ry’ - 24961,~~ + 17361~’ - 380~~ ~ 17Oy/’ + 112~ - 18) - 728~” + 3524~“’
- 9240~’
+ 900’//3 + 98I$ A further
computation
-
+ 14002’//’ - 1166Oy’
+ 3302~~
+ 2840~~
- 3008~”
112l/Y + 18)/‘J/h.
yields the result that g, < 0 on yl and g, > 0 on y2.
As a summary, consider Fig. 3(a) where the coordinate B of the equilibrium is plotted vs the parameters Da and ,u. The Hopf bifurcation diagram is shown on the @Da, ,u) surface. The projection of this diagram on to the Da, p plane can be seen in Fig. 3(b). At the critical value PC,,, = ((0 - l)cy)-’ of the third parameter p, a super- or subcritical Hopf bifurcation occurs according to the figure. On the curve g, = 0 the Hopf bifurcation is degenerate. The sign of g, determines the type of this degenerate Hopf bifurcation (see Fig. 2 as an illustration).
REFERENCES 1. KERT~SZ V., hlultilinear maps and bifurcation theory. To appear in Co//. math. Sot. Jdnos Solyur, Szeged (Hungary). 2. BAUTIN N. N., On the number of limit cycles which appear with the variation of coefficients from an equilibrium position of focus or center type, Mal. Sh. 30, (72) 181-196 (1952); Am. marh. Sot. Trunsl. No. 100 (1954). (In Russian.) 3. ANDRO~.OW A. A., LEONTOVI~~IE. A., GVRDON I. I. & MAIEK A. G., Qua/r/u/we Theory of Second Order D_vnarnicul S_vslems. Halsted Press, New York (1973). 4. MARSDEN J. E. & RICCRAC‘KEN bl., T/?eHopj’Bi~urcution und irs App/icotiom. Springer, New York (1976). 5. HASSARD B. & WAN Y. H.. Bifurcation formulae derived from center manifold theory, J. math. Ana/.vsis Applic. 63, 297-312 (1978). 6. SCHMIDT D. S., Hopf’s bifurcation theorem and the center theorem of Lyapunov with resonance cases, J. tnafh. .4na/.vsis Apphc. 63, 354-370 (I 978). 7. GIBBER F. 8; WILLAMOWSKI K. D., Lyapunov approach to multiple Hopf bifurcation, J. mar/z. Analysis Applrc. 71, 333-350 (1979). 8. NEGRINI P. & SALVADORI L.. Attractivity and Hopf bifurcation, Nonlinear Anu/_ysi;, 3, 87-99 (1979). 9. GOLLBITSKY hl. & LANWORD W. F., Clarsification and unfolding of degenerate Hopf bifurcation, J. d/yf. Eqns41, 357-415 (1981). 10. SHI SONGLIXG, A method of constructing cycles without contact around a weak focus, J. d/j__ Eyns 41, 301-312 (1981). 11. HASSERD B. D., KAZARINO~~ N. D. & \VAN Y-H., Theorj~andApplicaiion of HogfBifirrcarion. Cambridge Univertity Pres\, Cambridge (1981).
Degenerate
Hopf
bifurcation
in two dimensions
283
12. SHI SONGLING, On the structure of Poincav-Lyapunov constants for weak focus of Polynomial vector fields, J. diff. Eqns 52, 52-57 (1984). 13. UPPAL A., RAY W. H. & POORE A. B., On the dynamical behaviour of continuous stirred tank reactors, Chem. Engng Sri. 29, 967-985 (1974). 14. KEENER J. P., Infinite period bifurcation in simple chemical reactors, in Mode//kg qf Chemical Rem/ion Sysrems (Edited by K. H. EBERT, P. DEUFLHARD and W. JACER), pp. 126-137. Springer, Berlin (1981). 15. REMBNYIK., KERT~SZ V., Vii~bs L. & HORV~~THF., Dynamical behavior of coal char particle combustion. To appear in Cornbust. Flame.