The Numerical Computation of Homoclinic Orbits for Maps.
Transversal homoclinic orbits of maps are known to generate shift dynamics on a set with Cantor-like structure. In this paper a numerical method is developed for computation of the corresponding homoclinic orbits. They are approximated by finite-orbit segments subject to asymptotic boundary conditions. We provide a detailed error analysis including a shadowing-type result by which one can infer the existence of a transversal homoclinic orbit from a finite segment. This approach is applied to several examples. In some of them parameters appear and closed loops of homoclinic orbits are found by a path-following algorithm.
Introduction
One of the fundamental results on chaotic behavior in discrete dynamical systems is Smale's
Homoclinic Theorem, see [23], [21], [22]. For a more recent overview of homoclinic orbits,
their bifurcations and the history of their discovery we refer to [19].
Consider a (time-)discrete system
with a C l -diffeomorphism f assume that - 2 R k is a hyperbolic fixed point
of f , i.e. the Jacobian f 0 (-) has no eigenvalues on the unit circle. Further assume that x 0 is
a transversal homoclinic point, which means that the orbit generated by (1.1) satisfies
lim
and the stable and unstable manifolds of - intersect transversally at x 0 . Then the theorem
states that there exists a compact set M and an integer p such that the p-th iterate of f ,
denoted by f p , leaves M invariant and is topologically conjugate on M to the Bernoulli shift
on two symbols. It is remarkable that this chaotic behavior on a certain subset is created
by a homoclinic point that is transversal, a property that persists under the perturbation of
system (1.1).
The perturbation stability of transversal homoclinic points suggests that we should be
able to compute these points numerically in a robust and stable way. This is the topic of the
Supported by Sonderforschungsbereich 343 "Diskrete Strukturen in der Mathematik", Fakult?t f?r
Mathematik, Universit?t Bielefeld.
y Parts of this paper are based on the diploma thesis [17] of the second author.
present paper. Instead of homoclinic points our approach aims at computing the complete
homoclinic orbit by solving the 'boundary value problem' (1.1), (1.2). Any
direct method for the single homoclinic point x 0 implicitly tries to solve this boundary value
problem and hence is prone to the usual difficulties of shooting-type methods which are
caused by exponential divergence of trajectories.
Therefore, we propose to approximate the infinite orbit by a finite orbit segment x
which satisfies a 'finite boundary value problem'
Here (1.4) is a general set of boundary conditions defined by a smooth mapping b
R k . Together, (1.3) and (1.4) comprise a set of (n+ equations for the
same number of unknowns. We solve this system by Newton's method and take advantage
of the sparsity pattern of the corresponding Jacobian.
The most important examples for (1.4) are periodic boundary conditions
and projection boundary conditions
where the zero sets of b \Gamma and b + are linear approximations to the local unstable and stable
manifold of the fixed point.
The whole approach closely mimics the well known methods for approximating homoclinic
orbits of parametrized continuous systems -
on finite intervals [9]. [4], [2],
[11]. But there is one important difference to keep in mind. In the discrete case, stable and
unstable manifolds can, and in the generic sense do, intersect transversally at a homoclinic
point while in the continuous case both tangent spaces always contain the flow direction. If
the manifolds intersect at all in a continuous system they do so in a whole curve obtained
by shifting the phase of the orbit. Therefore, one parameter is needed for homoclinic orbits
to occur generically in a continuous system and an additional constraint is needed to fix the
arbitrary phase.
The paper is organized as follows. In section 2 we collect the necessary prerequisites from
the theory of exponential dichotomies for difference equations. This section is largely based
on the paper of Palmer [21] where exponential dichotomies are used to prove the Shadowing
Lemma and Smale's Homoclinic Theorem. Earlier references on dichotomies for difference
equations include [5] and [15].
Section 3 contains the basic results of this paper. In a neighborhood of a transversal
homoclinic orbit we find a unique solution of the equations (1.3), (1.4) provided b satisfies
a nondegeneracy condition and are taken sufficiently large. We also show error
estimates which are analogous to the continuous case.
In section 4 we treat several examples including the Poincar'e map for the Duffing system.
We also compute branches of discrete homoclinic orbits in cases where the system (1.1)
contains a parameter. Using a standard continuation procedure (see [1]) we can easily pass
through turning points of these branches. At such turning points the homoclinic orbits
actually become nontransversal.
In section 5 we take up the question of converse theorems which, given a solution of
the finite boundary value problem (1.3), (1.4), guarantee the existence of a transversal
homoclinic orbit. This is in the spirit of the recent shadowing type results as in [6], [7]
except that here we conclude the existence of an infinite sequence, rather than a continuous
orbit, from the knowledge of a finite orbit.
We emphasize that we do not intend a computational verification of the existence of a
homoclinic orbit as in [18] where interval analysis is used to verify an exact homoclinic orbit
by a pseudo orbit.
Exponential dichotomies for difference equations
In this section we review some basic tools from the theory of exponential dichotomies [21]
and we add a few results which are useful later on. Consider a homogeneous difference
equation in R k
with solution operator
I if
A
In the following let
be any interval in Z. If no confusion with real intervals arises we simply write
We also make frequent use of the Banach space of bounded sequences on J given by
Definition 2.1 The difference equation (2.1) has an exponential dichotomy on J if
there exist projectors Pn ; (n 2 J) in R k and constants K; ff ? 0 such that
and
For brevity we will say that (2.1) has an exponential dichotomy on J with data (K; ff; Pn ).
It is easy to see that this implies that the adjoint equation
has an exponential dichotomy on
where C is a constant such that
Notice that the solution operator of (2.2) is given by
\Phi   (n; m) := \Phi(m
On semifinite intervals the ranges of the projectors Pn are uniquely determined
and we can solve inhomogeneous equations of the type
according to the following lemma (cf. [21, Proposition 2.3, Lemma 2.7]).
Lemma 2.2 Let (2.1) have an exponential dichotomy on data (K; ff; Pn ).
Then for any n - n \Gamma we have
Suppose that (L; fi; Qn ) is another set of dichotomy data on J . Then R(Qn
and we have the estimate
For any - 2 R(Pn \Gamma ) and any bounded sequence r J 2 S J the system (2.3) has a unique
bounded solution u
Similarly, in the case
and the inhomogeneous equation
has a unique bounded solution which again satisfies (2.4).
In the next step we consider the case merely assume that (2.1) has dichotomies
on separately. The following lemma is the discrete analogue of [20, Lemma 4.2].
Lemma 2.3 Let have an exponential dichotomy on Z \Gamma with data
and on Z + with data (K
assume that An and A \Gamma1
are uniformly bounded.
Then the operator
is Fredholm of index rank(P
Let L   : SZ ! SZ be defined by
then dimN (L
Proof. From Lemma 2.2 it follows that any - 2 R(P
defines an element
and that this map is an isomorphism, hence
In a similar way we use the exponential dichotomy of L   on (\Gamma1; \Gamma1] and on [\Gamma1; 1) with
the projectors I
An application of Lemma 2.2 with
is an isomorphism from
onto N (L   ). Therefore
dimN (L
For the proof of (2.7) we first notice that
exists for all 'Z 2 N (L   ), r Z 2 SZ
since 'n decays exponentially as n ! \Sigma1. Assuming
For a proof of the converse in (2.7) we let u
be the solutions of the
initial value problem (2.3) and(2.5) with
respectively.
Next we show
We conclude that
is an element of N (L   ) and hence by our
assumptions
Thus (2.9) holds and there exists some - 2 R k such that
From this relation we find that the sequence
is consistently defined at Therefore r Z 2 R(L) holds.
From (2.7) we obtain dimN (L   be an orthogonal
basis of R(P
we know that ' i a basis of
N (L   ) and it is easy to verify with the help of (2.7) that the sequences u
form a basis in SZ =R(L).
Finally, we can compute the Fredholm index of L from (2.6) and (2.8),
In the special case rank(P +
0 ) the operator L is Fredholm of index 0. If in
addition, holds then for any r Z 2 SZ the inhomogeneous system
has a unique solution uZ 2 SZ . More precisely by [21, Proposition 2.6, Lemma 2.7] we have
Lemma 2.4 The difference equation (2.1) has an exponential dichotomy on Zif and only if
it has one on Z \Gamma and Z separately with projectors of equal rank and if (2.1) has no bounded
nontrivial solution. In case one of these equivalent assertions holds then the system (2.10)
has a unique solution uZ 2 SZ for each r Z 2 SZ which satisfies
Here (K; ff; Pn ) are the dichotomy data on Z.
Finally, we need the following consequence of the roughness theorem, see [21, Proposition
2.10, Remark 2.11].
Proposition 2.5 Consider two difference equations
with nonsingular matrices An ; Bn 2 R k;k such that
Assume that Anun has an exponential dichotomy on Z + with data (K; ff; Pn ) and
ff. Then the system Bnun has an exponential dichotomy on Z + with
suitable data (L; fi; Qn ) where Qn has the same rank as Pn and
Proof. Because of (2.11) we can choose N 0 such that the roughness theorem [21, Proposition
2.10, Remark 2.11] yields an exponential dichotomy for Bnun on any interval
well as an estimate
-N
We notice that the constants C 1 , K 1 and fi are independent of N , and hence the projectors
Pn and Q (N)
n are of equal rank if N 0 is taken sufficiently large. For a given " ? 0 we may
choose
An application of Lemma 2.2 with
Therefore, we can select N 2 - N 1 such that
By the triangle inequality we obtain
which proves (2.12). Finally, the exponential dichotomy on [N easily carries over to
after adjusting the constant K 1 (cf. [21, Remark 2.2]). \Xi
3 Approximation of transversal homoclinic orbits
k be a C 1 -diffeomorphism with a fixed point - 2 R k . We call
a homoclinic orbit with respect to - if the following holds
lim
We may write (3.1) as an operator equation
where the C 1 defined by
In our definition we have included the trivial case N. All the results
of the following sections also hold for this trivial case. Of course, in the applications we are
only interested in nontrivial orbits.
Assume that - is a hyperbolic fixed point and consider - -u such that each
stable eigenvalue - of f 0 (-) satisfies and each unstable one satisfies j-
be the corresponding decomposition into the stable and unstable subspaces
of f 0 (-) and let us choose the norm in R k such that the following holds (see [16])
This implies that has an exponential dichotomy on any interval J ae Zwith
solution operator
and data (1; ff; P s ) where
and P s is the projector onto X s along Xu .
The local stable and unstable manifolds of - are defined as
loc
W u
loc
where V ae R k is a sufficiently small neighborhood of 0. By definition these sets are positive
and negative invariant respectively. It is well known (see [16, Chapter 6]) that we may take
of the origin such that the local
stable and unstable manifold of - may be represented as graphs. We get
where the functions q
The global stable and unstable manifolds are then given by
loc
loc
with the differential structure induced by W s;u
loc via the mapping f .
Theorem 3.1 Let xZ be a homoclinic orbit of (3.1) with respect to a hyperbolic fixed point
-. Then the following assertions are equivalent.
(i) x 0 is a transversal homoclinic point, i.e. the tangent spaces to the stable and unstable
manifolds at x 0 satisfy
(ii) The linear system
has an exponential dichotomy on Z.
(iii) xZ is a regular zero of the
is a homeomorphism.
Proof. In [21, Definition 4.2] condition (ii) is used to define a transversal homoclinic point
and the equivalence with (3.8) is proved later on [21, Proposition 5.6]. The implication
(ii))(iii) is immediate from Lemma 2.4 since
Let us assume (iii). Then we have
and Proposition 2.5 shows that (3.9) has exponential dichotomies on both Z \Gamma and Z + with
projectors of the same rank. By our assumption equation (3.9) has no nontrivial bounded
solution on Z, so Lemma 2.4 proves assertion (ii). \Xi
Definition 3.2 A homoclinic orbit with respect to a hyperbolic fixed point - is
called transversal iff it satisfies one of the equivalent condition of Theorem 3.1.
Next we consider the approximate equations (1.3), (1.4) on a finite interval
and write them as
defined by
We treat equation (3.10) with the help of the following approximate inverse function
theorem (cf. [24, x3], [4, Lemma 3.1], [17, Lemma 4.7])
Proposition 3.3 Let F 2 C 1 (Y; Z) with Banach spaces Y , Z and assume that F 0 (y 0 ) is
a homeomorphism for some y 0 2 Y . Further, let constants ae; -; oe ? 0 be given with the
following properties
has a unique solution in B ae (y 0 ) and the following estimates hold
We can now state our main approximation theorem
Theorem 3.4 Let - xZ be a transversal homoclinic orbit with respect to a hyperbolic fixed
point - of a diffeomorphism f 2 C 1 (R k
and the map defined by
is nonsingular.
Then there exist constants ffi; C such that (3.10) has a unique solution
x J in
for all and the following estimates hold
Proof. We apply Proposition 3.3 with
We consider first
xn+ )un+ and show
for n \Sigma sufficiently large. This will be done in two steps.
Step 1: For any sequence z n , n 2 ~
exist un , n 2 J such that
and ku J k1 - Ckz ~
J k1 .
has a unique solution v J 2 S J for each r 2 R k and kv J k1 -
Ckrk.
Here C ? 0 denotes some generic constant independent of n \Sigma and the right hand sides. The
constants are chosen within the proof.
Suppose we have accomplished these two steps. Then for given z ~
J and r 2 R k we choose
u J as in step 1 and let v J solve
as in step 2. This implies
as well as
Thus we have shown the estimate (3.18).
For y sufficiently small we find
Taking
2C and ae = ffi in (3.11) we finally obtain from assumption (3.14)
for n \Sigma sufficiently large. Proposition 3.3 then yields the existence of the solution x J in
and the estimate (3.17) follows from (3.12) by setting y
Proof of step 1. Let (K; ff; Pn ) be the dichotomy data associated by Theorem 3.1 with
the linear system
We extend z ~
J by setting z
J . Then we use Theorem 3.1(iii) and solve
Obviously, un (n 2 J) solves (3.19) and satisfies
Proof of step 2. We notice that Proposition 2.5 implies
For jnj large we have kPn \GammaP s k ! 1, hence the matrices
are nonsingular with
and we find R
is bijective and satisfies
Similarly, bijective and
Using the solution operator \Phi(n; m) for equation (3.9) we may write v J as
are to be determined from the boundary conditions
By the exponential dichotomy we have
It is convenient to change coordinates in (3.24) via
We employ the linear map B from (3.15) and rewrite (3.24) as
where
Using (3.20), (3.21), (3.22) and (3.25), (3.26) we obtain
Therefore, assumption (3.15) implies that the linear equation (3.27) has a unique solution
with a constant C independent of n \Sigma . We then define by (3.26) and vn by (3.23).
From (3.20) and the exponential dichotomies we obtain the estimate of step 2. \Xi
Remark 3.5 Assuming that k\Gamma 0
moderate, the constant C 1 in (3.16) is the crucial
quantity which measures the condition of the nonlinear system (3.10). A more quantitative
analysis (see [17]) shows that C 1 depends linearly on 1
is the angle between the
stable and unstable manifold at -
x 0 . Therefore, this constant may also be taken as a measure
of the transversality of the orbit.
Inequality (3.17) shows that the approximation error k-x k1 is determined by the
boundary error
We say that the boundary conditions (1.4) are of order (p there exists a constant C
such that for all \Gamman \Gamma , n+ sufficiently large
In the hyperbolic case we have
are the spectral bounds for f 0 (-) as in the beginning of this section.
From (3.17), (3.28) and (3.29) we then obtain the estimate
Clearly, the periodic boundary conditions (1.5) are of order (1; 1) since
Moreover, they always satisfy conditions (3.14) and (3.15).
Approximations of second order are obtained from projection boundary conditions
given by
s
where the columns of Q s 2 R k;ks and Qu 2 R k;ku , k s a basis of the stable
and unstable subspace of f 0 (-) T respectively. More formally, we have a block diagonalization
s
!/
s
eigenvalues - of / L s 2 R ks;ks and 1
- of / L u 2 R ku ;ku . Transforming (3.31) into
where
s
we find for the projectors
Therefore, the projection boundary conditions require
which means that xn \Gamma , xn+ lie in the linear approximation of the stable and unstable
manifold at -.
The second order easily follows from (3.7). For f 2 C l (R k ; R k ) the functions q s , q u in
are C l -smooth. Hence for l - 2 we find for x 2 W u
lox by a Taylor expansion
and a similar estimate holds for x 2 W s
loc .
We notice that the projection boundary conditions also satisfy (3.14), (3.15) since
is nonsingular.
Corollary 3.6 Let -
xZ be a transversal homoclinic orbit as in Theorem 3.4. Then the conclusions
of Theorem 3.4 hold for periodic as well as for projection boundary conditions.
The corresponding finite orbits x -
For the second estimate we assume f 2 C 2 .
If, in the periodic case , we take n
N ) is a period-2 orbit of f N . It is
contained in the invariant set M on which f N exhibits shift dynamics according to Smale's
Theorem.
In conclusion of this section we outline the extension of our results to orbits that connect
two fixed points of a parametrized system
where f(\Delta; -) is a C 1 -diffeomorphism for all - 2 R p .
Assume that (-x Z ; -
R p is a solution of (3.32) such that
where - \Sigma are hyperbolic fixed points of f(\Delta; -
-) with stable and unstable dimensions k \Sigmas ,
. Introducing the operator
we find
where
By Proposition 2.5 we have exponential dichotomies for L on Z \Sigma with data (K \Sigma ; ff
and rankP \Sigma
2.3 the operator is Fredholm of index k+s \Gamma k \Gammas in the space
SZ and hence by the bordering lemma
-) is Fredholm of index
We want (-x Z ; -
-) to be a regular zero of \Gamma which leads us to the following definition.
Definition 3.7 The connecting orbit (-x Z ; -
-) is called transversal, if the following conditions
hold
Condition (i) yields that \Gamma 0 (-x Z ; -) is Fredholm of index 0 and (ii) implies that its null
space is trivial. We notice that this definition is in accordance with the homoclinic case
Assuming (i) we may express (ii) in Definition 3.7 equivalently by the following two
conditions
the matrix
is nonsingular where the columns of \Phi n 2 R k;p , (n 2 Z) form a basis of N (L   ).
It is clear that (ii) implies N hence (3.34) by (2.6). From (i) and Lemma 2.3
we then find
dimN (L
and (3.35) is a consequence of (ii) and (2.7). Conversely, assume (3.34), (3.35) and let
by (3.35). Therefore, xZ 2 N (L) holds and (3.34), (2.6) yield
For the extension of Theorem 3.4 to the general case (3.32), (3.33) we consider the
approximate system
Theorem 3.8 Let (-x Z ; -
-) be a transversal connecting orbit of the parametrized system (3.32)
with hyperbolic endpoints - \Sigma . Assume b 2 C 1 (R 2k+p
and the map
defined by
is nonsingular. Here
s is the stable subspace of f x
u is the unstable subspace
of f x (-
-).
Then there exist constants ffi , C ? 0 and N 2 N such that (3.36) has a unique solution
-)) for all and the following estimate holds
-)k:
The proof is quite similar to that of Theorem 3.4 and will be omitted. We notice that
B can be represented by a quadratic matrix, since by the transversality condition (i)
The projection boundary condition (3.30) now depends on the parameter -
Here - \Sigma (-) are hyperbolic fixed points of f(\Delta; -) and the columns of Q s (-) 2 R k;k \Gammas , Qu (-) 2
R k;k+u form bases of the stable subspace of f x (- \Gamma (-) T and the unstable subspace of
f x (-) T and all the quantities depend smoothly on -.
So far we considered the parameters - in (3.32), (3.33) to be part of the unknowns which
have to be determined by the numerical calculation, that is by Newton's method applied to
(3.36). If there are more than k \Gammas \Gamma k+s parameters in the system then we can use these
for the continuation of connecting orbits of type (k \Gammas ; k+s ). This will be done in the next
section for the continuation of homoclinic orbits.
4 Numerical implementation and applications
In the following examples we compute approximate homoclinic orbits by applying Newton's
method to the nonlinear system
This is justified by Theorem 3.4.
An iteration vector y J 2 S J for (4.1) is accepted if for a given tolerance " ? 0 the
inequality
holds. In each Newton step we use LU decomposition with absolute column pivoting and
take advantage of the sparsity pattern of \Gamma 0
J (y J ). With exception of Example 2 (Duffing
map) we use . The machine precision used in all computations is about 10 \Gamma16 .
We mainly use projection boundary conditions as in (3.30) except for the first example
where we also test the periodic boundary conditions (1.5).
In all examples to follow we have no explicitly given transversal homoclinic orbit. Actu-
ally, we do not know of any suitable example with explicit homoclinic orbits and an analytic
map f . Therefore, we take a zero of \Gamma [-n very large values of \Gamma -
orbit'. For the error analysis we then compare this orbit to various shorter orbits obtained
for moderate values of n \Gamma , n+ . We always use the Euclidean norm in R k .
-0.4
y
x
loc
W u
loc
x \Gamma2
x \Gamma3

Figure

1: Parts of the stable and unstable manifold for H'enon's map together with points
of the 'exact orbit'.
Example 1 (H'enon map)
For b 6= 0 the map
is a C 1 diffeomorphism and it has two fixed points for (b
In the following we take choose the parameters
We get fixed point and
as eigenvalues of f 0 (-). We set - 1:92373. For more details about the
H'enon map see [14, x15.2], [8, x2.9].

Figure

shows parts of the stable and unstable manifolds of -. They are obtained by
starting forward and backward iterations on the linear approximation of the unstable and
stable manifolds at - in some neighborhood of -. We show the orbit segment -
x \Gamma100;100 that
was taken as 'exact orbit' for comparison. It is calculated using the starting vector v [\Gamma100;100]
with
global
logarithmic
error
projection boundary conditions
periodic boundary conditions

Figure

2: Global error e 1 (\Gamman; n) with n varied using periodic and projection boundary
conditions for H'enon's map.
To analyze the approximation error we compute zeros y starting
ton's methods with the truncated 'exact orbit' - x . For these solutions let
be the local logarithmic error at
be the global logarithmic error.
In

Figure

2 we compare the global error e 1 (\Gamman; n) obtained for the projection and the
periodic boundary conditions. For the projection boundary conditions we get a slope of
nearly while the slope for periodic boundary conditions is nearly log
is predicted by Corollary 3.6. For the remaining tests we always use projection boundary
conditions.

Figure

3 shows the behavior of the local error along the approximate homoclinic orbits.
The right boundary is fixed and the left boundary is varied. We notice that the maximal
error occurs at the boundaries. A more detailed study of the local error at the boundaries
is given in Figure 4. The local error e(n; n; 4) at the left boundary has nearly the expected
slope of while the local error e(4; n; 4) at the right boundary is nearly constant.
So the global error cannot pass below this value e(4; n; 4). The dependence of the global
error on both variables n \Gamma and n+ is illustrated in Figure 5.
i scaled to [0,1]

Figure

3: Local error e(i; n; 4) with varied and i scaled to [0; 1] for H'enon's
map.
logarithmic
error
global error
local error e(n,n,4)
local error e(4,n,4)

Figure

4: Global error and local errors at the boundaries varying the left boundary for
H'enon's map.
logarithmic error

Figure

5: Global error with varied left and right boundary for H'enon's map.
Example 2 (Duffing map)
We consider the forced Duffing equation
with ! 6= 0 and define f as the period map, that is, f(x; y) is the value of the solution of
(4.4) at
with initial values (u(0); Because of the smoothness of (4.4)
the map f is a C 1 -diffeomorphism. We choose the parameters as
The initial value problems are solved with a code called 'ode' from NETLIB. This code is
an implementation of an Adams' method by Gordon and Shampine and its tolerance is set
to 10 \Gamma15 . The Jacobian of f is calculated by solving the variational equation of (4.4) and
the fixed point
of f is approximated using Newton's method with tolerance 10 \Gamma14 . For the eigenvalues of
we get the bounds
For more details about the Duffing map see [13, x2.2].
Due to the extreme values of the Floquet multipliers - s -u it is sufficient to take -
x [\Gamma15;15]
y
x
loc
W u
loc

Figure

Parts of the stable and unstable manifold together with points of the 'exact orbit'
for the Duffing map.
as 'exact orbit'. To compute it, Newton's method is started at v [\Gamma15;15] with
We set the tolerance to 10 \Gamma13 and use projection boundary conditions.
In

Figure

6 we show parts of the stable and unstable manifolds of - together with the
'exact orbit' -
x [\Gamma15;15] . The behavior of the approximation errors was quite similar to H'enon's

Figures

so we do not display them here.
We investigate for this example how the error depends on the perturbation of the boundary
condition. Let Q s , Q u and b + be as in (3.30). We set
s
cos ' sin '
sin ' cos '
and use the 'rotated' boundary condition
Calculating zeros of \Gamma [\Gamma2;5] for various values of ' gives the global errors shown in Figure 7.
We observe a sharp peak at clearly shows that the higher order convergence
from Section 3 only occurs at the accurate projection angle.
global
logarithmic
error

Figure

7: Global error e 1 (\Gamma2; 5) using a perturbed left projection boundary condition for
the Duffing map.
Example 3 (3D-H'enon map)
Changing variables x 7! ~
x
a , y 7! b~y
a for a 6= 0 H'enon's map is transformed to ~
f (~x; ~
x), see [14, Examples 15.11]. We define a 3D-H'enon map as
As the original H'enon map it is a C 1 -diffeomorphism for b 6= 0 and has two fixed points for
(b
We use the same parameter values as above, see (4.2), and consider the fixed point
This fixed point is hyperbolic with eigenvalues
so we use
For these calculations we use projection boundary conditions and the 'exact orbit' is computed
by starting at v [\Gamma200;200] defined as in (4.3). Some numbers are given in Table 1.
In

Figure

8 we show the global error and the local errors at the boundaries when both
boundaries are varied, like Figure 4 for the original H'enon. We notice that the slope of
e(\Gamman; \Gamman; n) is while the slope of e(n; \Gamman; n) is oscillating around
\Gamma6 0:807797 0:928589 0:856507
\Gamma5 1:004416 0:807797 0:928589
\Gamma4 0:669725 1:004416 0:807797
\Gamma3 1:193808 0:669725 1:004416
\Gamma2 0:276148 1:193808 0:669725
3 0:972826 0:507189 1:161985
5 0:908624 0:802205 0:972826
6 0:866250 0:908624 0:802205

Table

1: Some numbers of the 'exact orbit' (-x
-22
logarithmic
error
global error
local error e(-n,-n,n)
local error e(n,-n,n)

Figure

8: Local errors for the approximate homoclinic orbits, obtained with projection
boundary conditions for a 3D-H'enon map. The error behavior for the 3D-H'enon map varying
both boundaries.
amplitude

Figure

9: Following a branch for the H'enon map.
Continuation of homoclinic orbits
Let us assume that the map f in Theorem 3.1 depends smoothly on a parameter - 2 R (see
also
If this system has a transversal homoclinic orbit at some value of - then it is easy to show by
Theorem 3.1 and the implicit function theorem that a local branch of transversal homoclinic
orbits passes through the given one. Moreover, by an extension of Theorem 3.4 one can
show that there is a corresponding branch of approximating finite orbit segments.
To follow such branches numerically we employ a predictor-corrector method from [1,
Algorithm 3.3.7] and take care of the sparse matrix structure, see [1, Chapter 10]. The
corrector is of Gauss-Newton type using Moore-Penrose inverses. Its tolerance is set to
In the following let
be the amplitude of an orbit segment y
As a first example we revisit H'enon's map. We take the same parameters and starting
vector for the first corrector step as in Example 1. The parameter b is varied and projection
boundary conditions are used. Figure 9 shows the result of the numerical computation. In
an amplitude versus parameter diagram we get a closed loop of homoclinic orbits. A closer
look at the numerical values shows that each point of the homoclinic orbit describes a closed
loop in phase space upon variation of b.
scaled
amplitude

Figure

10: Closed loop of homoclinic orbits for Euler's map. The amplitude is scaled to
(Ampl(y
A remarkable feature is that nontransversal homoclinic orbits occur at the turning points
of the homoclinic loop. Due to the standard reparametrization during the continuation
procedure there are no problems to pass through these points. We notice however, that the
approximation properties of the finite orbits close to these turning points are not covered
by Theorem 3.4.
Example 4 (Euler map)
Consider the continuous system
Let f(x; be the result of one Euler step for (4.6) with steplength h starting at (x; y),
that is
In [3]the question was raised, how such a map behaves when the continuous system passes
through a homoclinic bifurcation. For the continuous system (4.6) there exists a homoclinic
orbit based at
The above question was largely solved in [10]. There it is shown that transversal homoclinic
orbits for the one-step map can occur for all sufficiently small step sizes h and they do
so at most within a wedge in (h; -plane that is exponentially small. We will explore this
wedge to some extent.
In the following numerical experiment we fix 0:4 and we start the branch of homoclinic
orbits for the Euler map (4.7) at the value suggested by [3, x2]. For
estimated
norm
of
the
inverse

Figure

11: The estimated norms of the inverses of the linearized operator obtained by fixing
the parameter - for the Euler map.
the numerical results shown in Figure 10 we used projection boundary conditions and orbit
segments with \Gamman the homoclinic orbits occur only in a very small
-interval we used
for the horizontal axis.
Here, as with the H'enon map we find a closed loop of homoclinic orbits. The two turning
points mark the boundary of the above mentioned wedge. Although the homoclinic orbits
have extremely weak transversality the continuation method has no problems since - is
never taken for parametrization.
In order to demonstrate the weak transversality we estimate (following [12, x3.5.4]) at
each point of the branch the norms of the inverses of the linearized operator obtained by
fixing the parameter -, that is
These estimates are shown in Figure 11. The norms are quite large and, as we expect, they
grow towards the turning points. This indicates the small angle between stable and unstable
manifolds, see Remark 3.5. In contrast to this, the norms of the inverses of the matrices
used to calculate the Moore-Penrose inverses which occur during the continuation are quite
moderate (- 20).
A remarkable feature of the current loop is that after one turning round with the continuation
procedure we get the same orbit in phase space but with the points y n being shifted
to y n+1 . In fact this is only true in a strict sense for the infinite orbit xZ of (4.5) while there
-5
x
y

Figure

12: The homoclinic crown: A curve followed by one point of a homoclinic orbit for the
Euler map after several turns of the closed loop. The marks indicate the initial homoclinic
orbit.
must be a deviation for the finite orbits obtained numerically. However, this deviation is so
small that it does not show up in the graphics in Figure 10.
A consequence of the shift after one turn is that all points of the homoclinic orbits lie
on a common curve in (x; This curve is shown in Figure 12. It was obtained by
plotting the points y \Gamma40 during several turns of the continuation method. If - is between the
two turning points any horizontal cut through this curve gives two orbits which exist for the
particular parameter value. For the sake of visualization in Figure 12 we have deactivated
the stepsize adaption during the continuation and we have drawn vertical lines at every 40th
point.
We finally remark that the same phenomenon of a closed loop with shifted orbits occurs
for H'enon's map when we choose parameter values different from those in (4.2), for example
5 A shadowing result for discrete homoclinic orbits
In this section we reverse the question considered in the theorems of section 3. We assume
that finite orbit segments which satisfy (3.10) are given and we want to conclude the existence
of a transversal homoclinic orbit.
One such result has been derived in [17, x4.4]. It assumes a sequence of solutions x J
which converges uniformly as J grows such that the linearizations about x J have uniformly
bounded inverses, see (3.16). In this way Theorem 3.4 can be turned into an equivalence
result.
Our aim here is to avoid these somewhat restrictive assumptions on an infinite sequence of
orbits and try to work with a single finite orbit instead. For that purpose we specialize (3.10)
to have separated boundary conditions as in (1.6)
We assume that ks are the dimensions
of the stable and unstable subspace X s and Xu of f 0 (-) and - is a hyperbolic fixed point of
the C 1 -diffeomorphism f . Moreover we assume
(b 0
(-) jXu are nonsingular. (5.4)
Assumptions (5.3) and (5.4) imply that the boundary conditions are at least of order (2; 2).
In particular, (5.2) - (5.4) are satisfied for the projection boundary conditions (3.30).
In order to avoid interrupting the argument we begin with a perturbation lemma for the
linear case which will be an important tool in the proof of our main result.
Lemma 5.1 Let A 2 R k;k be hyperbolic with stable and unstable projectors P s and P
and the norm k k be chosen so that (3.2), (3.3), (3.4) and (3.5) hold
with A instead of f 0 (-). Furthermore, let ae ? 0 and matrices An , n - \Gamma1 be given with
Then any bounded solution xZ \Gamma of
satisfies
and
Proof. By Lemma 2.2 with the equation
has a unique bounded solution y we have the estimate
The sequence z is bounded and satisfies z
Combining this with (5.6) gives the desired estimates. \Xi
Theorem 5.2 Let the above assumptions hold. Then for any constant C ? 0 there exists an
with the following conclusions.
is a solution of (5.1) which satisfies
k-x J k1 - C; (5.7)
then there exists a unique transversal homoclinic orbit xZ of f with respect to - satisfying
Remark 5.3 1. There are no explicit assumptions made on \Gamman \Gamma , n+ being large but, of
course, the inequality (5.9) requires this in practise. Moreover the homoclinic orbit xZ
may be trivial unless we know that k-x holds for at least one n 2 J .
2. As the proof will show, we can weaken the assumption \Gamma J (-x J
In this sense the theorem gives a shadowing result (see [21, Theorem 3.5]). It is of a
special type, however, since we shadow a finite orbit by an infinite one.
3. The assumptions (5.3), (5.4) will be crucial in our proof. It seems likely that condition
(5.3) can be dropped and more general boundary conditions can be allowed as in Theorem
3.4 and its converse in [17]. But we don't know of any proof or counterexample.
Proof. We assume C ? 0 to be given as in the theorem and we will collect the various
conditions on " during the proof. Let be such
that the local stable and unstable manifolds of - in the ball
have a representation as in (3.6) with
the projector onto X s along Xu and . By the choice of norms we have kP s k,
kPuk - 1.
Then an element x belongs to W u
loc iff it satisfies the boundary condition
Because of (5.4) this is equivalent to
With (5.3) we have
and (3.7) implies
Similarly, x 2 W s
loc is characterized by
and
holds.
We construct the homoclinic orbit by solving the 'boundary value problem'
~
where ~
defined by
with We apply Proposition 3.3 with
norm k k1 ,
". We have
We choose " 0 such that "
and such that the following estimates hold for all x 2
where C is the constant in (5.8). For by using (5.9)
and thus by (5.8) and the Banach Lemma
For use (5.9), (5.15) and find the estimate
Finally, we obtain
if
Notice that by (5.15) and (5.9)
From Proposition 3.3 we conclude the existence of a unique solution y J of (5.14) in B " (-x J ).
Since
the extended orbit yZ obtained by iterating
to -. Moreover, any homoclinic orbit xZ with the property (5.10) satisfies kxn
and hence satisfies the boundary conditions (5.11), (5.13). This proves the uniqueness.
To show that yZ is transversal let us assume that it is nontransversal. By Definition 3.2,
Theorem 3.1 and Lemma 2.4 there is a sequence xZ 2 SZ such that
Applying Lemma 5.1 to instead of xn we get
In a similar way
By the positive invariance of W s
loc and the negative invariance of W u
loc we have
Furthermore
~
and as a byproduct of Proposition 3.3 (see (3.13))
From (5.12), (5.15), (5.16) and (5.17) we obtain
Finally, we impose on " 0 the restriction
which yields kB 0
16C and similarly kB 0
16C . Using these estimates
and (5.19) in (5.18) we obtain
In particular, kPuxn
2 and the inequality (5.5) from Lemma 5.1 and (5.20) yield
In an analogous way we find kxn thus arrive at a contradiction to



--R

An Introduction
in Continuation and Bifurcation: Numerical Techniques and Applications
The effect of discretization on homoclinic orbits
The numerical computation of connecting orbits in dynamical systems

in Chaotic Numerics
Shadowing orbits of ordinary differential equations
An Introduction to Chaotic Dynamical Systems
Numerical computation of heteroclinic orbits
"invisible"
Numerical computation of invariant manifolds connecting fixed points
Matrix Computations
Dynamical Systems
Dynamics and Bifurcation
Geometric Theory of Semilinear Parabolic Equations
Smooth Dynamical Systems
Berechnung diskreter homokliner Orbits
Rigorous chaos verification in discrete dynamical systems
Hyperbolicity and Sensitive Chaotic Dynamics at Homoclinic Bifurcation

the shadowing lemma and transversal homoclinic points

Differentiable dynamical systems

--TR

--CTR
Y.-K. Zou , W.-J. Beyn, On manifolds of connecting orbits in discretizations of dynamical systems, Nonlinear Analysis: Theory, Methods & Applications, v.52 n.5, p.1499-1520, February
numerical methods;shadowing;homoclinic points for maps;dynamical systems