Computing Hopf Bifurcations I.
This paper addresses the problems of detecting Hopf bifurcations in systems of ordinary differential equations and following curves of Hopf points in two-parameter families of vector fields.  The established approach to this problem relies upon augmenting the equilibrium condition so that a Hopf bifurcation occurs at an isolated, regular point of the extended system. We propose two new methods of this type based on classical algebraic results regarding the roots of polynomial equations and properties of Kronecker products for matrices.  In addition to their utility as augmented systems for use with standard Newton-type continuation methods, they are also particularly well adapted for solution by computer algebra techniques for vector fields of small or moderate dimension.
Introduction
. Consider the n-dimensional system of ordinary differential
equations defined by
(1)
is a C 2 \Gammasmooth function, defined on a subset
U ae IR n , which depends upon a vector of parameters, ff 2 IR k . An equilibrium of
this system is a point x 2 IR n with the property that f(x; As the parameters
ff are varied, equilibrium points can undergo bifurcation. There are two types of elementary
bifurcations that occur in generic one-parameter families of systems: saddle
nodes and Hopf bifurcations. At a saddle-node bifurcation, the Jacobian derivative
evaluated at an equilibrium point possesses a simple zero eigenvalue. Similarly, a necessary
condition for Hopf bifurcation is the presence of a pure imaginary eigenvalue
pair in the spectrum of D x f . One would like robust algorithms for the calculation
of parameter values where one of these bifurcations occurs in parametrized families
of vector fields. In the case of saddle node bifurcations, one may obtain the bifurcation
locus by augmenting the equation f(x; with a procedure to calculate
equilibrium points where D x f has numerical rank (n \Gamma 1). An inflated system for the
detection of saddle nodes has been previously proposed using det(D x f) as the augmenting
equation [1, 28]. To treat Hopf bifurcations in a similar manner, one requires
explicit equations that determine whether the n-square matrix D x f has a pair of pure
imaginary eigenvalues. This paper examines procedures for locating Hopf bifurcations
based on the singularity of matrices obtained from algebraic transformations of the
Jacobian matrix at an equilibrium.
Previous investigators have proposed a variety of approaches to the determination
of Hopf bifurcation points, which may be divided into two broad classes. A typical
indirect method employs a numerical algorithm for computing the spectrum of the
Jacobian at each point along a path of equilibria and then interpolating to locate
parameter values at which a pair of eigenvalues cross the imaginary axis. For vector
fields of small dimension, the unsymmetric QR factorization algorithm is widely ad-
vocated, while for larger problems Krylov subspace techniques are fast and provide
Mathematics Department and Center for Applied Mathematics, Cornell University
y Center for Applied Mathematics, Cornell University
z Mathematics Department, Cornell University
adequate accuracy in few iterations [29, 13]. In a direct method, equations which
characterize the Hopf point augment the usual equilibrium condition and standard
techniques are used to solve for roots of the resulting system. Algorithms frequently
used in practice solve explicitly for the purely imaginary eigenvalues together with a
normalized basis of the corresponding eigenspace [17, 31, 32, 22]. Such an approach
may be viewed as simultaneously solving for an equilibrium point and an eigensub-
space by an iterative method. In this sense, transition from indirect to direct methods
has replaced the use of highly stable and robust linear algebra techniques for spectral
computations with root-finding algorithms which are only locally convergent. Such
a strategy has a drawback. The minimal dimension of the inflated system increases
to at least 2n + 1, which may decrease the performance expected from the algorithm
used to solve the resulting nonlinear system of equations and increase the demand for
good initial values and small stepsizes. A side effect of these methods is that they
produce more information than required - the detection of a pair of purely imaginary
eigenvalues - at a significant computational cost.
The purpose of this paper is to examine the algebraic structure of Hopf bifurcation
computations and use the insight gained to construct direct methods that rely
upon minimally augmented systems for the computation of Hopf bifurcations. Our
methods are based on the algebraic properties of polynomial resultants and matrix
symmetric products. This algebraic theory has been applied to many areas of applied
mathematics, in particular to develop criticality conditions in Lyapunov stability analysis
and root location [4]. The primary focus here is to develop a theoretical basis
for algorithms to study Hopf bifurcation in multiparameter vector fields of modest
dimensions where a complete and detailed description of the dynamics is a reasonable
objective. Our experience in this context suggests the approach may prove fruitful
for three reasons: First, the bifurcation analysis for a multiparameter dynamical system
typically proceeds through a number of phases, the earliest aimed at providing
a gross decomposition of the product space of phase and parameter variables into
regions of greater or lesser interest. Criteria for Hopf bifurcation that can be easily
evaluated based on the fixed point Jacobian alone, without requiring the details of its
spectrum, lead to fast algorithms for "scanning" the parameter space for candidate
regions of interest. Second, recent developments in computer algebra techniques and
the widespread availability of powerful computing machines require reassessment of
mathematical algorithms normally reserved exclusively for problems of very small di-
mension. Application of computer algebra methods to problems in dynamical systems
has already proved effective in normal form computation, center-manifold reduction
and perturbation analysis [30]. It is clear that approaches to computing Hopf bifurcations
which require the use of iterative algorithms or direct evaluation of the
eigenvalues are unsuitable for use with symbolic methods. If, instead, the Hopf condition
is expressed in terms of functions of the Jacobian then, for some problems,
application of symbolic methods become feasible. Third, steady progress in computer
hardware and software evolution promises to continue the trend towards very fast,
very robust methods for the solution of linear algebra problems. It seems desirable to
exploit reliable methods for matrix manipulation in exchange for reducing the size or
difficulty inherent in the solution of nonlinear equations, even if the required penalty
is an increase in the dimension of the matrices involved. For some problems, especially
those with special structure, this approach may yield algorithms that are more
effective than those in widespread use. A companion paper to this one [20] studies
the implementation of these methods on three examples. Detailed comparisons with
other algorithms are discussed there, though we make a few comments about these
comparisons in the final section of this paper.
Methods for detecting Hopf bifurcation based on algebraic techniques have appeared
in the literature, and the present paper complements and extends those results
in a number of ways. In particular, the program LINLBF by Khibnik and his coworkers
[26] incorporates a method for following curves of Hopf bifurcations based on the
resultant of two polynomials which is similar to that described here. Our subresul-
tant condition described in Section 2 provides a way of discriminating between Hopf
points and other resonant saddle points, and can be used by either method. We
describe three formally different forms of the resultant criteria that lead to matrix
formulations with distinctly different properties and prove that they are equivalent.
The development of augmented systems based on tensor products is new, although
the important classical role occupied by tensors in the study of stability criteria for
ordinary differential equations is very suggestive. Our proof that Hopf bifurcations
are isolated solutions to the augmented equations in a neighborhood that does not
contain a degenerate point applies to augmented systems based on either resultants
or tensor products.
In the next section we recall the algebraic preliminaries required to express the
necessary Hopf condition in terms of the singular set for certain algebraic functions of
either the Jacobian at an equilibrium point or its associated characteristic polynomial.
We present two new methods for following curves of Hopf bifurcation points in 2-
parameter families, and discuss the variations which arise due to the various ways
singularity of a matrix may be measured. Concluding remarks, comments about the
performance of our algorithms and a discussion of methods for extending this approach
to bifurcation points of greater degeneracy are collected together in Section 4.
2. Hopf Algorithms Using Polynomial Resultants.
2.1. Algebra of Matrix Resultants. We seek explicit criteria that specify
whether an n \Theta n matrix, with coefficients that may depend upon parameters, has a
pair of pure imaginary eigenvalues. In this section we develop the desired necessary
conditions expressed in terms of the corresponding characteristic polynomial for the
matrix.
Let J be the Jacobian matrix for f , and denote its characteristic polynomial by
Clearly, p has the non-zero root pair f-; \Gamma-g if and only if - is a common root of
the two equations p(-) + p(\Gamma-) and p(-) \Gamma p(\Gamma-). Making the substitution
and rearranging, we construct two new polynomials. If n is even, let
(2a)
n\Gamma2while if n is odd, set
Then p has a non-zero root pair f-; \Gamma-g if there exists a z that satisfies:
r
Two polynomials have a common root if and only if they share a common factor.
Since r e and r are univariate polynomials with the degree of r than or equal
to the degree of r e , the Euclidean algorithm may be used to answer this question.
We compute a sequence of polynomials with the property
that P i+1 is the remainder of dividing P i\Gamma1 by P i . Thus, one requires that the degrees
of the P i are strictly decreasing and that there are polynomials Q
that . The sequence terminates when P k is found that divides
. It is easy to verify by induction that the P i are all in the ideal generated by r e
and r o in the polynomial ring R[z] which implies that there are polynomials A i and
have no common factor, then the repeated
divisions will produce a P k that is a non-zero constant. On the other hand, if r e and
r do have a common factor, then their greatest common denominator appears as P k
which divides P k\Gamma1 . In this way the Euclidean algorithm may be used to determine
whether or not two polynomials share a common root.
The Euclidean algorithm is closely related to a collection of determinants constructed
from the coefficients of two polynomials. For n even, let the Sylvester matrix
of the pair of equations (2) be the (n \Gamma 1) \Theta (n \Gamma 1)-matrix given by
rows
nrows
while if n is odd, this matrix is defined to be
rows
rows
The determinant R is a polynomial in the coefficients of p referred to
as the Sylvester resultant of r e and r . The Sylvester matrix is singular (equivalently,
only if the polynomials r e and r share a common root. The
computation of the Sylvester resultant is usually performed through the construction
of a pseudo-remainder sequence for the pair of polynomials. Modulo units in the ring
of coefficients of the polynomials, the terms of the pseudo-remainder sequence are
the intermediate terms produced by the Euclidean algorithm. Explicit formulas for
the coefficients of these intermediate terms can be defined in terms of determinants of
appropriate submatrices of S, called subresultants. We shall use the subresultants that
correspond to a linear remainder term from the Euclidean algorithm. The coefficients
of the linear remainder may be obtained by considering determinants of two
3)-submatrices of S. For denote the matrix obtained from S
by deleting the rows 1 and nand the columns 1 and 2. Thus, for n even, S 0 is
obtained from the Sylvester matrix by removing the boxed entries shown below:
The relationship between the characteristic polynomial and its corresponding matrices
leads to the result:
Theorem 2.1. Let S be the Sylvester matrix for the polynomials r e and r o in
Equation 2. Then J has precisely one pair of pure imaginary eigenvalues if
then p(-) has no purely imaginary roots.
Proof: From the elementary properties of resultants [8, Chapter 3 - Prop. 8 and
share a common factor in IR[z]. Suppose that
. Then this common factor is linear and, by Loos [27], it has the
explicit form
for two polynomials IR[z]. The solution of det(S results
in a common complex root, -
z, of r e and r o which is unique, real and negative if
may have more than one pair of pure
imaginary roots, a case which may be resolved by computing higher subresultants.
When the conditions of the theorem are satisfied, the magnitude of the shared root -
z
is given by
which, in the case of Hopf bifurcation, can be related directly to the period of the
limit cycle created at the bifurcation point. An efficient method of evaluating all sub-
resultants (including the determinants of S, S 0 and S 1 ) is the generalized polynomial
remainder sequences algorithm due to Habicht, Brown and Collins. For a detailed
discussion see Loos [27].
In the case of Equation (2), the determinant of the Sylvester matrix can be
expressed explicitly in terms of the eigenvalues of J, according to the following result:
Theorem 2.2. Let - -n be the roots of the polynomial p(-) and let
be the (n \Gamma 1) \Theta (n \Gamma 1) Sylvester matrix of the associated pair (r e ; r
Y
1-i!j-n
Proof: We note that
Y
1-i!j-n
and the Sylvester resultant, R S , are symmetric polynomials in - . They can
be rewritten uniquely as polynomials in the coefficients c of the polynomial
Y
We denote these two polynomials by prod(c
spectively. Both polynomials have integer coefficients, and res is irreducible [34] in
The polynomials prod and res define the same hypersurface in the affine space C n .
Namely, a point lies in this hypersurface if and only if the polynomial
has a quadratic factor of the form - -. Hence there exists an integer, s, such that
res in Z[c
In order to show apply the classical algorithm for rewriting symmetric
functions in terms of elementary symmetric functions. We begin by noting
that the lexicographic leading term of
Y
1-i!j-n
. There is a unique monomial in the elementary symmetric
functions which has the same leading term, namely,
Therefore the monomial c 1 c 2 c 3 occurs with coefficient (\Gamma1) ( n
prod. How-
ever, the same monomial occurs with coefficient (\Gamma1) 1+2+:::+b(n\Gamma1)=2c in the expansion
of res = det(S), as can be verified from the definition of the matrix S. The desired
factor s is the quotient of these two coefficients. This quotient equals (\Gamma1) m and the
proof of Theorem 2.2 is complete.We note that Theorem 2.1 imposes two requirements on the characteristic polynomial
of J: a singularity condition for the resultant matrix and an inequality involving
subresultants. There are many classical results on root location in systems of polynomial
equations. For example, the well-known Routh-Hurwitz criterion for stability
matrices determines whether all the eigenvalues of a matrix have non-positive real
parts. We describe two other criteria for a polynomial to have a factor of the form
are appealing alternatives to the calculation of the Sylvester resultant.
For ease of exposition, we consider the case where n - 2 is even in the discussion that
similar results for the case where n is odd follow analogously.
Let M denote the space of all real n-square matrices and charpoly : M! IR n[x]
the map that associates an element of M to its characteristic polynomial. If N ae M
is the dense subset of matrices with simple eigenvalues, then the equivalence relation
defined by the condition that
induces a partition of N into similarity classes. For each equivalence class in this
partition, a useful representative element is a matrix that displays the coefficients
of the associated characteristic polynomial along a row or down a column; such an
element is called a companion matrix. Notice that for n even, r e is a monic polynomial.
Thus, we may take
C re =6 6 6 6 6 4
as the representative element of [C re ], since charpoly(C re . Moreover, for each
re there exists a (non-singular) matrix W such that
C re
Recall that a matrix polynomial of degree k with elements in M is an expression
of the form,
the coefficients a i are scalars and I is the unit matrix.
An alternative expression for R S is given by the following theorem [4, Section 1.2]:
Theorem 2.3. For r defined re is even and
r e (C ro is odd. If R C is used to denote det(C), then
Thus, the vanishing of R C provides an alternative to the condition R
in Theorem 2.1, but involves a matrix of smaller dimension. To develop a third
equivalent criterion, we consider the Bezout resultant. Again, for n even, consider the
two polynomials specified by Equation (2a). We define the brackets
nand we take c The Bezout matrix, B, corresponding
to the polynomial pair (r e ; r is an n-dimensional square, symmetric matrix with
entries constructed as sums of bracket products in the coefficients c i as follows: For
Then,
k=kmin
The only modification required in this definition for the case of n odd is that c n has
the value prescribed by the characteristic polynomial, p, and c n+1 is taken to be unity.
The Sylvester and Bezout matrices are intimately related (see Theorem 1.13 of
Barnett [4] and note that, in his notation, det(T 1 which leads to the following
theorem:
Theorem 2.4. For the Bezout matrix corresponding to the polynomial pair
det(B). Then we have
Moreover, the coefficients of the linear remainder term from the Euclidean algorithm
may be expressed as determinants of certain Bezout submatrices, in analogy to the
subresultants defined for the Sylvester form. Consider two submatrices, B i for
obtained by deleting the first column and the i th row of B. Then the Bezout
subresultants formed from the determinants of B 0 and B 1 , together with Equation 6,
may be used to specify the hypersurface of n-square matrices with a purely imaginary
eigenvalue pair. Since we have not found the particular relationship given in the
following theorem described in the literature, we provide a sketch of the proof:
Theorem 2.5. Let B be the Bezout matrix for the polynomials r e and r o in
Equation 2. Then J has precisely one pair of pure imaginary eigenvalues if
If RB 6= 0 or det(B 0 ) then p(-) has no purely imaginary roots.
Proof: The Bezout matrix B is a symmetric m \Theta m-matrix where
(Barnett [4], Section 1.5) it satisfies the following identity in the indeterminates w
and z:
Suppose that have a common
complex root, say w From Equation (7) we conclude
for all z 2 C. If we substitute m distinct values for z and invert the corresponding Vandermonde
matrix, then we see that Equation (8) implies that (1; w
is the zero vector. Now, by our assumptions, the Bezout matrix has rank m \Gamma 1. The
unique (up to scaling) vector in its kernel is computed by Cramer's rule. For some
non-zero constant c we find
This implies w negative if and only if
positive. Our proof is complete because w 0 is the only common
root of r e and r 1).Thus, we have three distinct - but equivalent - methods of computing the singularity
condition necessary to establish that the Jacobian has a pair of purely imaginary
eigenvalues: one (Sylvester) related directly to the sufficiency condition, one (Com-
panion) small in dimension compared to n and, finally, one (Bezout) which is not only
small but has symmetric structure. We illustrate the ideas presented above concerning
polynomial resultants with the following example.
Example: Consider the case Equation (2) becomes
r
The Sylvester matrix and the two relevant submatrices, S 0 and S 1 , are given by:
By Theorem 2.1, p has a pair of purely imaginary roots if and only if
1vanishes and the product
The singularity condition may be expressed in terms of the companion matrix for r e ,
through the matrix polynomial:
r re
Finally, the Bezout matrix associated with
and it is easy to verify . The Bezout subresultant condition
provided by the expression,
det
is necessarily positive when p has a pair of imaginary conjugate roots.

Table

1 provides a list of the resultant equality and subresultant inequality conditions,
as functions of the polynomial coefficients, for vector fields of dimension two to six.
3. Hopf Algorithms Using Symmetric Products.
3.1. Algebraic Properties of Symmetric Products. In practical applica-
tions, the computation of the characteristic polynomial coefficients is problematic
since known algorithms are numerically unstable. For each method (Danilewski, La
Verrier, etc.) otherwise unremarkable example matrices can be constructed for which
the technique yields arbitrarily inaccurate polynomial coefficients; for a detailed dis-
cussion, see Wilkinson [36]. Despite this fact, computational methods which employ
numerically-determined matrix characteristic polynomial coefficients are used extensively
in engineering linear systems analysis and control theory where empirical evidence
suggests they are effective. Thus, while it is clear that numerical calculation
of the explicit coefficients cannot be recommended as a general technique, the classes
of matrices for which these techniques break down and a complete understanding
of their failure modes is not well-established. For example, Wilkinson observes that
Jacobian matrices derived from damped oscillatory systems are frequently associated
with well-conditioned characteristic polynomials. In such circumstances, eigenvalue
methods based on explicit use of the coefficients appears to be fast and reliable, based
on numerical experiments with an early iterative solver, DEUCE [37].
To circumvent difficulties in explicitly determining the characteristic polynomial
coefficients, we seek a method to determine whether a square matrix has a pair of
eigenvalues whose sum is zero directly from the entries in the Jacobian matrix. A
simple procedure for doing this involves Kronecker or tensor products of matrices.
For the finite dimensional vector spaces V and W , let be a linear
operator with (n \Theta n) matrix representation in terms of the basis e i , and
let another linear operator with (m \Theta m) representation
terms of the basis f k . The tensor product T
V\Omega W is an (mn\Thetamn)
matrix with entries a ij b kl in terms of the basis e
V\Omega W . Moreover, eigenvalues
behave multiplicatively under tensor products: If - i and - k are eigenvalues of T 1 and
corresponding to eigenvectors u i and v k , then - i - k is an eigenvalue of T
in addition, we assume T 1 and T 2 are non-defective and dim(V then the
tensor sum (T
special action on its eigenspaces. Since
the spectrum of this operator consists of the n 2 pairwise sums (-
with the eigenvectors u
This tensor sum suggests an obvious candidate for an augmenting function: If J
is the Jacobian, then J has an eigenpair that sums to zero if and only if
det(J\Omega I n
n\Omega J)
vanishes where I n is the appropriate identity matrix. Notice that since each eigensum
of distinct eigenvalues occurs twice, the tensor sum will have corank-2 at a Hopf
bifurcation point, which has important numerical disadvantages.
To remove the twofold redundancy in the eigenvalues of the tensor product, one
can split
into the eigenspaces for \Sigma1 of the involution oe :
V\Omega V that
interchanges factors: oe(v
. The restriction to the n
dimensional
eigenspace of \Gamma1 for oe is an operator whose eigenvalues are (- i
commutes with oe, so it maps the
oe-eigenspaces to themselves. In what follows we construct the matrix corresponding
to the restriction to the (-1) eigenspace of oe. Using a different approach, the matrix
representation of this restricted operator was originally constructed directly from the
elements of the argument matrices by St'ephanos [33] and later by Fuller [12]:
Let A and B be n\Thetan matrices with entries (a ij ) and (b ij ), respectively,
1). Then the bialternate product (or biproduct) of A and
B, denoted AfiB, is an m\Thetam matrix whose rows are labeled pq for
columns labeled rs,
a pr a ps
ps
a qr a qs
oe
In the case we obtain from this an operator on an n
dimensional space
whose eigenvalues are the pairwise sums of the eigenvalues of A without repetition:
Theorem 3.1. Let A be an (n \Theta n) matrix with eigenvalues (-
(i) Afi A has eigenvalues - i
(ii) 2Afi I n has eigenvalues - i +- j ,
where I n is the n-square identity matrix and 1-j! i -n.
Notice that the subscripts used to compute the biproduct entries are indexed by pairs.
To avoid confusion, when the lexicographically ordered pair, or label, is used to refer
to an element of A fi B, we shall enclose it in braces, as in the definition above. When
the row/column matrix index is used, we employ the standard notation. Thus, the
index provides information about the position of the entry in the biproduct matrix
while the label indicates which elements of the arguments to the product are involved.
As an example, consider an arbitrary 3\Theta3 matrix 3. Then, the
conventions require that
For compactness in the index notation, when there is no possibility of confusion we
drop the comma separating the row and column indices. Thus, in the example, the
two products of A are given by:
a 22 a 21
a 12 a 11
a 13 a 11
a 23 a 22
a 13 a 12
a 12 a 11
a 13 a 11
a 13 a 12
a 22 a 21
a 23 a 21
a 23 a 22
I 3 =4 a 11 + a 22 a 23 \Gammaa 13
a a 33 a 12
\Gammaa 31 a 21 a 22 + a 33
Substituting I n into the definition of the bialternate product and solving for the
elements yields a simple formula for the entries. For the n
I n with rows, pq, and columns, rs, the entries are given by the formula,
ps if
(A) pr if r 6= p and
The algebraic properties of the bialternate product transformation confers a certain
structure upon the matrix Afi I n which may be used in the construction of Hopf
algorithms. For example, simple manipulation of the row and column labels yields
the following lemma:
Lemma 3.2. Let A be an (n \Theta n) matrix. If A is lower (upper) triangular then
I n is lower (upper) triangular. In particular, if A is diagonal so is 2Afi I n .
Proof: Suppose A is lower triangular and let (A) ij be non-zero. Notice that if
then it contributes only to the diagonal elements of 2Afi I n . Consider the case for
n. The first assignment in (9) may be written:
for k such that 1-j!k! i - n. By the row and column label convention,
l
l
l / kj
shows (2Afi I n ) fik;kjg is below the diagonal. Similar arguments hold for the second
and fourth assignments in (9); the fifth does not apply in the lower triangular case.
Thus, 2Afi I n is lower triangular. A similar argument holds for A upper triangular,
and the result for diagonal A follows immediately.Simple extension of the index counting arguments used above can be used to establish
how a non-zero off-diagonal element in the Jacobian is propagated to the product
Lemma 3.3. Suppose (A) ij is non-zero, i 6= j.
in entries of 2Afi I n . Moreover, if jg, each of these
entries is contained in a band of
i. The three assignments specified in Equation (9) which apply
may be written:
such that 1-j! i !k- n
shows the total number of off-diagonal entries in 2AfiIn contributed by (A) ij is given
by 2.
The non-zero entry (A) ij may contribute elements (2Afi I n ) fpq;rsg through each
of the three rules (13a-c) above. Let (i a
row
row
row
the row and column indices corresponding to the labels pq and rs under the rules
(13a), (13b) and (13c), respectively. To find the bandwidth we wish to maximize the
differences i
col where the asterisk is replaced by the rule identifier a,b or c.
Simple calculations show:
row
row
i a
Comparing these bounds establishes the result for (A) ij for j ! i. The case for
analogously.Finally, since the bialternate product was derived as the restriction of a tensor
product of matrices to an invariant subspace, it inherits several important properties
from tensor calculus. In particular, it is straightforward to verify:
Lemma 3.4. For ff a scalar and A;B in M depending on a parameter -, the
following properties hold:
(i) ff
(iv)
(v) @
\Theta @
@- A
\Theta @

(Provided the partial derivatives exist.)
(vi) (AB) fi I
3.2. Bialternate Product Algorithms. From the algebraic theory of symmetric
matrix products described above we have a simple necessary condition for a
Hopf point: If the point bifurcation point for -
then the
1)-dimensional system
det
vanishes. However, we have not found a condition that distinguishes purely imaginary
eigenvalues directly from the Jacobian and its bialternate products in analogy to the
subresultant criterion described earlier.
One can use bialternate products effectively in continuation method calculations
of Hopf bifurcation curves relying on the observation that transitions where eigenvalues
depart from the imaginary axis occur at degenerate bifurcations. As before, we
propose following curves of Hopf bifurcations in two parameter families by solving the
equations expressing the singularity of the bialternate product of the Jacobian, while
ensuring that the imaginary eigenvalues are bounded away from zero and that multiple
imaginary eigenpairs do not occur. We begin with the simple proof that Equation (14)
is regular on its zero set provided appropriate genericity and transversality conditions
hold:
Theorem 3.5. Suppose the system
has an equilibrium (x   ; ff   ) at which the following properties are satisfied:
single pair of eigenvalues -
whose sum is zero; ie., (-
dff
Then is an isolated non-singular solution of:
det (D x f fiI n )
Proof: Let It suffices to show that D y F (x   ; ff   ) is non-singular. By
Keller's bordering lemma [25],
is non-singular if d(x   ; ff   since, by assumption
k=3 be the complement of the critical
eigenvalues of J x (x   ; ff   ). The product
Y
Y
defines a smooth function in a neighborhood of the solution, which, by (E2), is non-zero
at . For simplicity of notation, let J
using Theorem 3.1,
@ff
since the second term in the sum vanishes. But it is also true that
dff
Differentiating the equilibrium condition f(x; evaluating the result at
dx
dff
Substituting Equation (18) into (17) we have
\GammaD x [det (J x fiI n )]
as well. Together with Equation (16) this implies
which establishes the result.To utilize either augmented system defined by the vanishing of a resultant or
biproduct in the Euler-Newton continuation framework we must specify how the Jacobian
of the extended system is to be computed since, at each corrector step, a linear
system of the form
must be solved. Here, and in the discussion that follows, we use to
denote the augmented column vector of independent variables of F . Referring to
Theorem 3.5, we will assume that the components D x f(x; ff; fi) and D ff f(x; ff; fi)
are readily available either analytically, by automatic differentiation (See [18] and
references therein) or, in the worst case, by finite-difference estimation. One way of
computing the derivatives
is to apply a forward or central difference formula to the scalar-valued function
for the (n partials. An alternative is suggested by an extension of a
lemma due to Halanai (published by Davidenko [11]):
Lemma 3.6. Let be a matrix, 1- whose entries are C 1
real-valued functions a ij
@
@x l
det
@x l
where tr() denotes the trace function and Adj() is the adjoint matrix of A. Moreover,
if A(x) is invertible
in\Omega , then
@
@x l
det
A
@x l
for all x
2\Omega .
Identifying the matrix A in the Lemma with either JfiIn or a resultant matrix
from Section 2.2 shows that the partial derivatives of the augmented equation may
be computed without resorting to differencing the determinant function. Moreover,
in those cases when second derivatives of f are known, these formulae indicate how
they may be used directly. The adjoint form is undesirable since the computation of
Adj(A(x)) is O(m 4 ); however, since the objective of the corrector step is, essentially,
to drive A to singularity we expect the calculation of A \Gamma1 to be increasingly unstable
near a solution point. Therefore, we require a similar formula valid near a generic
solution where we expect the rank of A to drop by one due to the vanishing of a
single eigenvalue. In the case that the singular values oe of A are well-
separated, the results of Chan [7] show that we can isolate the vanishing pivot in its
LU-decomposition using an appropriately chosen permutation matrix, P, to form
where ffl is on the order of the smallest singular value of A. Thus, we expect the
conditioning of L and U to be much improved compared with that of A, and using
block forms for the inverses of LA and UA together with the relation
ffl \Delta det(U), we obtain by rearrangement:
Lemma 3.7. Let A be a matrix defined as in lemma 3.6 and suppose that A is
invertible in some open
l -k and each x
2\Omega there exist matrices P, Q, U and L,
vectors v and w and ffl ? 0 such that
@
@x l
det
@x l
where P is a permutation matrix, Q is orthogonal, ffl is either O(oe n ) or O(- \Gamma1
and
Applying Lemmas 3.4 and 3.7 to the augmented system defined by Equation (14),
one obtains for 1- l - (n
@
@y l
det (J(y)fiI n
'- @
@y l
J(y)
As this formula shows, each partial of det(JfiI) is composed of weighted sums of entries
in the matrix derivatives of J. In particular, Z is the weight matrix and det(U) a scale
both dependent upon J but not on l. Thus, to compute the row of entries in
corresponding to the augmented singularity equation,

one need evaluate these factors only once per corrector step.
We conclude this section with a few remarks concerning the properties of bialter-
nate product matrices germane to the solution of the linear algebra problems which
arise in algorithms for Hopf bifurcation. Specifically, we wish to exploit the structure
of JfiIn in such a way as to mitigate the dimension increase of the linear systems
from n to n(n \Gamma 1). Foremost among these is the observation that the bialternate
product matrix of J is sparse - very sparse - even for n of modest dimension. Table
2 shows the fraction of non-zero entries in JfiIn for dense J as a function of n. The
sparsity of the bialternate product is a consequence of its relationship to the tensor
sum
I\Omega J which is composed of a dense block-diagonal band of bandwidth n,
off-diagonal blocks with (at most) n non-zero entries in each block. A
similar internal structure is inherited by the bialternate product matrix, a conclusion
which is implied by the Lemmas 3.2 and 3.3. For example, consider the fate of the
th lower subdiagonal of J under the bialternate product transformation; that is,
consider the elements of J for which
3.2, its image in JfiIn remains in the lower subtriangle. Inspection of the proof for
Lemma 3.3 shows the bound for the difference in row and column indices for elements
in the bialternate product matrix originating from (J) ij ,
is tight and monotonically increasing with i. Thus, the m th subdiagonal of J generates
a wedge or fan-like structure of non-zero entries in JfiI n which is narrow in the upper-left
corner of the bialternate product matrix and achieves its maximum width as a
result of the (n; n\Gammam) th entry of J. Finally, Equation (9) shows that while the product
matrix is not symmetric, in general, its basic sparsity pattern is. To illustrate these
observations, Figure 2 shows the sparsity pattern generated by a dense Jacobian
matrix of dimension
The observation that JfiIn is band-structured and sparse may be exploited by
Hopf path-following algorithms in a variety of ways. For example, the augmented
equation defined above depends upon the singularity of the bialternate product, a
property which is preserved under similarity transformation. Thus, Equation (14)
may be replaced by
~
det
Naturally, the most desirable choices for M will
be unitary as well. Reduction of J to Hessenberg form by chosing M a product of
Householder matrices is numerically stable and yields a bialternate product which
has block Hessenberg structure with a subdiagonal bandwidth (n \Gamma 2) and a sparsity
structure in the upper triangle as described above. If the Jacobian is reduced
to tridiagonal form, J fi I n is block tridiagonal. Matrix routines designed for these
banded structures can then be exploited to achieve a corresponding reduction
of computational work.
4. Concluding Remarks. The algorithms that we have described for computing
Hopf bifurcations are specifically designed to compute points at which there is a
simple pair of pure imaginary eigenvalues for the Jacobian of a system. We briefly
discuss the extension of the algorithms to ones which seek points of codimension two
bifurcation. There are three cases which arise from codimension two conditions on the
linearization of a vector field at an equilibrium, namely, Takens-Bogdanov bifurcation,
double Hopf bifurcation and simultaneous simple zero and pure imaginary eigenval-
ues. In each case, we seek minimally inflated systems of (n equations that locate
the codimension two bifurcation points. The case of the Takens-Bogdanov bifurcation
is easy. In terms of the characteristic polynomial of the Jacobian, one wants the constant
and linear coefficients to vanish. This condition is equivalent to the Jacobian
J having corank-1 and the square of the Jacobian J 2 having corank-2. Alternatively,
the Jacobian and the bialternate product of the Jacobian both vanish. These last
criteria are also satisfied at a point with a zero eigenvalue and a pair of imaginary
eigenvalues. Locating points of double Hopf bifurcation is a bit more complex. One
can calculate in terms of polynomial remainder sequences the presence of two pairs of
eigenvalues whose sum is zero. At such a parameter value, the two polynomials r e and
constructed from the characteristic polynomial have a common quadratic factor.
The coefficients of this polynomial can be computed as subresultants of the Sylvester
matrix of r e and r o . In order for the two pairs of roots to be imaginary, the common
quadratic factor of r e and r should have negative real roots. This is easily expressed
as inequalities on the coefficients of the common (monic) factor: it must have positive
coefficients and a positive discriminant. Bialternate product methods for computing
points of multiple Hopf bifurcation are discussed in a forthcoming paper of Govaerts,
Guckenheimer and Khibnik [16].
In a companion paper [20], we examine the issues relating to implementing Hopf
continuation using the resultant and biproduct formulations, and study their performance
on a suite of example problems of current research interest in neurobiology.
Several algorithms, including the ones based on the Bezout resultant and the bialter-
nate product, are applied to compute curves of Hopf bifurcations in a two dimensional
parameter space for the following six dimensional vector field:
\Gammag Na' 3
where,
We consider the accuracy and convergence of root finding algorithms as well as the
number of operations required to compute the curves of Hopf bifurcations. Summarizing
our findings, the number of floating point operations required for continuation
of a curve of Hopf bifurcations using an augmentation function based upon the Bezout
resultant required slightly fewer operations (664777 flops) than the use of an
algorithm described by Jepson, Griewank and Reddien [17] based upon a 2n dimensional
augmentation function (701244 flops). However, in a test the convergence
of root finding from "naturally" chosen initial values for parameter values near a point
of double Hopf bifurcation, the method based upon use of the Bezout resultant as an
augmentation function gave more consistent results. Use of the the determinant of
the biproduct as an augmentation function was substantially slower (2786163 flops),
but these biproduct calculations did not exploit the sparsity of the biproduct matrix
in computing its determinant.
Acknowledgments

This research of John Guckenheimer and Mark Myers was partially
supported by grants from the National Science Foundation and the Department
of Energy. Many of the calculations were performed with the software package DsTool
[3] developed at Cornell University. Bernd Sturmfels was partially supported by NSF
grants DMS-9201453 and DMS-9258547 (NYI), and a David and Lucile Packard Fellowship



--R

'An Efficient Algorithm for the Determination of Certain Bifurcation Points,' Journal of Computation and Applied Mathematics
'Lectures on Bifurcations in Versal Families,' Russian Mathematical Surveys
'dstool: Computer Assisted Exploration of Dynamical Systems,' AMS Notices
Polynomials and Linear Control Theory
'A Generalized Gauss-Newton Procedure for Multi-Response Parameter Estimation,' SIAM Journal of Scientific and Statistical Computing

'On the Existence and Computation of LU-Factorizations with Small Pivots,' Mathematics of Computation


'Inversion of Matrices by the Method of Variation of Parameters,' Soviet Mathematics
'The Evaluation of Determinants by the Method of Variation of Parameters,' Soviet Mathematics
'Conditions for a Matrix to Have Only Characteristic Roots With Negative Real Parts,' Journal of Mathematical Analysis and Applications
'Two Methods for the Numerical Detection of Hopf Bifurcations,' International Series of Numerical Mathematics
Matrix Computations
'Classification and Unfoldings of Degenerate Hopf Bifurca- tions,' Journal of Differential Equations
in preparation.
'The Calculation of Hopf Points by a Direct Method,' IMA Journal of Numerical Analysis
Automatic Differentiation of Algorithms: Theory
Dynamical Systems
'Computing Hopf Bifurcations II: Three Examples From Neurophysiology.' Preprint
'On a Four Parameter Family of Planar Vector Fields,' Archive for Rational Mechanics and Analysis
'New Algorithms for the Evaluation of Complex Bifurcation Points in Ordinary Differential Equations. A comparative Numerical Study,' Applied Mathematics and Computation

'On Eigenvalues of Matrices Dependent Upon a Parameter,' Numerische Math- ematik
'Numerical Solution of Bifurcation and Nonlinear Eigenvalue Problems,' in Applications of Bifurcations Theory
'LINLBF: A Program for Continuation and Bifurcation Analysis of Equilibria Up to Codimension Three,' Continuation and Bifurcation: Numerical Techniques and Applications
'Generalized Polynomial Remainder Sequences,' in Computer Algebra - Symbolic and Algebraic Computation
'A Comparison of Methods for Determining Turning Points of Nonlinear Equa- tions,' Computing
'The Numerical Detection of Hopf Bifurcation Points,' Continuation and Bifurcation: Numerical Techniques and Applications
Bifurcation Theory and Computer Algebra
'A Direct Method for the Computation of Hopf Bifurcation Points,' SIAM Journal of Applied Mathematics
'An Algorithm for the Computation of Hopf Bifurcation Points in Comparison With Other Methods,' Journal of Computational and Applied Mathematics


'Computation of Hopf Branches Bifurcating From Takens- Bogdanov Points for Problems With Symmetries,' International Series of Numerical Math- ematics
The Algebraic Eigenvalue Problem
'The Evaluation of the Zeros of Ill-Conditioned Polynomials. Part II,' Numerische Mathematik
--TR

--CTR
W. Govaerts , Yu. A. Kuznetsov , B. Sijnave, Implementation of Hopf and double-Hopf continuation using bordering methods, ACM Transactions on Mathematical Software (TOMS), v.24 n.4, p.418-436, Dec. 1998
Randall D. Beer, Parameter space structure of continuous-time recurrent neural networks, Neural Computation, v.18 n.12, p.3009-3051, December 2006
bialternate product;hopf bifurcation;resultant