0 t!?1 t!1 t! - CiteSeerX
COMPUTATION AND CONTINUATION OF HOMOCLINIC AND HETEROCLINIC ORBITS WITH ARCLENGTH PARAMETERIZATION LIXIN LIUy , GERALD MOOREz , AND ROBERT D. RUSSELLx
Abstract. In this paper, we study a numerical method for the computation and continuation of homoclinic and heteroclinic orbits based upon the arclength parameterization of the orbits. Unlike most other methods, this method utilizes the geometric structure of the homoclinic and heteroclinic orbits and does not require solving a boundary value problem on an in nite interval. However, the boundary value problem formulated by this method can have a singularity at the end of the domain, and thus we introduce a special collocation method to handle such a singularity. We discuss the convergence properties of our collocation method and the implementation of the method which uses the software AUTO. For several examples we show that the arclength parameterization compares very favorably with the other numerical methods, although there are some limitations in the Sil'nikov case. Key words. homoclinicorbits, heteroclinic orbits, collocationmethods, numerical continuation. AMS subject classi cations. 34C37, 65L10, 65L60.
1. Introduction. Numerical computation of homoclinic and heteroclinic orbits is of interest in a variety of contexts. For example, structural changes in dynamical systems are often related to the appearance or disappearance of trajectories connecting one or several stationary points. Homoclinic orbits often arise as the limiting case of periodic solutions while heteroclinic orbits can represent travelling wave solutions of parabolic partial di erential equations. In the global bifurcation analysis, the occurrence of homoclinic orbits is closely tied to chaotic behavior in dynamical systems [18]. Consider the parameterized ODE system (1.1) u0(t) = f (u(t); ); ?1 < t < 1 where u 2 Rn and 2 Rp . For a given , if there exists a non-constant solution u(t) of (1.1) such that (1.2) u? = t!?1 lim u(t); u+ = tlim !1 u(t); f (u ; ) = 0; then the pair (u(t); ) is called a connecting orbit between stationary solutions or simply a connecting orbit. Moreover, if u? = u+ , the orbit is called a homoclinic orbit; otherwise, it is called a heteroclinic orbit. We assume both u+ and u? are hyperbolic xed points. Let MS and MU be the stable and unstable manifolds of these two stationary solutions de ned by t lim t(u) = u g; MS = fu 2 Rn : tlim MU = fu 2 Rn : t!?1 !1 (u) = u g; Received by the editors July 1, 1995; accepted for publication (in revised form) December 20, 1995. y Department of Computer Science, Concordia University, Montreal, Qecbec, Canada H3G 1M8. Current address: Nortel Technology, P.O. Box 3511, Station C, Ottawa, ON, Canada K1Y 4H7 ([email protected]). z Department of Mathematics, Imperial College, Queen's Gate, London, UK SW7 2BZ ([email protected]). x Department of Mathematics & Statistics, Simon Fraser University, Burnaby, BC, Canada V5A 1S6 ([email protected]). 1
2
L. LIU, G. MOORE, AND R. D. RUSSELL
with dimensions dim(MU? ) = n? ; dim(MS? ) = n ? n? ; dim(MS+ ) = n+ ; dim(MU+ ) = n ? n+ ; where t is the ow corresponding to (1.1). Note that when MU? and MS+ intersect, MU? \ MS+ is at least one dimensional. Moreover, if n? + n+ = n + 1, then generally there exists exactly one orbit between u? and u+ . When n? + n+ > n + 1, dim(MU? \ MS+ ) > 1, and further conditions are required to obtain the unique connecting orbit (e.g. see [18]). When n? + n+ < n + 1, the connecting orbit is not structurally stable, and we need to add some free parameters to establish the stable connection between u? and u+ . In this case, the number of free parameters is (1.3)
p = n + 1 ? (n? + n+ ):
In particular, for homoclinic orbits p = 1, or n? +n+ = n. A connecting orbit (u(t); ) between stationary solutions u? and u+ is non-degenerate if both are hyperbolic xed points, A = t!1 lim f u (u(t); ) = f u (u ; ); p = n + 1 ? (n? + n+ ); and the ODE system
v0 = f u(u; )v + f (u; ) has only the trivial solution = 0, v = cu0 for some c 2 R [7]. Throughout the
(1.4)
paper, we assume the connecting orbits are non-degenerate. In section 2, we review the methods developed in [6], [7], [14], [15], [16] and [17] for the computation of connecting orbits. These are all based upon solving boundary value problems (BVPs) in a truncated nite interval [T?; T+ ] (and usually this interval is rescaled to [0; 1]), so one can apply standard BVP solvers such as collocation methods or multiple shooting methods to compute the solutions. One drawback is that the derivative of the solution is usually quite large in a narrow region, and the BVP is similar to a singular perturbation problem. In section 3, we discuss an arclength parameterization method for the computation of homoclinic and heteroclinic orbits and their solution branches using continuation. The main advantage of this method is that the BVP is de ned on a nite interval using this natural \phase space scaling" (the total arclength). It thus avoids truncating and rescaling the in nite interval. Moreover, the solution is smooth over the whole domain, except possibly at the end points. When f u (u; ) at u? =u+ has dominant real positive/negative eigenvalues, this method has a high rate of convergence and has distinct computational advantages over the previous methods (allowing, e.g., the possibility to trace saddle node homoclinic orbits of very large length). In the case of Sil'nikov-type problems, the simple use of arclength implemented here has some distinct limitations, although we discuss more sophisticated alternative implementations which warrant pursuing in the future. In [25] the arclength parametrization method has been implemented in the well-known software package AUTO [14]. In section 4 we give a brief discussion of this implementation, which requires several signi cant changes to AUTO itself. Numerical examples which demonstrate the robustness of the approach are presented in section 5.
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
3
2. Numerical Methods for Computing Connecting Orbits. The strategy used in numerical methods for computing connecting orbits can usually be summarized as solving the boundary value problem (2.1) u0 ? f (u; ) = 0; (2.2) B(u(T? ); u(T+ ); ) = 0; (2.3) J (u; ) = 0; on an interval J = [T? ; T+ ], T? < 0, T+ > 0. Here, B(u(T? ); u(T+ ); ) are the boundary conditions and J (u; ) is a so-called phase condition. The methods are generally distinguished from one another by the way in which di erent types of boundary conditions and phase conditions are introduced. For convenience, we consider the BVP in the normalized interval t 2 [0; 1], so (2.4) u0 ? T f (u; ) = 0; (2.5) B(u(0); u(1); ) = 0; (2.6) (u; ) = 0; where T = T+ ? T? . In general, a robust choice of the phase condition is the integral phase condition (2.7)
(u; ) :=
or (2.8)
(u; ) :=
1
Z
0
1
Z
0
[u(t) ? u0 (t)]T u0 (t)dt = 0
[f (u; ) ? f (u0 ; 0)]T f u (u; )f (u; )dt = 0
where u0(t) is a nearby known solution (usually the solution obtained from the previous continuation step). These conditions can be obtained by minimizingthe functional (2.9)
E() :=
or (2.10)
G() :=
1
Z
0 Z
0
1
ku~ (t + ) ? u0 (t)k2dt ku~ 0 (t + ) ? u00(t)k2 dt
against the phase shift . The phase condition (2.7) is used in [6], [7] and [14] while (2.8) is used in [15]{[17]. The boundary conditions (2.5) can be dealt with in a variety of ways. 2.1. Periodic Boundary Conditions. The periodic boundary conditions are given by (2.11) u(0) = u(1): In many applications, homoclinic orbits occur as limiting cases of periodic solutions. It is natural to use these periodic solutions to approach these connecting orbits, and it can in fact be shown that there exists a branch of periodic solutions near a homoclinic orbit [6]. The computation of a connecting orbit this way is usually achieved
4
L. LIU, G. MOORE, AND R. D. RUSSELL
by the continuation of a branch of periodic solutions where the period T tends to in nity when the connecting orbit is approached. This method is basically available for computing homoclinic orbits in the software AUTO. There is one free parameter in the system, which means that p = 1 and n? +n+ = n. This condition always holds for homoclinic orbits. One disadvantage of the periodic boundary condition approach is that it is not generally suitable for the computation of heteroclinic orbits. 2.2. Eigenvector Boundary Conditions. The eigenvector boundary condition is introduced in [15]{[17]. If u are stationary solutions of (1.1),
f (u?; ) = 0;
f (u+ ; ) = 0: ? , and (u ; ) has If f u(u? ; ) has n? eigenvalues with positive real parts, f?igni=1 fu + n + n+ eigenvalues with negative real parts, f+i gi=1 , satisfying (2.13) f u(u; )vi = ivi ; kvik = 1; i = 1; : : :; n;
(2.12)
then the eigenvector boundary conditions can be de ned as (2.14) u(0) ? u? = ?
n
? X
i=1
c?i v?i ; u(1) ? u+ = +
n
+ X
i=1
c+i v+i ;
n
X
i=1
c2i = 1
are constants. where and fci gni=1 The geometric interpretation of the eigenvector boundary conditions is to require that the endpoint u(0)=u(1) of the orbit u(t) lie on the unstable/stable tangent space of u? =u+ and be a certain distance from the xed point u? =u+ . These tangent ? , and fv gn+ , subspaces, denoted as E?U () and E+S (), are spanned by fv?i gni=1 +i i=1 respectively. When the orbit is non-degenerate, there are (n + 2)(nn? + n+ ) + p + 2 = , and ), n (n + 1)(n? + n+ ) + n + 3 unknown parameters (fi; vi; ci gi=1 di erential equations and (n+1)(n? +n+ )+2n+3 constraints. A more straightforward as orthonormal bases for the tangent implementation would be to choose fvigni=1 S U subspaces E? () and E+ (). 2.3. Projection Boundary Conditions. An alternative (but mathematically equivalent) approach uses the projection boundary conditions (also known as asymptotic boundary conditions) [6], [7]. The boundary conditions again take the form
(2.15)
u(0) ? u? 2 E?U (); u(1) ? u+ 2 E+S ();
but for numerical calculations the projection boundary conditions are often presented in the form (2.16) B? (u(0); ) := P? ()[u(0) ? u? ] = 0; B+ (u(1); ) := P+ ()[u(1) ? u+ ] = 0 where the rows of the matrices P? () 2 R(n?n?)n and P+ () 2 R(n?n+)n span the stable subspace of f Tu (u? ; ) and the unstable subspace of f Tu (u+ ; ), respectively. Note that (2.15) and (2.16) are equivalent when u are both hyperbolic xed points. Equation (2.16) de nes 2n ? (n? + n+ ) boundary conditions. The BVP (2.4){(2.6) is well-posed if the connecting orbit is non-degenerate, as the number of free parameters is p = n + 1 ? (n? + n+ ). In general, one cannot handle non-hyperbolic xed points without some modi cation. However, Schecter [30], [31] shows that under certain conditions the projection
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
5
boundary conditions can be modi ed to handle certain types of connecting orbits with semi-hyperbolic xed points (quadratic turning points). The equivalent treatment for the eigenvector boundary condition approach is in [17]. The computation of fully non-hyperbolic homoclinic orbits is studied in [4], [10]. Based on the projection boundary conditions, Champneys and Kuznetsov [11] develop a method for the detection and continuation of several di erent types of codimension-two homoclinic bifurcation points. For hyperbolic xed points, these types of boundary conditions all have the property that convergence to the connecting orbit is exponential with respect to the length of the truncated interval. The eigenvector and projection boundary conditions have a better rate of convergence than the periodic boundary conditions, and they can be used to handle general types of heteroclinic orbits (p 6= 1). Both approaches are similar to the one used in [24] for solving BVPs over in nite intervals, where the boundary conditions are obtained from linear approximations of the stable and unstable subspaces of the xed points. However, because of their di erent implementations, each has its own advantages and disadvantages. For the eigenvector boundary condition approach, the dimension of the original problem is increased because of the introduction of the unknown parameters fi; vi ; cig and . Each requires certain knowledge of the invariant subspace structure at the xed points, one a smooth basis and the other a smooth projection. A smooth QR factorization is desirable in such a case (e.g. see [7]). Finally, to perform numerical continuation for the connecting orbit branches, we add a pseudo-arclength continuation equation [22] (2.17)
2
u
Z
0
1
[u(t) ? u0 (t)]T u_ 0 (t)dt + 2 ( ? 0)T _ 0 + 2 ( ? 0 )_ 0 = s
to (2.4){(2.6), where u , and are scaling parameters, (u0 ; 0; 0 ) is the solution at a previous continuation step, (u_ 0 ; _ 0; _ 0 ) denotes the derivative along the solution branch, and s is the continuation stepsize. The number of free parameters in such a situation is n + 2 ? (n1 + n+ ) (an extra parameter is introduced for continuation purposes). The continuation equation (2.17) is used for solving BVPs in AUTO. 3. Arclength Parameterization. All three methods discussed in the previous section have the disadvantage that when T becomes very large, u0 (t) varies fast (O(T)) in only a small part of the domain [0; 1], and on the rest of the domain u0 (t) 0. In fact, the BVP system (2.4){(2.6) is similar to a singular perturbation problem, where the existence of an interior layer or boundary layer can cause great numerical diculty. However, in phase space the connecting solution itself often remains bounded and smooth except at the xed points. This suggests that a geometrical formulation could lead to a superior method. Moore [28] has recently introduced such a geometrical approximation based upon an arclength parameterization. 3.1. Parameterization and Continuation. Let L be the total arclength of the solution u(t) in phase space, i.e. (3.1)
L=
Z
1
?1
ku0(t)kdt;
and let be the scaled arclength parameterization variable Z t 1 (3.2) ku0( )kd : (t) = L ?1
6
L. LIU, G. MOORE, AND R. D. RUSSELL
Then, if t() denotes the inverse of (t), (3.3) v() := u(t()) satis es the arclength parameterized di erential equation (); ) ; 0 < < 1: (3.4) v0() = L kff ((vv(); )k Note that the endpoints are now singular points for (3.4), and away from them kv0()k = L. Assuming that the connecting orbit is not of Sil'nikov type, i.e., the dominant eigenvalue of f u (u ; ) is real (the dominant eigenvalue for f u (u? ; ) is that with smallest positive real part, while the dominant eigenvalue for f u (u+ ; ) is that with largest negative real part), 0 lim v0 () and lim !0 !1 v ()
exist and lie in the dominant eigenspaces of u? and u+ , respectively. The continuity conditions require (3.5) kv0 (0)k = kv0 (1)k = L; and the boundary conditions (3.6) v(0) = t!?1 lim u(t) = u? ; v(1) = tlim !1 u(t) = u+ ensure that the xed points u? and u+ are the endpoint values for the solution. In this formulation a phase condition is not needed, as the endpoints are xed and L is bounded. The projection boundary conditions (2.16) become (3.7) P?()v0 (0) = 0; P+ ()v0 (1) = 0 where the projection matrices P? () 2 R(n?n? )n and P+ () 2 R(n?n+)n are de ned as in section 2. Moore proves that the BVP system (3.4){(3.7) is well-posed when the connecting orbit is non-degenerate [28]. A pseudo-arclength continuation equation (3.8) v2
Z
0
1
[v() ? v0 ()]T v_ 0()d + 2 ( ? 0)T _ 0 + L2 (L ? L0 )L_ 0 = s
may be used to perform numerical continuation for a branch of connecting orbits. The total number of free parameters (; L) for the system (3.4){(3.8) is n+2 ? (n? +n+ ). 3.2. Numerical Discretization. The BVP (3.4){(3.8) consists of n ODEs, one integral equation, 4n+2 ? (n? +n+ ) algebraic constraints, and p+1 parameters. Thus, this system has 2n extra conditions. Also, the BVP has singularities at = 0 and 1, where the right hand side of (3.4) is unde ned. Consequently, a standard boundary value solver is not suitable for use without some modi cation. Moore [28] recommends solving the BVP (3.4){(3.7) by a natural Gauss-Lobatto collocation method: Use a C 1 piecewise polynomial of degree m and a mesh 0 = 0 < < N = 1; Collocate (3.4) at all Lobatto points apart from = 0; 1; Require (3.5): kv0(0)k = kv0 (1)k = L; Require (3.6): v(0) = u? (); v(1) = u+ ();
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
Require (3.7): P? ()v0(0) = 0; P+ ()v0(1) = 0.
7
One of our goals has been to modify this algorithm for the software AUTO, which uses Lagrange polynomials as the basis functions and Gauss collocation points (and consequently involves di erent continuity conditions). Its modi cation to use Lobatto points would be dicult, so our implementation of the arclength formulation in AUTO uses Gauss collocation points in the interior subintervals, and for consistency between the equations and unknowns, Radau collocation points on the rst and last subintervals (see [25], [28]). The resulting numerical algorithm is the Gauss-Radau collocation method: Use a C 0 piecewise polynomial of degree m and a mesh 0 = 0 < < N = 1; Collocate (3.4) at m Gauss points fik gmk=1 in [i; i+1] for i = 1; :::; N ? 2; Collocate (3.4) at the left m ? 1 Radau points f0kgmk=2 in the rst subinterval [0; 1]; ?1 in the last subin Collocate (3.4) at the right m ? 1 Radau points fN ?1;k gmk=1 terval [N ?1 ; N ]; Require (3.5): kv0(0)k = kv0 (1)k = L; Require (3.6) or equivalently: f (v(0); ) = 0; f (v(1); ) = 0; Require (3.7): P? ()v0(0) = 0; P+ ()v0(1) = 0. For numerical continuation, the algorithm is supplemented with the continuation equation (3.8). A more detailed description of the implementation for AUTO is in the next section. 3.3. Convergence Results. The system (3.4){(3.7) is a nonlinear boundary value problem with a singularity of the rst kind at the endpoints. This type of singular BVP has been studied by de Hoog and Weiss [19], [20], and Weinmuller [33]{[36]. We now discuss the convergence properties of the Gauss-Radau collocation algorithm for solving (3.4){(3.7). The variational equation for (3.4) at the solution (v(); L; ) is dy ? L[I ? P()] f u (v; ) y = 0; 0 < < 1 (3.9) d kf (v; )k where 0()v 0()T f f T (3.10) P() = vv0() T v0 () = f T f is an orthogonal projection onto v0 (). If we assume that the dominant eigenvalues of the matrices f u (u ; ) are real and
v0 (0) := lim v0 () and v0 (1) := lim v0() !0+ !1?
exist and lie in their dominant eigenspaces, then we can rewrite (3.9) as y d M M 0 1 y (3.11) d ? + 1 ? ? B() = 0 where L[I ? P()]f u = 1 [I ? P(0)] (u ; ); (3.12a) M0 = lim fu ? !0 kf k ~1
8
L. LIU, G. MOORE, AND R. D. RUSSELL
(3.12b)
~1 = v0 (0)T f u (u? ; )v0(0)=L2 ;
(3.13a)
(1 ? )L[I ? P()]f u = 1 [I ? P(1)] (u ; ); M1 = lim fu + !1 kk ~
(3.13b) and
f
1
~1 = v0(1)T f u (u+ ; )v0 (1)=L2; B() = ?L[I ? P()] kf u((vv;;))k + M0 + 1M?1 :
f
We approximate y() by p(), a C 0 piecewise polynomial of degree m, satisfying the collocation equations 0 + M1 ? B( ) p( ) = 0; p0 (ik ) ? M (3.14) ik ik ik 1 ? ik where 0 i N ? 1; 1 k m, and ik = i + k hi are the collocation points in the interval [i ; i+1]. Note that for the Gauss-Radau collocation algorithm, (3.14) is not well-de ned at 0;1 = 0 and N ?1;m = 1, and these equations are replaced by the projection boundary conditions (3.15) P? ()p(0) = 0; P+ ()p(1) = 0: To simplify our analysis for the BVP (3.11), (3.15), let 2 0 21 = 2(1 ? ) 21 1 and ( ) 0 21 ; z( ) = y0( ) ; B() = B0 ( ) 0 21 : y() = yy0( ) 1 1 y1( ) B1 ( ) 12 1 1 2 Then (3.11) becomes dz 1 (3.16) d ? M z + D( )z = 0 where # " M1 ? B0 ( ) M 0 2 ? 2 M= ?M1 ; D( ) = ? 2M?0 + B12( ) : From (3.15), the boundary conditions become P? () (3.17) z (0) = 0 ; I ? I z(1) = 0: P+ () Our next result follows immediately from [19] Theorem 3.5. Lemma 3.1. Let D( ) 2 C p[0; 1], p 0. Then the BVP (3.16) and (3.17) has a unique solution z( ) such that
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
9
(i) z( ) 2 C p+1 (0; 1]. (ii) z( ) 2 C p+1 [0; 1] when all eigenvalues of M0 have non-positive real parts and all eigenvalues of M1 have non-negative real parts. (iii) Let be the eigenvalue of M0 with smallest positive real part, be the eigenvalue of M1 with largest negative real part and
= minfRe( ); ? Re( )g: Then
kz(p+1) ( )k const ?p?1 (j log jd ?1 + 1) kz(p+1) ( )k const (j log jd + 1) z( ) 2 C p+1 [0; 1]
p< < p+1
= p+1
> p+1 where d is the size of the largest Jordan block of M associated with or : The eigenvalues of M0 , M1 in (3.11) and the constants , in Lemma 3.1 have a close relationship with the eigenvalues of f u (u ; ). In fact, ordering the eigenvalues of f u(u? ; ) with positive real parts figki=1 such that 1 Re(2 ) Re(3 ) Re(k ), v0(0) is the approximate dominant eigenvector of f u (u? ; ), so from (3.12), ~1 1 . Moreover, v0(0) is an eigenvector of M0 with associated eigenvalue 1(M0 ) = 0. Other eigenvalues of M0 with positive real parts are i (M0 ) = i =~1 i =1 for 1 < i k, which implies 2=1 . Similarly, 1(M1 ) = 0, i (M1 ) ?i =1 for 1 < i l, and ?2 =1, where the eigenvalues of f u (u+ ; ) with negative real parts are ordered such that 1 Re(2) Re(3 ) Re(l ). Therefore,
min Re 2 ; Re 2 ; 1 1 and from Lemma 3.1, the smoothness of the solution z() depends on the size of . The analogue of (3.14) for solving (3.11) is the collocation method (3.18) p0 ( ik ) ? 1 M p( ik ) + D( ik )p( ik ) = 0; 0 i N ? 1; 1 k m ik for solving (3.16) where ik = i + k hi are the collocation points in [ i; i+1 ]. Singularity only occurs at = 0, and the collocation equations at this point are replaced by P?() p(1) = 0: p (3.19) (0) = 0 ; I ? I P+ () Theorem 3.2. Let D( ) 2 C m [0; 1]. For h = max0iN ?1fhi g suciently small, the solution p( ) of (3.18) and (3.19) satis es 8 const h j loghjd ?1; if 1 < < m; > > < m j loghjd ?1 ; const h if m < m + 1; kz( ) ? p( )k > const hm j loghjd ; if = m + 1; > : const hm ; if > m + 1 or there is no positive : Proof. Here we give the proof for the uniform mesh case; the extension to the non-uniform case can be obtained following the approach in [35]. (i) Let 1 < < m + 1. From Lemma 3.1, kz(m+1) ( )k const ?m?1 (1 + j log jd ?1 ):
10
L. LIU, G. MOORE, AND R. D. RUSSELL
De ne a degree m piecewise polynomialerror function e( ) satisfying the interpolation conditions e0( ik ) = z0 ( ik ) ? p0( ik ): Then, for 2 [ i ; i+1], m X ? ik 0 z0( ik ) ? p0 ( ) e ( ) = Lk h k=1 m X = z0 ( ) ? p0 ( ) + Lk ?h ik z0( ik ) ? z0( ) k=1 0 0 p = z ( ) ? ( ) + O(hm ?m?1 (1 + (log )d ?1)) where the Lk are the Lagrange basis functions. After integration, e( ) = z( ) ? p( ) + O(hm ?m (1 + (log )d ?1 )): It follows that e( ) satis es (3.20) e0( ik ) ? 1 M e( ik ) + D( ik )e( ik ) = O(hm ik ?m?1 (log ik )d ?1) ik for 0 i N ? 1; 1 k m, and the boundary conditions P ( ) ? (3.21) I ?I e(1) = 0: P+ () e(0) = 0; In a similar manner as that in the proof of Lemma 3.8 in [35], we have kz( ) ? p( )k const hm j loghjd ?1 ; = m; kz( ) ? p( )k const hm j loghjd ?1 j1 ? ?m j const hm j loghjd ?1 > m; const h j loghjd ?1
< m: (ii) If = m + 1, then kz(m+1) ( )k const (1 + j log jd ); and the right hand side of (3.20) is replaced by O(hm (1 + (log )d ), which gives kz( ) ? p( )k const hm j loghjd : (iii) If > m + 1, then z 2 C m+1 [0; 1] and the right hand side of (3.20) is replaced by O(hm ). Therefore, kz( ) ? p( )k const hm : Finally, the proof for the non-positive case is given in Theorem 4.5 of [20]. One advantage of the collocation method when applied to standard boundary value problems (without singular points) is the property of superconvergence for suitably chosen 1; ; m . If the solution is suciently smooth and (3.22)
Z
0
1
sk w(s)ds = 0;
k = 0; ; < m;
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
then
ky(i ) ? p(i )k const hm+ +1 ;
11
i = 0; ; N: In particular, for the Gauss points, the rate of convergence is O(h2m ); for the Lobatto points, it is O(h2m?2); and for the Radau points, the superconvergence rate is O(h2m?1). In fact, applying the similar technique from [20] Theorem 5.1 and [35] Theorem 4.1, the rate of convergence in Theorem 3.2 can be improved to O(hm+1 ) provided D( ) 2 C m+2 [0; 1], > m + 2 and (3.22) holds. This rate is most likely optimal (see examples given in [20], [35]). It should be noticed that the above convergence results are mainly dependent upon the continuity conditions of the solution z( ), which is in turn determined by the ratio of the two smallest eigenvalues (in magnitude) of f u (u ; ). In general, when z( ) does not satisfy the continuity conditions, the optimal rate of convergence cannot be obtained. In [9], de Boor shows that a non-smooth function such as f(t) = t with 2 (0; 1) can be approximated by a piecewise polynomial of order k with accuracy O(hk ) if jf (k) (t)j decreases monotonically in (0; 1] and the mesh points are given by the equidistribution method Z 1 Z ti+1 1 ( k ) 1 =k jf (t)j dt = N jf (k) (t)j1=kdt; i = 0; 1; ; N ? 1: (3.23) ti 0 In AUTO, the equidistribution method (3.23) is used for k = m + 1 where m is the number of collocation points per subinterval. In general, each jzi(m+1) ()j may not be monotonically decreasing on (0; 1). However, for the non-Sil'nikov case, all jzi(m+1) ()j decrease monotonically in (0; ] and increase monotonically in [1 ? ; 1) for suciently small . Therefore, with the proper mesh selection scheme (base on (3.23)), the collocation method can be expected to have a rate of convergence at least O(hm+1 ) regardless of the continuity of the solution z() at the end-point. Our numerical examples presented in section 5 indicate that this is the case. 3.4. S il'nikov-type Connecting Orbits. For the Sil'nikov case v0(0) and/or 0 v (1) do not exist. As a result, the arclength parametrization is more complicated in principle, and practical results using a straightforward implementation have been less good. However, we believe that even then, using arclength away from the xed points is an advantageous procedure. Therefore, in this section, we show how the arclength idea may be adapted to this more dicult situation. For de niteness we assume that the Sil'nikov behaviour occurs at u? but not at u+ , although the algorithms below generalize easily. In order to avoid the singularity at u? (), we introduce a boundary condition involving the unstable manifold of u? () rather than u? () itself. Thus we have the di erential equation (3.4), i.e. (); ) ; 0 < < 1; (3.24) v0() = L kff ((vv(); )k with boundary conditions (3.25) v(0) 2 MU? at a \distance" ? from u? (); v(1) = u+ () and end conditions (3.26) P+ v0(1) = 0; kv0(1)k = L:
12
L. LIU, G. MOORE, AND R. D. RUSSELL
Here ? > 0 is a small user-chosen constant, which is analogous to the choice of T? in section 2. This set of equations may be discretized by Gauss collocation on all subintervals apart from the last, where Gauss-Radau collocation should be used with (3.26) replacing the di erential equation at = 1. The remaining question is how to make the boundary condition at = 0 practical and, in particular, specify what we mean by \distance". There are two possibilities. 1) Using the unstable subspace E?U (). Here we simply mimic (2.14) and take v(0) to be in the unstable subspace of u? (), i.e. (3.27)
v(0) ? u?() = ?
n? X i=1
c?iv?i ;
n? X i=1
c2?i = 1:
The di erence is that now ? is known and only c?1 ; : : :; c?n? are to be determined. The size of ? controls the O(2?) error introduced by approximatingthe manifold with ? should be chosen the subspace. Moreover, to avoid possible ill-conditioning, fv?i gni=1 as a well-conditioned basis for E?U () for which the restriction fu(u? (); ) to E?U () is positive-de nite, as simply choosing an orthonormal basis is usually unacceptable [29]. 2) Approximating the manifold more accurately. In [29] a number of algorithms are described for computing the unstable manifold of a given xed point. For our connecting orbit problem, we think the potentially most useful approach is to solve the system (3.28) w0( ) = 1 f (w( ); ); 0 < < 1 ? with boundary condition
w(0) = u?();
(3.29) end condition
P? ()w0 (0) = 0;
(3.30) and integral side-constraint (3.31)
Z
1 kf (w( ); )k
d = ? :
? Now ? measures distance in terms of arclength, and ? > 0 is a user-chosen constant which corresponds to the explicit exponential change-of-variable 0
= e? t ; w( ) u(t) (see [29]). Hence, choosing ? < min1in? fRe(?i )g (notation as in (2.13)) forces w0 (0) to be zero. Equation (3.28) may be discretized by Gauss collocation on all subintervals except the rst, where Gauss-Radau is used with (3.30) replacing the di erential equation. This is also the natural quadrature rule for (3.31), making use of the fact that the integrand is zero at = 0. Note that (3.28) should not require many subintervals for accurate approximation. Finally (3.28){(3.31) are coupled to (3.24){(3.26) through the continuity condition (3.32) v(0) = w(1):
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
13
To conclude, we have two techniques for dealing with Sil'nikov behaviour or, more generally, with v() having a low degree of smoothness at an endpoint. The rst is a direct analogue of those described in sections 2.2 and 2.3 (which could, in fact, be used at all endpoints). The second, although requiring further work to establish the optimal choices of ?; ? and the mesh for (3.28), warrants further investigation as a more sophisticated approach to solve this dicult class of problems. 4. Software Implementation. We now describe our implementation of the Gauss-Radau collocation algorithm of section 3 in the software AUTO. For this algorithm there are only (m ? 1)n collocation equations in the rst and the last subintervals, so n equations are needed to replace the collocation equations at = 0 in the rst subinterval and at = 1 in the last subinterval. To do this, we use the equation (4.1) f (v(0); ) = 0 to replace the collocation equation at = 0, and (4.2) f (v(1); ) = 0 at = 1. Equations (3.5), (3.7) and (3.8) are treated as boundary conditions, and these are placed in the bottom rows of the collocation matrix. Here, v0 (0) and v0(1) are approximated by mth order Lagrange interpolation for non-Sil'nikov case. The projection matrices P?() and P+ () are obtained by using the LAPACK [1] subroutine DGEES. The projection matrices are not unique, and care must be taken to maintain their smoothness. The nal structure of the linear system is shown in Fig. 4.1(a) for the case of 2 di erential equations, 4 mesh intervals, 3 collocation points and 3 free parameters. When solving a linear system arising from discretization of a BVP, the rst step in AUTO is to use a local condensation method to eliminate the unknowns at non-mesh points. This process uses no row pivoting, so in our case reordering of the equations in the rst block is necessary to avoid breakdown. We place the end conditions (4.1) at the bottom of the rst block, and the coecient matrix pro le is then as shown in Fig. 4.1(b) [25]. The condensed linear system has the (almost) block bi-diagonal structure 3 2 A1 C1 D1 6 A2 C2 D2 77 6 6 .. 77 . . .. .. (4.3) 6 . 7 6 4 AN CN DN 5 B0 B1 BN ?1 BN E where the blocks C1 and AN are zero matrices. Numerical methods for solving linear systems like (4.3) are studied in [26]; in general, either a QR factorization or LU decomposition with row pivoting can be applied to eliminate the interior mesh points. After such an elimination process, the coecient matrix is reduced to the form 2 0 R1 F1 H1 3 6 0 R2 F2 H2 77 6 6 .. 77 . . . .. .. .. 6 . 7: (4.4) 6 6 GN ?1 0 F H N ? 1 N ?1 7 6 7 4 0 C~N D~ N 5 B0 B1 BN ?1 BN E
14
L. LIU, G. MOORE, AND R. D. RUSSELL r r r r r r
r r r r r
r r r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r r r r r r r
r r r r r
r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r r r
r r r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r
r r r r r r r r r r
r r r r
r r r r
r r r r
r r r r
r r r r r r r r r r r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r r r r r r r
r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
(a) r r r r r r
r r r r r
r r r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r
r r r r r r
r r r r r r
r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r r r
r r r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r
r r r r r r r r r r
r r r r
r r r r
r r r r
r r r r
r r r r r r r r r r r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r r r r r r r
(b) Fig. 4.1. Linear system structure: (a) before reordering the equations; (b) after reordering the
equations
The blocks B1 ; ; BN ?2 can be eliminated using R1; ; RN ?2. However, if not further modi ed, the linear system solver of AUTO will fail when eliminating the block BN ?1 . It is sucient to reorder the equations and unknowns in such a way that the linear system has the form 3 2 A2 C2 D2 6 A3 C3 D3 77 6 6 .. 77 . . .. .. 6 . 7: (4.5) 6 6 0 C D N N 7 7 6 4 0 A1 D1 5 B1 B2 BN B0 E Because v(0) and v(1) are hyperbolic xed points, the matrices A1 and CN are
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
15
non-singular, and the linear system can be solved by the QR factorization or LU decomposition as before. The nal reduced linear system is 2 3 G2 R2 F2 H2 6 G3 R3 F3 H3 77 6 6 . .. 77 . . .. .. 6 .. . 7: (4.6) 6 6 0 R 0 H N N 7 6 7 4 0 A N D1 5 B1 B2 BN B0 E In order to be able to use AUTO with the arclength implementation, we have also modi ed various parts of the code to allow users to provide discrete solution data for the starting solution, and we have implemented the arclength formulation for continuation of periodic solution branches. More speci cally, the users can switch from the \time" formulation to the \arclength" formulation directly. However, several quantities, such as the solution norm, then have di erent interpretations, and the output data must be interpreted judiciously. Also, we should mention that our current software implementation is not suitable for solving large dimensional problems (n > 20). Some major changes in the basic data structure and subroutine interdependency structure in AUTO are necessary before this implementation can be realized. 5. Numerical Examples. In this section, we consider several numerical examples and show the eciency of the arclength parameterization method by comparing it to the methods given in section 2. We refer to the arclength parameterization method as simply the arclength method, and we call the (time parametrized) implementations of periodic boundary conditions, eigenvector boundary conditions and the projection boundary conditions as the periodic method, eigenvector method and projection method, respectively. In all of our computations, we use the software AUTO, having made the necessary modi cations of it for each method. Therefore, the BVP is solved by a spline collocation method with Lagrange basis functions. In all cases, we choose 4 collocation points on each mesh subinterval (the default for AUTO). The numerical tests are performed on a Silicon Graphics Indigo running IRIX 4.0.5F with a MIPS FORTRAN 77 compiler with the \?g" option for debugging purposes. Example 5.1. The Nagumo Equations u01 = u2; u02 = cu2 ? u1 (1 ? u1 )(u1 ? a);
0 < a < 1:
This system has been used in [7], [15] and [27] as a test problem for computing heteroclinic orbits (a travelling wave solution branch) between the xed points (0; 0) and (1; 0). The eigenvector method and the arclength method are tested for this problem. For the eigenvector method, we use 10 mesh intervals (NTST=10). As in [15], we start with the exact solution p exp( 22pt) ; u2(t) = u01(t) u1(t) = 2 1 + exp( 2 t)
at (a; c) = ( 12 ; 0) and T = 5 and follow the solution by increasing T to 1000. Fixing T at this point and using a as the primary continuation parameter and c is thepsecondary continuation parameter, we compute both solution branches for c = 2(a ? 12 ).
16
L. LIU, G. MOORE, AND R. D. RUSSELL Table 5.1
Relative error for b at c = 14
NTST 12 25 50 100
Projection method Periodic method Arclength method { { 8.05e-05 1.42e-05 1.69e-05 4.59e-07 8.31e-09 1.36e-09 6.09e-09 < 1.0e-11 1.58e-09 8.08e-11
The eigenvector method fails after only 6 steps. After increasing the number of mesh subintervals to NTST=15, we successfully compute both branches. The absolute error for the parameter c in the lower half branch varies from O(10?8) to O(10?4). For the arclength method, we successfully compute both branches with only 10 subintervals (NTST=10), and the absolute error for c is O(10?8) for all solutions (see Fig. 5.1(a)). In this computation, the tolerance for the Newton iteration is set to 10?6. The distribution of the mesh points for the lower half branch solution is shown in Figs. 5.1(b) and (c) for the arclength and eigenvector method, respectively. Note that there are 41 points for the arclength method (10 mesh subintervals and 4 collocation points per interval) and 61 points for the eigenvector method (15 mesh subintervals and 4 collocation points per interval). In this and the other examples, the size of NTST is a good measure of the relative speed of the methods. Since both xed points for the heteroclinic orbits are saddle points and there is one positive and one negative eigenvalue in the Jacobian at (0; 0) and (1; 0), the solution v() 2 C 1 [0; 1], and the mesh points are smoothly distributed on [0; 1] for the arclength method. Example 5.2. The Lorenz equations x0 = a(y ? x); y0 = bx ? y ? xz; z 0 = xy ? cz:
At the Lorenz value (a; b; c) = (10; 13:92656; 38 ), there exists a homoclinic orbit with L = 64:89938. Starting with the Hopf bifurcation point at b = 24:73684, we use AUTO to trace out the periodic solution branch for (a; c) = (10; 83 ) until very close to the homoclinic orbit. Then we use several methods to continue the homoclinic orbit branch. The tolerance for the Newton iteration is set to 10?8. The primary and secondary continuation parameters are c and b. Their relation is shown in Fig. 5.2(a), and the solution in phase space is shown in Fig. 5.2(b). All homoclinic orbits on this branch have the stationary point (0; 0; 0). The eigenvalues of the Jacobian at the xed points are i h p ?c; ? 21 (a + 1) (a + 1)2 + 4a(b ? 1) : In Fig. 5.2(c) we plot the ratio between two negative eigenvalues of the Jacobian at (0; 0; 0). For the arclength parameterization method, using only 25 mesh intervals (NTST =25), we successfully continue the homoclinic orbit branch until c = 14:376. The nal parameter value for b is 95040:54 and the total arclength L = 592401:13. When we use NTST=50, we have no trouble in continuing the homoclinic branch to c = 14:455. Next, we compute this branch of homoclinic orbits with the projection method and the periodic method. We start both computations from a periodic solution with the xed period T = 1000 and use c and b as the primary and the secondary continuation parameters. Using 25 mesh subintervals and a tolerance of 10?8, we are unable to
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
17
1e-3 Arclength Method (NTST=10) Eigenvector Method (NTST=15)
Absolute error for c
1e-4 1e-5 1e-6 1e-7 1e-8 1e-9
0
0.2
0.4
0.6
0.8
1
0.6
0.8
1
0.6
0.8
1
a
(a) u1 1 0.8 0.6 0.4 0.2 0 0
0.2
0.4
σ
(b) u1 1 0.8 0.6 0.4 0.2 0 0
0.2
0.4
t
(c) Fig. 5.1. The Nagumo equations: (a) absolute error for c; (b) distribution of the mesh points for the arclength method; (c) distribution of the mesh points for the eigenvector method.
18
L. LIU, G. MOORE, AND R. D. RUSSELL b 1e+6
1e+5
1e+4
1e+3
1e+2
1e+1
2
4
6
8
10
12
14
16
c
(a) z 300 250 200 150 100 8/3 14/3 20/3 26/3 32/3
50 0
0
10
20
30
40
50
60
70
10
12
14
16
x
(b) λ2/λ1 20 18 16 14 12 10 8 6 4 2
2
4
6
8
c
(c) Fig. 5.2. Homoclinic orbits in the Lorenz equations: (a) solution branch in parameter c ? b space; (b) some homoclinic orbits in the x ? z plane on the branch; (c) ratio between two negative eigenvalues of the Jacobian at (0; 0; 0).
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
19
Table 5.2
Absolute error for b at c = 8:3
NTST 10 20 40 80 160
m=2 9.69e-02 2.65e-02 1.92e-03 1.21e-04 7.42e-06
m=3 5.90e-03 2.48e-04 1.53e-05 8.96e-07 5.58e-08
m=4 2.61e-05 9.45e-07 1.53e-08 2.30e-10
m=5 m=6 3.49e-06 1.88e-07 7.29e-08 5.60e-10 1.51e-09 3.00e-11
Table 5.3
Convergence rate for b at c = 8:3
NTST m = 2 10!20 1.87 20!40 3.78 40!80 3.99 80!160 4.03
m=3 m=4 m=5 m=6 4.57 4.79 5.58 8.39 4.02 5.94 5.59 4.09 6.06 5.65 4.00
perform continuation for either method. Taking the larger tolerance 10?6, we can continue the solution branch with both methods. The projection method fails at c = 14:328 and the periodic method fails at c = 13:920. However, the projection method cannot locate the solution at c = 14 during this continuation. In order to obtain the solution at c = 14, we must further raise the tolerance to 10?5. Both methods are then able to obtain the required solution. In our numerical results, the \exact" value of b and L is obtained by both arclength and projection methods using NTST=200 (they agree to the digits printed). In Table 5.1 we give the relative errors in the computed values of b for c = 14 using various values of NTST. Superconvergence is clearly observed for the arclength method and the projection method. For NTST=25 or less, the arclength method is far more accurate than the projection and periodic methods. For the periodic method, the accuracy does not improve when going from NTST=50 to 100, presumably because T = 1000 is not suciently large for this method. In Figs. 5.3{5.4, the distribution of the mesh points for c = 14 is shown for all three methods with NTST=50. Note, in Figs. 5.3(a) and 5.4(a), because of the lack of a higher order derivative for v(), there are more mesh points near = 1 to maintain the higher rate of convergence. To verify the rate of convergence for our collocation method, we compute the homoclinic orbit for the Lorenz equation at (a; c) = (10; 8:3). At this value, b = 60:5547266, the total arclength L = 303:0482482, and as we observe from Fig. 5.2, the ratio between the two negative eigenvalues of the Jacobian at the xed point (0; 0; 0) (1 = ?8:3, 2 = ?30:51594) is near minimal ( = 12 3:6766). Therefore, from our analysis in section 3.3, v() is at most C 3[0; 1]. In Tables 5.2{5.5, we show the absolute error and corresponding approximate rate of convergence for the parameters b and L. The rate of convergence is calculated by log2 (N =2N ) where N is the absolute error obtained with N mesh intervals. For the odd and even number of collocation points, it is roughly O(hm+1 ) and O(hm+2 ), respectively. In all cases, the rate of convergence is higher than the degree of continuity of v(). Clearly, the mesh selection scheme (3.23) plays the key role for the higher rate of convergence (better than O(h3)) in this example. Examples 5.1 and 5.2 clearly show that the arclength parameterization method can be very ecient for computing connecting orbits. However, recall that special care is needed if the dominant eigenvalues at the xed points are not real, i.e., the orbit
20
L. LIU, G. MOORE, AND R. D. RUSSELL x 500 450 400 350 300 250 200 150 100 50 0
0
0.2
0.4
0.6
0.8
1
σ
(a) x 500
400
300
200
100
0 0.05
0.055
0.06
0.065
0.07
0.075
0.08
0.085
0.09
0.095
0.1
t
(b) x 500
400
300
200
100
0 0.995
0.996
0.997
0.998
0.999
1
t
(c) Fig. 5.3. Distribution of the mesh points in Example 5.2 for x: (a) the arclength method; (b) the projection method; (c) the periodic method.
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
21
y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
300
350
400
450
500
300
350
400
450
500
300
350
400
450
500
x
(a) y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
x
(b) y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
x
(c) Fig. 5.4. Distribution of the mesh points in Example 5.2 in x ? y plane: (a) the arclength method; (b) the projection method; (c) the periodic method.
22
L. LIU, G. MOORE, AND R. D. RUSSELL Table 5.4
Absolute error for L at c = 8:3
NTST 10 20 40 80 160
m=2 2.16e-01 1.10e-01 7.78e-03 4.89e-04 2.94e-05
m=3 2.25e-02 8.88e-04 5.57e-05 3.24e-06 2.01e-07
m=4 4.35e-05 3.32e-06 5.51e-08 8.00e-10
m=5 m=6 1.25e-05 6.67e-07 2.62e-07 2.00e-09 5.50e-09 1.00e-10
Table 5.5
Convergence rate for L at c = 8:3
NTST m = 2 10!20 0.97 20!40 3.82 40!80 3.99 80!160 4.05
m=3 m=4 m=5 m=6 4.67 3.71 5.57 8.38 3.99 5.91 5.57 4.11 6.11 5.78 4.01
is Sil'nikov-type. To obtain the high order convergence in such a case, it would be necessary to use a more accurate approximation such as that discussed in section 3.4. Since this would require a major change in AUTO, we use a simple implementation for which the independent variable of equation (2.4) is transformed from \time" (t) to \arclength" () directly, i.e., (); ) ; 0 < < 1; (5.1) v0() = L kff ((vv(); )k where Z 1 Z t T 0 L = T ku (t)kdt; (t) = L ku0 ( )kd ; 0 0 and the eigenvector boundary conditions (2.12){(2.14) are used. Since (5.1) has no singularity in [0; 1], the system can be solved by a standard Gauss collocation method. Note that we do not use a phase condition here but rather keep the value of xed. Example 5.3. Chua's Electronic Circuit [23] x0 = a(y ? (x)); y0 = x ? y + z; z 0 = ?by where (x) = 61 (x3 ? x). For a > 1:20245, this system has a branch of Sil'nikov-type homoclinic orbits with a xed point at the origin. The stable manifold is two dimensional. We use AUTO to trace out the periodic branch starting from the Hopf bifurcation point (1; 0; ?1) at (a; b) = (3:474937; 5) and locate a periodic solution at (a; b) = (4:35164; 5) with period T = 200. We use both the periodic method and the eigenvector method to compute the homoclinic branch with xed T = 200. With NTST=50 and the Newton tolerance set to 10?8 both methods successfully continue this branch to a > 50. For the eigenvector method, we need to increase the tolerance for the parameters since small parameters ( ) are sensitive to the relative error tolerance. The estimated absolute error for b at a = 6 is less than 10?8 in both cases. When we use the equations (5.1) and (2.12){(2.14), the continuation of the homoclinic branch is successful to around a = 10. As before, the estimated absolute error for b at a = 6 is O(10?8). In this computation, NTST=50 and the error tolerances for parameters and solutions are 10?3 and 10?8. The mesh distributions for
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
23
x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2
0
0.2
0.4
0.6
0.8
1
σ
(a) x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.7
0.75
0.8
t
(b) x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.5
0.55
0.6
0.65
t
(c) Fig. 5.5. Distribution of the mesh points in Example 5.3 for x: (a) the arclength method (with eigenvector boundary conditions); (b) the eigenvector method; (c) the periodic method.
24
L. LIU, G. MOORE, AND R. D. RUSSELL
all three successful methods at a = 6 are shown in Fig. 5.5. The early failure of the eigenvector method with arclength is due to the qualitative change of the solution. When a better approximate solution at larger value of a is provided, the continuation can be carried out successfully. An interesting question is whether or not there is another phase condition or a di erent treatment near the endpoints which leads to more robustness of this simpli ed arclength approach for this Sil'nikov case, or is one of the more sophisticated approaches in section 3.4 necessary. REFERENCES [1] E. Anderson, et al., LAPACK Users' Guide, SIAM, Philadelphia, 1992. [2] U. Ascher, J. Christiansen, and R. D. Russell, Collocation software for boundary value ODE's, ACM Trans. Math. Software, 7 (1981), pp. 209{222. [3] U. Ascher, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Di erential Equations, Prentice-Hall, Englewood Cli s, NJ, 1988. [4] F. Bai, A. R. Champneys, Numerical detection and continuation of saddle-node homoclinic bifurcations of codimension one and two, Mathematics Preprint 94-04, University of Bath, 1994. [5] F. Bai, A. Spence, and A. M. Stuart, The numerical computation of heteroclinic connections in systems of gradient partial di erential equations, SIAM J. Appl. Math., 53 (1993), pp. 743{769. [6] W.-J. Beyn, Global bifurcations and their numerical computation, Continuation and Bifurcations: Numerical Techniques and Applications, D. Rossed, B. D. Dier and A. Spence eds., Kluwer, Dordrecht, 1990, pp. 169{181. [7] , The numerical computation of connecting orbits in dynamical systems, IMA J. Numer. Anal., 9 (1991), pp. 379{405. [8] , Numerical methods for dynamical systems, Advances in Numerical Analysis, W. Light, ed. Clarendon Press, Oxford, 1991-1992, pp. 175{236. [9] C. de Boor, Good approximation by splines with variable knots, Spline Functions and Approximation Theory, A. Meir and A. Sharma ed., ISNM Vol. 21, Birkhauser Verlag, Basel (1973), pp. 57{72. [10] V. Canale, The Computation of Paths of Homoclinic Orbits, Ph.D. Thesis, Department of Computer Science, University of Toronto, 1994. [11] A. R. Champneys and Yu. A. Kuznetsov, Numerical detection and continuation of codimension-two homoclinic bifurcations, Int. J. Bifurcation and Chaos, 4 (1994), pp. 785{ 822. [12] L. Dieci, J. Lorenz, and R. D. Russell, Numerical calculation of invariant tori, SIAM J. Sci. Stat. Comput., 26 (1991), pp. 607{647. [13] E. J. Doedel, AUTO: a program for the automatic bifurcation analysis of autonomous systems, Congr. Numer., 30 (1981), pp. 265{284. [14] E. J. Doedel and J. P. Kernevez, AUTO: software for continuation and bifurcation problems in ordinary di erential equations, Technical Report, Applied Mathematics, California Institute of Technology, Pasadena, CA, 1986. [15] E. J. Doedel and M. J. Friedman, Numerical computation of heteroclinic orbits, J. Comp. Appl. Math., 26 (1989), pp. 155{170. [16] M. J. Friedman and E. J. Doedel, Numerical computation and continuation of invariant manifolds connecting xed points, SIAM J. Numer. Anal., 26 (1991), pp. 789{808. [17] , Computational methods for global analysis of homoclinic and heteroclinic orbits: a case study, J. Dynamics Di . Eqns., 5 (1993), pp. 59{87. [18] J. Guckenheimer and P. Holmes, Nonlinear Oscillations Dynamical Systems, and Bifurcations of Vector Fields, Applied Mathematic Sciences Vol. 42, Springer, New York, NY, 1983. [19] F. de Hoog and R. Weiss, Di erence methods for boundary value problems with a singularity of rst kind, SIAM J. Numer. Anal., 13 (1976), pp. 775{813. [20] , Collocation methods for singular boundary value problems, SIAM J. Numer. Anal., 15 (1978), pp. 198{217. [21] , An approximation theory for boundary value problems on in nite intervals, Computing, 24 (1980), pp. 227{239. [22] H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, Proc. Applications of Bifurcation Theory, P. H. Rabinowitz, ed., Academic Press, New York,
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
25
1977, pp. 359{384. [23] A. I. Khibnik, D. Roose and L. O. Chua, On periodic orbits and homoclinic bifurcations in Chua's circuit with smooth nonlinearity, Int. J. Bifurcation and Chaos, 3 (1993), pp. 363{ 384. [24] M. Lentini and H. B. Keller, Boundary value problems over semi-in nite intervals and their numerical solution, SIAM J. Numer. Anal., 17 (1980), pp. 577{604. [25] L. Liu, Numerical methods and software for solving boundary value ODEs for bifurcation analysis, Ph.D. Thesis, Department of Mathematics & Statistics, Simon Fraser University, 1993. [26] L. Liu and R. D. Russell, Linear system solvers for boundary value ODEs, J. Comp. App. Math., 45 (1993), pp. 1{15. [27] Y. Liu, L. Liu, and T. Tang, The numerical computation of connecting orbits in dynamical systems: a rational spectral approach, J. Comp. Phys., 111 (1994), pp. 373{380. [28] G. Moore, Computation and parametrisation of periodic and connecting orbits, IMA J. Numer. Anal., 15 (1995), pp. 245{263. [29] G. Moore and E. Hubert, Algorithms for computing stable manifolds of stationary solutions, submitted to IMA J. Numer. Anal. [30] S. Schecter, Numerical computation of saddle-node homoclinic bifurcation points, SIAM J. Numer. Anal., 30 (1993), pp. 1155{1178. [31] , Rate of convergence of numerical approximations to homoclinic bifurcation points, IMA J. Numer. Anal., 15 (1995), pp. 23{60. [32] L. P. Sil'nikov, A case of the existence of a denumerable set of periodic motions, Sov. Math. Dokl., 6 (1965), pp 163{166. [33] E. Weinmuller, A di erence method for a singular boundary value problems of second order, Math. Comp., 42 (1984), pp. 441{464. [34] , On the boundary value problem for systems of ordinary second order di erential equations with a singularity of the rst kind, SIAM J. Math. Anal., 15 (1984), pp. 287{307. [35] , Collocation for singular boundary value problems of second order, SIAM J. Numer. Anal., 23 (1986), pp. 1062{1095. [36] , Stability of singular boundary value problems and their discretization by nite di erences, SIAM J. Numer. Anal., 26 (1989), pp. 180{213.
Abstract. In this paper, we study a numerical method for the computation and continuation of homoclinic and heteroclinic orbits based upon the arclength parameterization of the orbits. Unlike most other methods, this method utilizes the geometric structure of the homoclinic and heteroclinic orbits and does not require solving a boundary value problem on an in nite interval. However, the boundary value problem formulated by this method can have a singularity at the end of the domain, and thus we introduce a special collocation method to handle such a singularity. We discuss the convergence properties of our collocation method and the implementation of the method which uses the software AUTO. For several examples we show that the arclength parameterization compares very favorably with the other numerical methods, although there are some limitations in the Sil'nikov case. Key words. homoclinicorbits, heteroclinic orbits, collocationmethods, numerical continuation. AMS subject classi cations. 34C37, 65L10, 65L60.
1. Introduction. Numerical computation of homoclinic and heteroclinic orbits is of interest in a variety of contexts. For example, structural changes in dynamical systems are often related to the appearance or disappearance of trajectories connecting one or several stationary points. Homoclinic orbits often arise as the limiting case of periodic solutions while heteroclinic orbits can represent travelling wave solutions of parabolic partial di erential equations. In the global bifurcation analysis, the occurrence of homoclinic orbits is closely tied to chaotic behavior in dynamical systems [18]. Consider the parameterized ODE system (1.1) u0(t) = f (u(t); ); ?1 < t < 1 where u 2 Rn and 2 Rp . For a given , if there exists a non-constant solution u(t) of (1.1) such that (1.2) u? = t!?1 lim u(t); u+ = tlim !1 u(t); f (u ; ) = 0; then the pair (u(t); ) is called a connecting orbit between stationary solutions or simply a connecting orbit. Moreover, if u? = u+ , the orbit is called a homoclinic orbit; otherwise, it is called a heteroclinic orbit. We assume both u+ and u? are hyperbolic xed points. Let MS and MU be the stable and unstable manifolds of these two stationary solutions de ned by t lim t(u) = u g; MS = fu 2 Rn : tlim MU = fu 2 Rn : t!?1 !1 (u) = u g; Received by the editors July 1, 1995; accepted for publication (in revised form) December 20, 1995. y Department of Computer Science, Concordia University, Montreal, Qecbec, Canada H3G 1M8. Current address: Nortel Technology, P.O. Box 3511, Station C, Ottawa, ON, Canada K1Y 4H7 ([email protected]). z Department of Mathematics, Imperial College, Queen's Gate, London, UK SW7 2BZ ([email protected]). x Department of Mathematics & Statistics, Simon Fraser University, Burnaby, BC, Canada V5A 1S6 ([email protected]). 1
2
L. LIU, G. MOORE, AND R. D. RUSSELL
with dimensions dim(MU? ) = n? ; dim(MS? ) = n ? n? ; dim(MS+ ) = n+ ; dim(MU+ ) = n ? n+ ; where t is the ow corresponding to (1.1). Note that when MU? and MS+ intersect, MU? \ MS+ is at least one dimensional. Moreover, if n? + n+ = n + 1, then generally there exists exactly one orbit between u? and u+ . When n? + n+ > n + 1, dim(MU? \ MS+ ) > 1, and further conditions are required to obtain the unique connecting orbit (e.g. see [18]). When n? + n+ < n + 1, the connecting orbit is not structurally stable, and we need to add some free parameters to establish the stable connection between u? and u+ . In this case, the number of free parameters is (1.3)
p = n + 1 ? (n? + n+ ):
In particular, for homoclinic orbits p = 1, or n? +n+ = n. A connecting orbit (u(t); ) between stationary solutions u? and u+ is non-degenerate if both are hyperbolic xed points, A = t!1 lim f u (u(t); ) = f u (u ; ); p = n + 1 ? (n? + n+ ); and the ODE system
v0 = f u(u; )v + f (u; ) has only the trivial solution = 0, v = cu0 for some c 2 R [7]. Throughout the
(1.4)
paper, we assume the connecting orbits are non-degenerate. In section 2, we review the methods developed in [6], [7], [14], [15], [16] and [17] for the computation of connecting orbits. These are all based upon solving boundary value problems (BVPs) in a truncated nite interval [T?; T+ ] (and usually this interval is rescaled to [0; 1]), so one can apply standard BVP solvers such as collocation methods or multiple shooting methods to compute the solutions. One drawback is that the derivative of the solution is usually quite large in a narrow region, and the BVP is similar to a singular perturbation problem. In section 3, we discuss an arclength parameterization method for the computation of homoclinic and heteroclinic orbits and their solution branches using continuation. The main advantage of this method is that the BVP is de ned on a nite interval using this natural \phase space scaling" (the total arclength). It thus avoids truncating and rescaling the in nite interval. Moreover, the solution is smooth over the whole domain, except possibly at the end points. When f u (u; ) at u? =u+ has dominant real positive/negative eigenvalues, this method has a high rate of convergence and has distinct computational advantages over the previous methods (allowing, e.g., the possibility to trace saddle node homoclinic orbits of very large length). In the case of Sil'nikov-type problems, the simple use of arclength implemented here has some distinct limitations, although we discuss more sophisticated alternative implementations which warrant pursuing in the future. In [25] the arclength parametrization method has been implemented in the well-known software package AUTO [14]. In section 4 we give a brief discussion of this implementation, which requires several signi cant changes to AUTO itself. Numerical examples which demonstrate the robustness of the approach are presented in section 5.
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
3
2. Numerical Methods for Computing Connecting Orbits. The strategy used in numerical methods for computing connecting orbits can usually be summarized as solving the boundary value problem (2.1) u0 ? f (u; ) = 0; (2.2) B(u(T? ); u(T+ ); ) = 0; (2.3) J (u; ) = 0; on an interval J = [T? ; T+ ], T? < 0, T+ > 0. Here, B(u(T? ); u(T+ ); ) are the boundary conditions and J (u; ) is a so-called phase condition. The methods are generally distinguished from one another by the way in which di erent types of boundary conditions and phase conditions are introduced. For convenience, we consider the BVP in the normalized interval t 2 [0; 1], so (2.4) u0 ? T f (u; ) = 0; (2.5) B(u(0); u(1); ) = 0; (2.6) (u; ) = 0; where T = T+ ? T? . In general, a robust choice of the phase condition is the integral phase condition (2.7)
(u; ) :=
or (2.8)
(u; ) :=
1
Z
0
1
Z
0
[u(t) ? u0 (t)]T u0 (t)dt = 0
[f (u; ) ? f (u0 ; 0)]T f u (u; )f (u; )dt = 0
where u0(t) is a nearby known solution (usually the solution obtained from the previous continuation step). These conditions can be obtained by minimizingthe functional (2.9)
E() :=
or (2.10)
G() :=
1
Z
0 Z
0
1
ku~ (t + ) ? u0 (t)k2dt ku~ 0 (t + ) ? u00(t)k2 dt
against the phase shift . The phase condition (2.7) is used in [6], [7] and [14] while (2.8) is used in [15]{[17]. The boundary conditions (2.5) can be dealt with in a variety of ways. 2.1. Periodic Boundary Conditions. The periodic boundary conditions are given by (2.11) u(0) = u(1): In many applications, homoclinic orbits occur as limiting cases of periodic solutions. It is natural to use these periodic solutions to approach these connecting orbits, and it can in fact be shown that there exists a branch of periodic solutions near a homoclinic orbit [6]. The computation of a connecting orbit this way is usually achieved
4
L. LIU, G. MOORE, AND R. D. RUSSELL
by the continuation of a branch of periodic solutions where the period T tends to in nity when the connecting orbit is approached. This method is basically available for computing homoclinic orbits in the software AUTO. There is one free parameter in the system, which means that p = 1 and n? +n+ = n. This condition always holds for homoclinic orbits. One disadvantage of the periodic boundary condition approach is that it is not generally suitable for the computation of heteroclinic orbits. 2.2. Eigenvector Boundary Conditions. The eigenvector boundary condition is introduced in [15]{[17]. If u are stationary solutions of (1.1),
f (u?; ) = 0;
f (u+ ; ) = 0: ? , and (u ; ) has If f u(u? ; ) has n? eigenvalues with positive real parts, f?igni=1 fu + n + n+ eigenvalues with negative real parts, f+i gi=1 , satisfying (2.13) f u(u; )vi = ivi ; kvik = 1; i = 1; : : :; n;
(2.12)
then the eigenvector boundary conditions can be de ned as (2.14) u(0) ? u? = ?
n
? X
i=1
c?i v?i ; u(1) ? u+ = +
n
+ X
i=1
c+i v+i ;
n
X
i=1
c2i = 1
are constants. where and fci gni=1 The geometric interpretation of the eigenvector boundary conditions is to require that the endpoint u(0)=u(1) of the orbit u(t) lie on the unstable/stable tangent space of u? =u+ and be a certain distance from the xed point u? =u+ . These tangent ? , and fv gn+ , subspaces, denoted as E?U () and E+S (), are spanned by fv?i gni=1 +i i=1 respectively. When the orbit is non-degenerate, there are (n + 2)(nn? + n+ ) + p + 2 = , and ), n (n + 1)(n? + n+ ) + n + 3 unknown parameters (fi; vi; ci gi=1 di erential equations and (n+1)(n? +n+ )+2n+3 constraints. A more straightforward as orthonormal bases for the tangent implementation would be to choose fvigni=1 S U subspaces E? () and E+ (). 2.3. Projection Boundary Conditions. An alternative (but mathematically equivalent) approach uses the projection boundary conditions (also known as asymptotic boundary conditions) [6], [7]. The boundary conditions again take the form
(2.15)
u(0) ? u? 2 E?U (); u(1) ? u+ 2 E+S ();
but for numerical calculations the projection boundary conditions are often presented in the form (2.16) B? (u(0); ) := P? ()[u(0) ? u? ] = 0; B+ (u(1); ) := P+ ()[u(1) ? u+ ] = 0 where the rows of the matrices P? () 2 R(n?n?)n and P+ () 2 R(n?n+)n span the stable subspace of f Tu (u? ; ) and the unstable subspace of f Tu (u+ ; ), respectively. Note that (2.15) and (2.16) are equivalent when u are both hyperbolic xed points. Equation (2.16) de nes 2n ? (n? + n+ ) boundary conditions. The BVP (2.4){(2.6) is well-posed if the connecting orbit is non-degenerate, as the number of free parameters is p = n + 1 ? (n? + n+ ). In general, one cannot handle non-hyperbolic xed points without some modi cation. However, Schecter [30], [31] shows that under certain conditions the projection
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
5
boundary conditions can be modi ed to handle certain types of connecting orbits with semi-hyperbolic xed points (quadratic turning points). The equivalent treatment for the eigenvector boundary condition approach is in [17]. The computation of fully non-hyperbolic homoclinic orbits is studied in [4], [10]. Based on the projection boundary conditions, Champneys and Kuznetsov [11] develop a method for the detection and continuation of several di erent types of codimension-two homoclinic bifurcation points. For hyperbolic xed points, these types of boundary conditions all have the property that convergence to the connecting orbit is exponential with respect to the length of the truncated interval. The eigenvector and projection boundary conditions have a better rate of convergence than the periodic boundary conditions, and they can be used to handle general types of heteroclinic orbits (p 6= 1). Both approaches are similar to the one used in [24] for solving BVPs over in nite intervals, where the boundary conditions are obtained from linear approximations of the stable and unstable subspaces of the xed points. However, because of their di erent implementations, each has its own advantages and disadvantages. For the eigenvector boundary condition approach, the dimension of the original problem is increased because of the introduction of the unknown parameters fi; vi ; cig and . Each requires certain knowledge of the invariant subspace structure at the xed points, one a smooth basis and the other a smooth projection. A smooth QR factorization is desirable in such a case (e.g. see [7]). Finally, to perform numerical continuation for the connecting orbit branches, we add a pseudo-arclength continuation equation [22] (2.17)
2
u
Z
0
1
[u(t) ? u0 (t)]T u_ 0 (t)dt + 2 ( ? 0)T _ 0 + 2 ( ? 0 )_ 0 = s
to (2.4){(2.6), where u , and are scaling parameters, (u0 ; 0; 0 ) is the solution at a previous continuation step, (u_ 0 ; _ 0; _ 0 ) denotes the derivative along the solution branch, and s is the continuation stepsize. The number of free parameters in such a situation is n + 2 ? (n1 + n+ ) (an extra parameter is introduced for continuation purposes). The continuation equation (2.17) is used for solving BVPs in AUTO. 3. Arclength Parameterization. All three methods discussed in the previous section have the disadvantage that when T becomes very large, u0 (t) varies fast (O(T)) in only a small part of the domain [0; 1], and on the rest of the domain u0 (t) 0. In fact, the BVP system (2.4){(2.6) is similar to a singular perturbation problem, where the existence of an interior layer or boundary layer can cause great numerical diculty. However, in phase space the connecting solution itself often remains bounded and smooth except at the xed points. This suggests that a geometrical formulation could lead to a superior method. Moore [28] has recently introduced such a geometrical approximation based upon an arclength parameterization. 3.1. Parameterization and Continuation. Let L be the total arclength of the solution u(t) in phase space, i.e. (3.1)
L=
Z
1
?1
ku0(t)kdt;
and let be the scaled arclength parameterization variable Z t 1 (3.2) ku0( )kd : (t) = L ?1
6
L. LIU, G. MOORE, AND R. D. RUSSELL
Then, if t() denotes the inverse of (t), (3.3) v() := u(t()) satis es the arclength parameterized di erential equation (); ) ; 0 < < 1: (3.4) v0() = L kff ((vv(); )k Note that the endpoints are now singular points for (3.4), and away from them kv0()k = L. Assuming that the connecting orbit is not of Sil'nikov type, i.e., the dominant eigenvalue of f u (u ; ) is real (the dominant eigenvalue for f u (u? ; ) is that with smallest positive real part, while the dominant eigenvalue for f u (u+ ; ) is that with largest negative real part), 0 lim v0 () and lim !0 !1 v ()
exist and lie in the dominant eigenspaces of u? and u+ , respectively. The continuity conditions require (3.5) kv0 (0)k = kv0 (1)k = L; and the boundary conditions (3.6) v(0) = t!?1 lim u(t) = u? ; v(1) = tlim !1 u(t) = u+ ensure that the xed points u? and u+ are the endpoint values for the solution. In this formulation a phase condition is not needed, as the endpoints are xed and L is bounded. The projection boundary conditions (2.16) become (3.7) P?()v0 (0) = 0; P+ ()v0 (1) = 0 where the projection matrices P? () 2 R(n?n? )n and P+ () 2 R(n?n+)n are de ned as in section 2. Moore proves that the BVP system (3.4){(3.7) is well-posed when the connecting orbit is non-degenerate [28]. A pseudo-arclength continuation equation (3.8) v2
Z
0
1
[v() ? v0 ()]T v_ 0()d + 2 ( ? 0)T _ 0 + L2 (L ? L0 )L_ 0 = s
may be used to perform numerical continuation for a branch of connecting orbits. The total number of free parameters (; L) for the system (3.4){(3.8) is n+2 ? (n? +n+ ). 3.2. Numerical Discretization. The BVP (3.4){(3.8) consists of n ODEs, one integral equation, 4n+2 ? (n? +n+ ) algebraic constraints, and p+1 parameters. Thus, this system has 2n extra conditions. Also, the BVP has singularities at = 0 and 1, where the right hand side of (3.4) is unde ned. Consequently, a standard boundary value solver is not suitable for use without some modi cation. Moore [28] recommends solving the BVP (3.4){(3.7) by a natural Gauss-Lobatto collocation method: Use a C 1 piecewise polynomial of degree m and a mesh 0 = 0 < < N = 1; Collocate (3.4) at all Lobatto points apart from = 0; 1; Require (3.5): kv0(0)k = kv0 (1)k = L; Require (3.6): v(0) = u? (); v(1) = u+ ();
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
Require (3.7): P? ()v0(0) = 0; P+ ()v0(1) = 0.
7
One of our goals has been to modify this algorithm for the software AUTO, which uses Lagrange polynomials as the basis functions and Gauss collocation points (and consequently involves di erent continuity conditions). Its modi cation to use Lobatto points would be dicult, so our implementation of the arclength formulation in AUTO uses Gauss collocation points in the interior subintervals, and for consistency between the equations and unknowns, Radau collocation points on the rst and last subintervals (see [25], [28]). The resulting numerical algorithm is the Gauss-Radau collocation method: Use a C 0 piecewise polynomial of degree m and a mesh 0 = 0 < < N = 1; Collocate (3.4) at m Gauss points fik gmk=1 in [i; i+1] for i = 1; :::; N ? 2; Collocate (3.4) at the left m ? 1 Radau points f0kgmk=2 in the rst subinterval [0; 1]; ?1 in the last subin Collocate (3.4) at the right m ? 1 Radau points fN ?1;k gmk=1 terval [N ?1 ; N ]; Require (3.5): kv0(0)k = kv0 (1)k = L; Require (3.6) or equivalently: f (v(0); ) = 0; f (v(1); ) = 0; Require (3.7): P? ()v0(0) = 0; P+ ()v0(1) = 0. For numerical continuation, the algorithm is supplemented with the continuation equation (3.8). A more detailed description of the implementation for AUTO is in the next section. 3.3. Convergence Results. The system (3.4){(3.7) is a nonlinear boundary value problem with a singularity of the rst kind at the endpoints. This type of singular BVP has been studied by de Hoog and Weiss [19], [20], and Weinmuller [33]{[36]. We now discuss the convergence properties of the Gauss-Radau collocation algorithm for solving (3.4){(3.7). The variational equation for (3.4) at the solution (v(); L; ) is dy ? L[I ? P()] f u (v; ) y = 0; 0 < < 1 (3.9) d kf (v; )k where 0()v 0()T f f T (3.10) P() = vv0() T v0 () = f T f is an orthogonal projection onto v0 (). If we assume that the dominant eigenvalues of the matrices f u (u ; ) are real and
v0 (0) := lim v0 () and v0 (1) := lim v0() !0+ !1?
exist and lie in their dominant eigenspaces, then we can rewrite (3.9) as y d M M 0 1 y (3.11) d ? + 1 ? ? B() = 0 where L[I ? P()]f u = 1 [I ? P(0)] (u ; ); (3.12a) M0 = lim fu ? !0 kf k ~1
8
L. LIU, G. MOORE, AND R. D. RUSSELL
(3.12b)
~1 = v0 (0)T f u (u? ; )v0(0)=L2 ;
(3.13a)
(1 ? )L[I ? P()]f u = 1 [I ? P(1)] (u ; ); M1 = lim fu + !1 kk ~
(3.13b) and
f
1
~1 = v0(1)T f u (u+ ; )v0 (1)=L2; B() = ?L[I ? P()] kf u((vv;;))k + M0 + 1M?1 :
f
We approximate y() by p(), a C 0 piecewise polynomial of degree m, satisfying the collocation equations 0 + M1 ? B( ) p( ) = 0; p0 (ik ) ? M (3.14) ik ik ik 1 ? ik where 0 i N ? 1; 1 k m, and ik = i + k hi are the collocation points in the interval [i ; i+1]. Note that for the Gauss-Radau collocation algorithm, (3.14) is not well-de ned at 0;1 = 0 and N ?1;m = 1, and these equations are replaced by the projection boundary conditions (3.15) P? ()p(0) = 0; P+ ()p(1) = 0: To simplify our analysis for the BVP (3.11), (3.15), let 2 0 21 = 2(1 ? ) 21 1 and ( ) 0 21 ; z( ) = y0( ) ; B() = B0 ( ) 0 21 : y() = yy0( ) 1 1 y1( ) B1 ( ) 12 1 1 2 Then (3.11) becomes dz 1 (3.16) d ? M z + D( )z = 0 where # " M1 ? B0 ( ) M 0 2 ? 2 M= ?M1 ; D( ) = ? 2M?0 + B12( ) : From (3.15), the boundary conditions become P? () (3.17) z (0) = 0 ; I ? I z(1) = 0: P+ () Our next result follows immediately from [19] Theorem 3.5. Lemma 3.1. Let D( ) 2 C p[0; 1], p 0. Then the BVP (3.16) and (3.17) has a unique solution z( ) such that
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
9
(i) z( ) 2 C p+1 (0; 1]. (ii) z( ) 2 C p+1 [0; 1] when all eigenvalues of M0 have non-positive real parts and all eigenvalues of M1 have non-negative real parts. (iii) Let be the eigenvalue of M0 with smallest positive real part, be the eigenvalue of M1 with largest negative real part and
= minfRe( ); ? Re( )g: Then
kz(p+1) ( )k const ?p?1 (j log jd ?1 + 1) kz(p+1) ( )k const (j log jd + 1) z( ) 2 C p+1 [0; 1]
p< < p+1
= p+1
> p+1 where d is the size of the largest Jordan block of M associated with or : The eigenvalues of M0 , M1 in (3.11) and the constants , in Lemma 3.1 have a close relationship with the eigenvalues of f u (u ; ). In fact, ordering the eigenvalues of f u(u? ; ) with positive real parts figki=1 such that 1 Re(2 ) Re(3 ) Re(k ), v0(0) is the approximate dominant eigenvector of f u (u? ; ), so from (3.12), ~1 1 . Moreover, v0(0) is an eigenvector of M0 with associated eigenvalue 1(M0 ) = 0. Other eigenvalues of M0 with positive real parts are i (M0 ) = i =~1 i =1 for 1 < i k, which implies 2=1 . Similarly, 1(M1 ) = 0, i (M1 ) ?i =1 for 1 < i l, and ?2 =1, where the eigenvalues of f u (u+ ; ) with negative real parts are ordered such that 1 Re(2) Re(3 ) Re(l ). Therefore,
min Re 2 ; Re 2 ; 1 1 and from Lemma 3.1, the smoothness of the solution z() depends on the size of . The analogue of (3.14) for solving (3.11) is the collocation method (3.18) p0 ( ik ) ? 1 M p( ik ) + D( ik )p( ik ) = 0; 0 i N ? 1; 1 k m ik for solving (3.16) where ik = i + k hi are the collocation points in [ i; i+1 ]. Singularity only occurs at = 0, and the collocation equations at this point are replaced by P?() p(1) = 0: p (3.19) (0) = 0 ; I ? I P+ () Theorem 3.2. Let D( ) 2 C m [0; 1]. For h = max0iN ?1fhi g suciently small, the solution p( ) of (3.18) and (3.19) satis es 8 const h j loghjd ?1; if 1 < < m; > > < m j loghjd ?1 ; const h if m < m + 1; kz( ) ? p( )k > const hm j loghjd ; if = m + 1; > : const hm ; if > m + 1 or there is no positive : Proof. Here we give the proof for the uniform mesh case; the extension to the non-uniform case can be obtained following the approach in [35]. (i) Let 1 < < m + 1. From Lemma 3.1, kz(m+1) ( )k const ?m?1 (1 + j log jd ?1 ):
10
L. LIU, G. MOORE, AND R. D. RUSSELL
De ne a degree m piecewise polynomialerror function e( ) satisfying the interpolation conditions e0( ik ) = z0 ( ik ) ? p0( ik ): Then, for 2 [ i ; i+1], m X ? ik 0 z0( ik ) ? p0 ( ) e ( ) = Lk h k=1 m X = z0 ( ) ? p0 ( ) + Lk ?h ik z0( ik ) ? z0( ) k=1 0 0 p = z ( ) ? ( ) + O(hm ?m?1 (1 + (log )d ?1)) where the Lk are the Lagrange basis functions. After integration, e( ) = z( ) ? p( ) + O(hm ?m (1 + (log )d ?1 )): It follows that e( ) satis es (3.20) e0( ik ) ? 1 M e( ik ) + D( ik )e( ik ) = O(hm ik ?m?1 (log ik )d ?1) ik for 0 i N ? 1; 1 k m, and the boundary conditions P ( ) ? (3.21) I ?I e(1) = 0: P+ () e(0) = 0; In a similar manner as that in the proof of Lemma 3.8 in [35], we have kz( ) ? p( )k const hm j loghjd ?1 ; = m; kz( ) ? p( )k const hm j loghjd ?1 j1 ? ?m j const hm j loghjd ?1 > m; const h j loghjd ?1
< m: (ii) If = m + 1, then kz(m+1) ( )k const (1 + j log jd ); and the right hand side of (3.20) is replaced by O(hm (1 + (log )d ), which gives kz( ) ? p( )k const hm j loghjd : (iii) If > m + 1, then z 2 C m+1 [0; 1] and the right hand side of (3.20) is replaced by O(hm ). Therefore, kz( ) ? p( )k const hm : Finally, the proof for the non-positive case is given in Theorem 4.5 of [20]. One advantage of the collocation method when applied to standard boundary value problems (without singular points) is the property of superconvergence for suitably chosen 1; ; m . If the solution is suciently smooth and (3.22)
Z
0
1
sk w(s)ds = 0;
k = 0; ; < m;
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
then
ky(i ) ? p(i )k const hm+ +1 ;
11
i = 0; ; N: In particular, for the Gauss points, the rate of convergence is O(h2m ); for the Lobatto points, it is O(h2m?2); and for the Radau points, the superconvergence rate is O(h2m?1). In fact, applying the similar technique from [20] Theorem 5.1 and [35] Theorem 4.1, the rate of convergence in Theorem 3.2 can be improved to O(hm+1 ) provided D( ) 2 C m+2 [0; 1], > m + 2 and (3.22) holds. This rate is most likely optimal (see examples given in [20], [35]). It should be noticed that the above convergence results are mainly dependent upon the continuity conditions of the solution z( ), which is in turn determined by the ratio of the two smallest eigenvalues (in magnitude) of f u (u ; ). In general, when z( ) does not satisfy the continuity conditions, the optimal rate of convergence cannot be obtained. In [9], de Boor shows that a non-smooth function such as f(t) = t with 2 (0; 1) can be approximated by a piecewise polynomial of order k with accuracy O(hk ) if jf (k) (t)j decreases monotonically in (0; 1] and the mesh points are given by the equidistribution method Z 1 Z ti+1 1 ( k ) 1 =k jf (t)j dt = N jf (k) (t)j1=kdt; i = 0; 1; ; N ? 1: (3.23) ti 0 In AUTO, the equidistribution method (3.23) is used for k = m + 1 where m is the number of collocation points per subinterval. In general, each jzi(m+1) ()j may not be monotonically decreasing on (0; 1). However, for the non-Sil'nikov case, all jzi(m+1) ()j decrease monotonically in (0; ] and increase monotonically in [1 ? ; 1) for suciently small . Therefore, with the proper mesh selection scheme (base on (3.23)), the collocation method can be expected to have a rate of convergence at least O(hm+1 ) regardless of the continuity of the solution z() at the end-point. Our numerical examples presented in section 5 indicate that this is the case. 3.4. S il'nikov-type Connecting Orbits. For the Sil'nikov case v0(0) and/or 0 v (1) do not exist. As a result, the arclength parametrization is more complicated in principle, and practical results using a straightforward implementation have been less good. However, we believe that even then, using arclength away from the xed points is an advantageous procedure. Therefore, in this section, we show how the arclength idea may be adapted to this more dicult situation. For de niteness we assume that the Sil'nikov behaviour occurs at u? but not at u+ , although the algorithms below generalize easily. In order to avoid the singularity at u? (), we introduce a boundary condition involving the unstable manifold of u? () rather than u? () itself. Thus we have the di erential equation (3.4), i.e. (); ) ; 0 < < 1; (3.24) v0() = L kff ((vv(); )k with boundary conditions (3.25) v(0) 2 MU? at a \distance" ? from u? (); v(1) = u+ () and end conditions (3.26) P+ v0(1) = 0; kv0(1)k = L:
12
L. LIU, G. MOORE, AND R. D. RUSSELL
Here ? > 0 is a small user-chosen constant, which is analogous to the choice of T? in section 2. This set of equations may be discretized by Gauss collocation on all subintervals apart from the last, where Gauss-Radau collocation should be used with (3.26) replacing the di erential equation at = 1. The remaining question is how to make the boundary condition at = 0 practical and, in particular, specify what we mean by \distance". There are two possibilities. 1) Using the unstable subspace E?U (). Here we simply mimic (2.14) and take v(0) to be in the unstable subspace of u? (), i.e. (3.27)
v(0) ? u?() = ?
n? X i=1
c?iv?i ;
n? X i=1
c2?i = 1:
The di erence is that now ? is known and only c?1 ; : : :; c?n? are to be determined. The size of ? controls the O(2?) error introduced by approximatingthe manifold with ? should be chosen the subspace. Moreover, to avoid possible ill-conditioning, fv?i gni=1 as a well-conditioned basis for E?U () for which the restriction fu(u? (); ) to E?U () is positive-de nite, as simply choosing an orthonormal basis is usually unacceptable [29]. 2) Approximating the manifold more accurately. In [29] a number of algorithms are described for computing the unstable manifold of a given xed point. For our connecting orbit problem, we think the potentially most useful approach is to solve the system (3.28) w0( ) = 1 f (w( ); ); 0 < < 1 ? with boundary condition
w(0) = u?();
(3.29) end condition
P? ()w0 (0) = 0;
(3.30) and integral side-constraint (3.31)
Z
1 kf (w( ); )k
d = ? :
? Now ? measures distance in terms of arclength, and ? > 0 is a user-chosen constant which corresponds to the explicit exponential change-of-variable 0
= e? t ; w( ) u(t) (see [29]). Hence, choosing ? < min1in? fRe(?i )g (notation as in (2.13)) forces w0 (0) to be zero. Equation (3.28) may be discretized by Gauss collocation on all subintervals except the rst, where Gauss-Radau is used with (3.30) replacing the di erential equation. This is also the natural quadrature rule for (3.31), making use of the fact that the integrand is zero at = 0. Note that (3.28) should not require many subintervals for accurate approximation. Finally (3.28){(3.31) are coupled to (3.24){(3.26) through the continuity condition (3.32) v(0) = w(1):
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
13
To conclude, we have two techniques for dealing with Sil'nikov behaviour or, more generally, with v() having a low degree of smoothness at an endpoint. The rst is a direct analogue of those described in sections 2.2 and 2.3 (which could, in fact, be used at all endpoints). The second, although requiring further work to establish the optimal choices of ?; ? and the mesh for (3.28), warrants further investigation as a more sophisticated approach to solve this dicult class of problems. 4. Software Implementation. We now describe our implementation of the Gauss-Radau collocation algorithm of section 3 in the software AUTO. For this algorithm there are only (m ? 1)n collocation equations in the rst and the last subintervals, so n equations are needed to replace the collocation equations at = 0 in the rst subinterval and at = 1 in the last subinterval. To do this, we use the equation (4.1) f (v(0); ) = 0 to replace the collocation equation at = 0, and (4.2) f (v(1); ) = 0 at = 1. Equations (3.5), (3.7) and (3.8) are treated as boundary conditions, and these are placed in the bottom rows of the collocation matrix. Here, v0 (0) and v0(1) are approximated by mth order Lagrange interpolation for non-Sil'nikov case. The projection matrices P?() and P+ () are obtained by using the LAPACK [1] subroutine DGEES. The projection matrices are not unique, and care must be taken to maintain their smoothness. The nal structure of the linear system is shown in Fig. 4.1(a) for the case of 2 di erential equations, 4 mesh intervals, 3 collocation points and 3 free parameters. When solving a linear system arising from discretization of a BVP, the rst step in AUTO is to use a local condensation method to eliminate the unknowns at non-mesh points. This process uses no row pivoting, so in our case reordering of the equations in the rst block is necessary to avoid breakdown. We place the end conditions (4.1) at the bottom of the rst block, and the coecient matrix pro le is then as shown in Fig. 4.1(b) [25]. The condensed linear system has the (almost) block bi-diagonal structure 3 2 A1 C1 D1 6 A2 C2 D2 77 6 6 .. 77 . . .. .. (4.3) 6 . 7 6 4 AN CN DN 5 B0 B1 BN ?1 BN E where the blocks C1 and AN are zero matrices. Numerical methods for solving linear systems like (4.3) are studied in [26]; in general, either a QR factorization or LU decomposition with row pivoting can be applied to eliminate the interior mesh points. After such an elimination process, the coecient matrix is reduced to the form 2 0 R1 F1 H1 3 6 0 R2 F2 H2 77 6 6 .. 77 . . . .. .. .. 6 . 7: (4.4) 6 6 GN ?1 0 F H N ? 1 N ?1 7 6 7 4 0 C~N D~ N 5 B0 B1 BN ?1 BN E
14
L. LIU, G. MOORE, AND R. D. RUSSELL r r r r r r
r r r r r
r r r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r r r r r r r
r r r r r
r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r r r
r r r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r
r r r r r r r r r r
r r r r
r r r r
r r r r
r r r r
r r r r r r r r r r r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r r r r r r r
r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
r r r r r r r r r r r r r r r r r r r r r r r r r r r r r
(a) r r r r r r
r r r r r
r r r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r r
r r r r
r r r r
r r r r r r
r r r r r r
r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r r r
r r r r r r r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r
r r r r r r r r r r
r r r r r r r r r r
r r r r
r r r r
r r r r
r r r r
r r r r r r r r r r r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r
r r r r r r r r r r r
(b) Fig. 4.1. Linear system structure: (a) before reordering the equations; (b) after reordering the
equations
The blocks B1 ; ; BN ?2 can be eliminated using R1; ; RN ?2. However, if not further modi ed, the linear system solver of AUTO will fail when eliminating the block BN ?1 . It is sucient to reorder the equations and unknowns in such a way that the linear system has the form 3 2 A2 C2 D2 6 A3 C3 D3 77 6 6 .. 77 . . .. .. 6 . 7: (4.5) 6 6 0 C D N N 7 7 6 4 0 A1 D1 5 B1 B2 BN B0 E Because v(0) and v(1) are hyperbolic xed points, the matrices A1 and CN are
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
15
non-singular, and the linear system can be solved by the QR factorization or LU decomposition as before. The nal reduced linear system is 2 3 G2 R2 F2 H2 6 G3 R3 F3 H3 77 6 6 . .. 77 . . .. .. 6 .. . 7: (4.6) 6 6 0 R 0 H N N 7 6 7 4 0 A N D1 5 B1 B2 BN B0 E In order to be able to use AUTO with the arclength implementation, we have also modi ed various parts of the code to allow users to provide discrete solution data for the starting solution, and we have implemented the arclength formulation for continuation of periodic solution branches. More speci cally, the users can switch from the \time" formulation to the \arclength" formulation directly. However, several quantities, such as the solution norm, then have di erent interpretations, and the output data must be interpreted judiciously. Also, we should mention that our current software implementation is not suitable for solving large dimensional problems (n > 20). Some major changes in the basic data structure and subroutine interdependency structure in AUTO are necessary before this implementation can be realized. 5. Numerical Examples. In this section, we consider several numerical examples and show the eciency of the arclength parameterization method by comparing it to the methods given in section 2. We refer to the arclength parameterization method as simply the arclength method, and we call the (time parametrized) implementations of periodic boundary conditions, eigenvector boundary conditions and the projection boundary conditions as the periodic method, eigenvector method and projection method, respectively. In all of our computations, we use the software AUTO, having made the necessary modi cations of it for each method. Therefore, the BVP is solved by a spline collocation method with Lagrange basis functions. In all cases, we choose 4 collocation points on each mesh subinterval (the default for AUTO). The numerical tests are performed on a Silicon Graphics Indigo running IRIX 4.0.5F with a MIPS FORTRAN 77 compiler with the \?g" option for debugging purposes. Example 5.1. The Nagumo Equations u01 = u2; u02 = cu2 ? u1 (1 ? u1 )(u1 ? a);
0 < a < 1:
This system has been used in [7], [15] and [27] as a test problem for computing heteroclinic orbits (a travelling wave solution branch) between the xed points (0; 0) and (1; 0). The eigenvector method and the arclength method are tested for this problem. For the eigenvector method, we use 10 mesh intervals (NTST=10). As in [15], we start with the exact solution p exp( 22pt) ; u2(t) = u01(t) u1(t) = 2 1 + exp( 2 t)
at (a; c) = ( 12 ; 0) and T = 5 and follow the solution by increasing T to 1000. Fixing T at this point and using a as the primary continuation parameter and c is thepsecondary continuation parameter, we compute both solution branches for c = 2(a ? 12 ).
16
L. LIU, G. MOORE, AND R. D. RUSSELL Table 5.1
Relative error for b at c = 14
NTST 12 25 50 100
Projection method Periodic method Arclength method { { 8.05e-05 1.42e-05 1.69e-05 4.59e-07 8.31e-09 1.36e-09 6.09e-09 < 1.0e-11 1.58e-09 8.08e-11
The eigenvector method fails after only 6 steps. After increasing the number of mesh subintervals to NTST=15, we successfully compute both branches. The absolute error for the parameter c in the lower half branch varies from O(10?8) to O(10?4). For the arclength method, we successfully compute both branches with only 10 subintervals (NTST=10), and the absolute error for c is O(10?8) for all solutions (see Fig. 5.1(a)). In this computation, the tolerance for the Newton iteration is set to 10?6. The distribution of the mesh points for the lower half branch solution is shown in Figs. 5.1(b) and (c) for the arclength and eigenvector method, respectively. Note that there are 41 points for the arclength method (10 mesh subintervals and 4 collocation points per interval) and 61 points for the eigenvector method (15 mesh subintervals and 4 collocation points per interval). In this and the other examples, the size of NTST is a good measure of the relative speed of the methods. Since both xed points for the heteroclinic orbits are saddle points and there is one positive and one negative eigenvalue in the Jacobian at (0; 0) and (1; 0), the solution v() 2 C 1 [0; 1], and the mesh points are smoothly distributed on [0; 1] for the arclength method. Example 5.2. The Lorenz equations x0 = a(y ? x); y0 = bx ? y ? xz; z 0 = xy ? cz:
At the Lorenz value (a; b; c) = (10; 13:92656; 38 ), there exists a homoclinic orbit with L = 64:89938. Starting with the Hopf bifurcation point at b = 24:73684, we use AUTO to trace out the periodic solution branch for (a; c) = (10; 83 ) until very close to the homoclinic orbit. Then we use several methods to continue the homoclinic orbit branch. The tolerance for the Newton iteration is set to 10?8. The primary and secondary continuation parameters are c and b. Their relation is shown in Fig. 5.2(a), and the solution in phase space is shown in Fig. 5.2(b). All homoclinic orbits on this branch have the stationary point (0; 0; 0). The eigenvalues of the Jacobian at the xed points are i h p ?c; ? 21 (a + 1) (a + 1)2 + 4a(b ? 1) : In Fig. 5.2(c) we plot the ratio between two negative eigenvalues of the Jacobian at (0; 0; 0). For the arclength parameterization method, using only 25 mesh intervals (NTST =25), we successfully continue the homoclinic orbit branch until c = 14:376. The nal parameter value for b is 95040:54 and the total arclength L = 592401:13. When we use NTST=50, we have no trouble in continuing the homoclinic branch to c = 14:455. Next, we compute this branch of homoclinic orbits with the projection method and the periodic method. We start both computations from a periodic solution with the xed period T = 1000 and use c and b as the primary and the secondary continuation parameters. Using 25 mesh subintervals and a tolerance of 10?8, we are unable to
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
17
1e-3 Arclength Method (NTST=10) Eigenvector Method (NTST=15)
Absolute error for c
1e-4 1e-5 1e-6 1e-7 1e-8 1e-9
0
0.2
0.4
0.6
0.8
1
0.6
0.8
1
0.6
0.8
1
a
(a) u1 1 0.8 0.6 0.4 0.2 0 0
0.2
0.4
σ
(b) u1 1 0.8 0.6 0.4 0.2 0 0
0.2
0.4
t
(c) Fig. 5.1. The Nagumo equations: (a) absolute error for c; (b) distribution of the mesh points for the arclength method; (c) distribution of the mesh points for the eigenvector method.
18
L. LIU, G. MOORE, AND R. D. RUSSELL b 1e+6
1e+5
1e+4
1e+3
1e+2
1e+1
2
4
6
8
10
12
14
16
c
(a) z 300 250 200 150 100 8/3 14/3 20/3 26/3 32/3
50 0
0
10
20
30
40
50
60
70
10
12
14
16
x
(b) λ2/λ1 20 18 16 14 12 10 8 6 4 2
2
4
6
8
c
(c) Fig. 5.2. Homoclinic orbits in the Lorenz equations: (a) solution branch in parameter c ? b space; (b) some homoclinic orbits in the x ? z plane on the branch; (c) ratio between two negative eigenvalues of the Jacobian at (0; 0; 0).
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
19
Table 5.2
Absolute error for b at c = 8:3
NTST 10 20 40 80 160
m=2 9.69e-02 2.65e-02 1.92e-03 1.21e-04 7.42e-06
m=3 5.90e-03 2.48e-04 1.53e-05 8.96e-07 5.58e-08
m=4 2.61e-05 9.45e-07 1.53e-08 2.30e-10
m=5 m=6 3.49e-06 1.88e-07 7.29e-08 5.60e-10 1.51e-09 3.00e-11
Table 5.3
Convergence rate for b at c = 8:3
NTST m = 2 10!20 1.87 20!40 3.78 40!80 3.99 80!160 4.03
m=3 m=4 m=5 m=6 4.57 4.79 5.58 8.39 4.02 5.94 5.59 4.09 6.06 5.65 4.00
perform continuation for either method. Taking the larger tolerance 10?6, we can continue the solution branch with both methods. The projection method fails at c = 14:328 and the periodic method fails at c = 13:920. However, the projection method cannot locate the solution at c = 14 during this continuation. In order to obtain the solution at c = 14, we must further raise the tolerance to 10?5. Both methods are then able to obtain the required solution. In our numerical results, the \exact" value of b and L is obtained by both arclength and projection methods using NTST=200 (they agree to the digits printed). In Table 5.1 we give the relative errors in the computed values of b for c = 14 using various values of NTST. Superconvergence is clearly observed for the arclength method and the projection method. For NTST=25 or less, the arclength method is far more accurate than the projection and periodic methods. For the periodic method, the accuracy does not improve when going from NTST=50 to 100, presumably because T = 1000 is not suciently large for this method. In Figs. 5.3{5.4, the distribution of the mesh points for c = 14 is shown for all three methods with NTST=50. Note, in Figs. 5.3(a) and 5.4(a), because of the lack of a higher order derivative for v(), there are more mesh points near = 1 to maintain the higher rate of convergence. To verify the rate of convergence for our collocation method, we compute the homoclinic orbit for the Lorenz equation at (a; c) = (10; 8:3). At this value, b = 60:5547266, the total arclength L = 303:0482482, and as we observe from Fig. 5.2, the ratio between the two negative eigenvalues of the Jacobian at the xed point (0; 0; 0) (1 = ?8:3, 2 = ?30:51594) is near minimal ( = 12 3:6766). Therefore, from our analysis in section 3.3, v() is at most C 3[0; 1]. In Tables 5.2{5.5, we show the absolute error and corresponding approximate rate of convergence for the parameters b and L. The rate of convergence is calculated by log2 (N =2N ) where N is the absolute error obtained with N mesh intervals. For the odd and even number of collocation points, it is roughly O(hm+1 ) and O(hm+2 ), respectively. In all cases, the rate of convergence is higher than the degree of continuity of v(). Clearly, the mesh selection scheme (3.23) plays the key role for the higher rate of convergence (better than O(h3)) in this example. Examples 5.1 and 5.2 clearly show that the arclength parameterization method can be very ecient for computing connecting orbits. However, recall that special care is needed if the dominant eigenvalues at the xed points are not real, i.e., the orbit
20
L. LIU, G. MOORE, AND R. D. RUSSELL x 500 450 400 350 300 250 200 150 100 50 0
0
0.2
0.4
0.6
0.8
1
σ
(a) x 500
400
300
200
100
0 0.05
0.055
0.06
0.065
0.07
0.075
0.08
0.085
0.09
0.095
0.1
t
(b) x 500
400
300
200
100
0 0.995
0.996
0.997
0.998
0.999
1
t
(c) Fig. 5.3. Distribution of the mesh points in Example 5.2 for x: (a) the arclength method; (b) the projection method; (c) the periodic method.
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
21
y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
300
350
400
450
500
300
350
400
450
500
300
350
400
450
500
x
(a) y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
x
(b) y 8000 6000 4000 2000 0 -2000 -4000 -6000
0
50
100
150
200
250
x
(c) Fig. 5.4. Distribution of the mesh points in Example 5.2 in x ? y plane: (a) the arclength method; (b) the projection method; (c) the periodic method.
22
L. LIU, G. MOORE, AND R. D. RUSSELL Table 5.4
Absolute error for L at c = 8:3
NTST 10 20 40 80 160
m=2 2.16e-01 1.10e-01 7.78e-03 4.89e-04 2.94e-05
m=3 2.25e-02 8.88e-04 5.57e-05 3.24e-06 2.01e-07
m=4 4.35e-05 3.32e-06 5.51e-08 8.00e-10
m=5 m=6 1.25e-05 6.67e-07 2.62e-07 2.00e-09 5.50e-09 1.00e-10
Table 5.5
Convergence rate for L at c = 8:3
NTST m = 2 10!20 0.97 20!40 3.82 40!80 3.99 80!160 4.05
m=3 m=4 m=5 m=6 4.67 3.71 5.57 8.38 3.99 5.91 5.57 4.11 6.11 5.78 4.01
is Sil'nikov-type. To obtain the high order convergence in such a case, it would be necessary to use a more accurate approximation such as that discussed in section 3.4. Since this would require a major change in AUTO, we use a simple implementation for which the independent variable of equation (2.4) is transformed from \time" (t) to \arclength" () directly, i.e., (); ) ; 0 < < 1; (5.1) v0() = L kff ((vv(); )k where Z 1 Z t T 0 L = T ku (t)kdt; (t) = L ku0 ( )kd ; 0 0 and the eigenvector boundary conditions (2.12){(2.14) are used. Since (5.1) has no singularity in [0; 1], the system can be solved by a standard Gauss collocation method. Note that we do not use a phase condition here but rather keep the value of xed. Example 5.3. Chua's Electronic Circuit [23] x0 = a(y ? (x)); y0 = x ? y + z; z 0 = ?by where (x) = 61 (x3 ? x). For a > 1:20245, this system has a branch of Sil'nikov-type homoclinic orbits with a xed point at the origin. The stable manifold is two dimensional. We use AUTO to trace out the periodic branch starting from the Hopf bifurcation point (1; 0; ?1) at (a; b) = (3:474937; 5) and locate a periodic solution at (a; b) = (4:35164; 5) with period T = 200. We use both the periodic method and the eigenvector method to compute the homoclinic branch with xed T = 200. With NTST=50 and the Newton tolerance set to 10?8 both methods successfully continue this branch to a > 50. For the eigenvector method, we need to increase the tolerance for the parameters since small parameters ( ) are sensitive to the relative error tolerance. The estimated absolute error for b at a = 6 is less than 10?8 in both cases. When we use the equations (5.1) and (2.12){(2.14), the continuation of the homoclinic branch is successful to around a = 10. As before, the estimated absolute error for b at a = 6 is O(10?8). In this computation, NTST=50 and the error tolerances for parameters and solutions are 10?3 and 10?8. The mesh distributions for
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
23
x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2
0
0.2
0.4
0.6
0.8
1
σ
(a) x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.5
0.55
0.6
0.65
0.7
0.75
0.8
0.7
0.75
0.8
t
(b) x 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0.5
0.55
0.6
0.65
t
(c) Fig. 5.5. Distribution of the mesh points in Example 5.3 for x: (a) the arclength method (with eigenvector boundary conditions); (b) the eigenvector method; (c) the periodic method.
24
L. LIU, G. MOORE, AND R. D. RUSSELL
all three successful methods at a = 6 are shown in Fig. 5.5. The early failure of the eigenvector method with arclength is due to the qualitative change of the solution. When a better approximate solution at larger value of a is provided, the continuation can be carried out successfully. An interesting question is whether or not there is another phase condition or a di erent treatment near the endpoints which leads to more robustness of this simpli ed arclength approach for this Sil'nikov case, or is one of the more sophisticated approaches in section 3.4 necessary. REFERENCES [1] E. Anderson, et al., LAPACK Users' Guide, SIAM, Philadelphia, 1992. [2] U. Ascher, J. Christiansen, and R. D. Russell, Collocation software for boundary value ODE's, ACM Trans. Math. Software, 7 (1981), pp. 209{222. [3] U. Ascher, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Di erential Equations, Prentice-Hall, Englewood Cli s, NJ, 1988. [4] F. Bai, A. R. Champneys, Numerical detection and continuation of saddle-node homoclinic bifurcations of codimension one and two, Mathematics Preprint 94-04, University of Bath, 1994. [5] F. Bai, A. Spence, and A. M. Stuart, The numerical computation of heteroclinic connections in systems of gradient partial di erential equations, SIAM J. Appl. Math., 53 (1993), pp. 743{769. [6] W.-J. Beyn, Global bifurcations and their numerical computation, Continuation and Bifurcations: Numerical Techniques and Applications, D. Rossed, B. D. Dier and A. Spence eds., Kluwer, Dordrecht, 1990, pp. 169{181. [7] , The numerical computation of connecting orbits in dynamical systems, IMA J. Numer. Anal., 9 (1991), pp. 379{405. [8] , Numerical methods for dynamical systems, Advances in Numerical Analysis, W. Light, ed. Clarendon Press, Oxford, 1991-1992, pp. 175{236. [9] C. de Boor, Good approximation by splines with variable knots, Spline Functions and Approximation Theory, A. Meir and A. Sharma ed., ISNM Vol. 21, Birkhauser Verlag, Basel (1973), pp. 57{72. [10] V. Canale, The Computation of Paths of Homoclinic Orbits, Ph.D. Thesis, Department of Computer Science, University of Toronto, 1994. [11] A. R. Champneys and Yu. A. Kuznetsov, Numerical detection and continuation of codimension-two homoclinic bifurcations, Int. J. Bifurcation and Chaos, 4 (1994), pp. 785{ 822. [12] L. Dieci, J. Lorenz, and R. D. Russell, Numerical calculation of invariant tori, SIAM J. Sci. Stat. Comput., 26 (1991), pp. 607{647. [13] E. J. Doedel, AUTO: a program for the automatic bifurcation analysis of autonomous systems, Congr. Numer., 30 (1981), pp. 265{284. [14] E. J. Doedel and J. P. Kernevez, AUTO: software for continuation and bifurcation problems in ordinary di erential equations, Technical Report, Applied Mathematics, California Institute of Technology, Pasadena, CA, 1986. [15] E. J. Doedel and M. J. Friedman, Numerical computation of heteroclinic orbits, J. Comp. Appl. Math., 26 (1989), pp. 155{170. [16] M. J. Friedman and E. J. Doedel, Numerical computation and continuation of invariant manifolds connecting xed points, SIAM J. Numer. Anal., 26 (1991), pp. 789{808. [17] , Computational methods for global analysis of homoclinic and heteroclinic orbits: a case study, J. Dynamics Di . Eqns., 5 (1993), pp. 59{87. [18] J. Guckenheimer and P. Holmes, Nonlinear Oscillations Dynamical Systems, and Bifurcations of Vector Fields, Applied Mathematic Sciences Vol. 42, Springer, New York, NY, 1983. [19] F. de Hoog and R. Weiss, Di erence methods for boundary value problems with a singularity of rst kind, SIAM J. Numer. Anal., 13 (1976), pp. 775{813. [20] , Collocation methods for singular boundary value problems, SIAM J. Numer. Anal., 15 (1978), pp. 198{217. [21] , An approximation theory for boundary value problems on in nite intervals, Computing, 24 (1980), pp. 227{239. [22] H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems, Proc. Applications of Bifurcation Theory, P. H. Rabinowitz, ed., Academic Press, New York,
COMPUTATION OF HOMOCLINIC AND HETEROCLINIC ORBITS
25
1977, pp. 359{384. [23] A. I. Khibnik, D. Roose and L. O. Chua, On periodic orbits and homoclinic bifurcations in Chua's circuit with smooth nonlinearity, Int. J. Bifurcation and Chaos, 3 (1993), pp. 363{ 384. [24] M. Lentini and H. B. Keller, Boundary value problems over semi-in nite intervals and their numerical solution, SIAM J. Numer. Anal., 17 (1980), pp. 577{604. [25] L. Liu, Numerical methods and software for solving boundary value ODEs for bifurcation analysis, Ph.D. Thesis, Department of Mathematics & Statistics, Simon Fraser University, 1993. [26] L. Liu and R. D. Russell, Linear system solvers for boundary value ODEs, J. Comp. App. Math., 45 (1993), pp. 1{15. [27] Y. Liu, L. Liu, and T. Tang, The numerical computation of connecting orbits in dynamical systems: a rational spectral approach, J. Comp. Phys., 111 (1994), pp. 373{380. [28] G. Moore, Computation and parametrisation of periodic and connecting orbits, IMA J. Numer. Anal., 15 (1995), pp. 245{263. [29] G. Moore and E. Hubert, Algorithms for computing stable manifolds of stationary solutions, submitted to IMA J. Numer. Anal. [30] S. Schecter, Numerical computation of saddle-node homoclinic bifurcation points, SIAM J. Numer. Anal., 30 (1993), pp. 1155{1178. [31] , Rate of convergence of numerical approximations to homoclinic bifurcation points, IMA J. Numer. Anal., 15 (1995), pp. 23{60. [32] L. P. Sil'nikov, A case of the existence of a denumerable set of periodic motions, Sov. Math. Dokl., 6 (1965), pp 163{166. [33] E. Weinmuller, A di erence method for a singular boundary value problems of second order, Math. Comp., 42 (1984), pp. 441{464. [34] , On the boundary value problem for systems of ordinary second order di erential equations with a singularity of the rst kind, SIAM J. Math. Anal., 15 (1984), pp. 287{307. [35] , Collocation for singular boundary value problems of second order, SIAM J. Numer. Anal., 23 (1986), pp. 1062{1095. [36] , Stability of singular boundary value problems and their discretization by nite di erences, SIAM J. Numer. Anal., 26 (1989), pp. 180{213.
