Traveling Wave Solutions in Coupled Chua's Circuits, Part I: Periodic Solutions *

We study a singularly perturbed system of partial differential equations that models a one-dimensional array of coupled Chua's circuits. The PDE system is a natural generalization to the FitzHugh-Nagumo equation. In part I of the paper, we show that similar to the FitzHugh-Nagumo equation , the system has periodic traveling wave solutions formed alternatively by fast and slow flows. First, asymptotic method is used on the singular limit of the fast/slow systems to construct a formal periodic solution. Then, dy-namical systems method is used to obtain an exact solution near the formal periodic soluion. In part II, we show that the system can have more complicated periodic and chaotic traveling wave solutions that do not exist in the FitzHugh-Nagumos equation.


Introduction and Preliminary
Coupled oscillators are important models for many engineering and biological systems such as neural networks [3] and cardiovascular systems [15,16].
Among all the networks, electronic networks are easy to simulate analogously. The equations can be written down precisely and the parameters of the systems can be measured easily. In this paper we consider an array of infinitely many Chua's circuits as in [3,23]. Each identical unit of the circuits consists of a parallel oscillator L and C 2 connected to a semiconductor N R through R and C 1 . The circuits are interconnected by resistors R 1 . See Figure 1.
Using k as the index for the kth circuit, the system of equations can be written as:

Figure 1. A one dimensional coupled Chua's Circuits
where G is the conductance of Chua's diode. The funcion G has negative slope in the middle of its domain, providing energy to the network to compensate the loss of energy due to the resistances of R and R 1 . If we combine V k C1 + RG(V k C1 ) into h(u), the above system can be transformed into the dimensionless form: The function h considered in [23] was a piecewise linear function. In this paper we assume that h is an N -shaped function with odd symmetry. That is: (1) h ∈ C ∞ with odd symmetry. It has three roots h(−c) = h(0) = h(c) for some c > 0; (2) h has a local minimum at u = u m , a local maximum at u = −u m , where 0 < u m < c; (3) h(u) is strictly monotone decreasing for |u| < u m and strictly monotone increasing for |u| > u m . h(u) → ±∞ as u → ±∞.
A typical example is h(u) = mu(u + c)(u − c), c > 0, m > 0. For numerical works in this paper, our test function will be H(u) = (1/30)u(u 2 − 9). Let = 1/α, ∆x = √ , u k (t) = u(t, k∆x). If ∆x is small, (1.1) can be approximated by the following PDE: u t = (y − h(u)) + 2 Du xx , y t = u − y + z, z t = −βy. (1.2) Our system is one of the simplest generalizations of the diffusive FitzHugh-Nagumo equation, which is a second order bistable PDE coupled with a linear first order ODE. In FitzHugh-Nagumo's system the reduced 1D slow equation has an equilibrium with only one real eigenvalue [13,14]. However, in our systemthe reduced 2D slow system of ODEs has an equilibrium with two complex eigenvalues if β > 0; or one positive, one negative eigevalues if β < 0. The case β < 0 was studied by Acosta in his 1998 thesis [1], while in this paper we only consider the physically meaningful case β > 0.
Denote the traveling wave solution of (1.2) by w = w(x − st) := w(ξ), where w = (u, v, y, z). We switch the independent variable ξ to t so the traveling wave becomes w(t) andẇ = dw/dt. The system obtained is called the slow system where t is the slow time scale:   For y < −y m or y > y m , only the solution u − (y) or u + (y) exists. Observe that Dh(u ± (y)) > 0 and Dh(u 0 (y)) < 0.
The slow manifold S consists of equilibrium points of (1.6). We are interested in the most left and right submanifolds S − and S + . Given any 0 < y 0 < y m , we show that for y = −y 0 , there exists a unique wave speed s = s 0 > 0 such that (1.6) has an increasing heteroclinic orbit w 0 1 connecting S − to S + . For the same wave speed s 0 but y = y 0 , (1.6) has a decreasing heteroclinic orbits w 0 3 connecting S + to S − . To form a closed loop, the layer solutions (w 0 1 , w 0 3 ) must be connected by two orbits of (1.5) on the slow manifolds S − and S + . The orbit w 0 2 on S + connects u = u + (−y 0 ) to u = u + (y 0 ) while the orbit w 0 4 on S − connects u = u − (y 0 ) to u = u − (−y 0 ). For > 0 and small, under some generic conditions, the existence of an exact solution near the singular limit approximation has been proved in singular perturbation theory. In some papers, the proof is based on the graph transformation theory called the Exchange Lemma, similar to the Inclination Lemma or Lambda Lemma in dynamical systems [7-9, 11, 25]. Such geometric method works better if the boundary conditions at the two ends are isolated, such as the Dirichlet or Norman boundary conditions, but should also work with the periodic boundary conditions by an iteration scheme as in [4].
In this paper, we use an analytical method to show the existence of the exact solution. We look for error correction terms to the singular orbits to form the exact solution for > 0. Rather than tracking the nonlinar manifolds, the analytic method tracks the linear manifolds under the flow, called the invariant stable and unstable subspaces of exponential dichotomies [5,22]. We also use an iteration method, but the iteration is on a system of linear equations. To solve the system of linear equations, we use a combination of methods from the classic singular perturbation theory and the method of heteroclinic bifurcations as in [10,19,20]. If the linear system for the error corrections has been solved, the nonlinear system can be solved by the contraction mapping principle.
We present notations used in this paper and an outline of our approach in §2. In §3, we construct the 0th order approximation for the singular periodic traveling wave. The singular orbits in internal layers and outer layers are constructed in §3.1 and §3.2 respectively. In §3.3, we combine the approximations in regular and singular layers to form the 0th approximation with small residual and jump errors. By checking some generic conditions as in [19], we conclude in §4 that there is an exact periodic traveling wave for > 0 and small, near the singular one for = 0. The outline of the proof is given in §5. After deriving the linear variational systems for the correction terms in singular and regular layer, in §5.1 we derive the jump conditions between the regular and singular layers. In §5.2 we solve the boundary value problems that determine the correction terms (u i , v i ), i = 1, 2, 3, 4. In §5.3 we obtain the solutions (y i , z i ), i = 2, 4 on regular layers. Finally, in §5.4 we show how to use the correction terms to show the existence of the exact periodic traveling wave for > 0 and small. This research constitutes part of the doctoral dissertation of Ming Jiang who is grateful to Professor X.-B. Lin and faculty of the Mathematics Department of NCSU for their advise and support.

Notation of this paper and an outline of the approach
Consider a linear ODE in R 2 where I is a bounded or unbounded interval and A(τ ) is continuous and uniformly bounded on I. Let T (τ, σ) be the solution matrix for the system. Definition 2.1. We say (2.1) has an exponential dichotomy in I if there exist projections P s (τ ) + P u (τ ) = I, constant K ≥ 1 and exponent α > 0 such that For example, if A(τ ) = A 0 is a constant matrix with eigenvalues λ 1 < 0 < λ 2 , then x = A 0 x has an exponential dichotomy, and the projections are spectral projections to stable and unstable eigenspaces. If A(τ ) depends on τ , then the existence of exponential dichotomies can be proved by the perturbation theory (or the roughness) of exponential dichotomies, [5,22]. Two such cases is given in the next lemma.
(2) Assume that for each fixed τ ∈ I, the matrix A(τ ) has eigenvalues λ 1 (τ ) < 0 < λ 2 (τ ) where |λ i (τ )| ≥ c > 0 uniformly for t ∈ I. If the system is slowly varying in the sense that |A (τ )| << c, then the time dependent system (2.1) has an exponential dichotomy in I. by gluing the local functions together, where is called the concatenation operator defined as follows. Assume x(0) = x 0 (0) to fix the phase. If a 0 ≤ t ≤ b 0 , then x(t) = x 0 (t); if t > b 0 , then there exists some j ≥ 1 such that and if t < a 0 then there exists some j < 0 such that If j = −1 then the summation j+1 i=−1 is ignored. If j < −1, the indices in the summations must count down from i = −1 to j or j + 1 by step size −1. We do not require b i = a i+1 so the coordinate on I i is local. This gives us the flexibility of changing the end points a i , b i of I i locally without changing the definitions of x i on other intervals.
, we can still define the concatenated function x(t) which has a jump of size J i,i+1 between two consecutive intervals.
The exact periodic solution of (1.3), = 0, will be denoted by w ex . Using local time coordinates, we look for w ex i (t), > 0 on intervals For even i, w ex i satisfies (1.3). For odd i, w ex i satisfies (1.4) and we will change the variable t → τ so the function is w ex i (τ ) defined on [a i / .b i / ]. The exact solution will be obtained as the limit of a sequence of approximations denoted by w j , j = 0, . . . , ∞, so that in some suitable topology (not in L ∞ (R), since the period of each iteration will be different): The first term w 0 in the sequence is the singular periodic orbit, obtained by solutions of (1.5), (1.6) where = 0. The approximation w j+1 is constructed from w j by an iteration process to reduce the error of approximation by a factor 0 < c < 1.
Similar local time coordinates also apply to the jth approximation w j which is the union of fast and slow layers defined on the intervals I j i , i = 1, 2, 3, 4. For i = 2, 4, we use the slow time t for the solutions so that w j i = w j i (t). For i = 1, 3, we change t to the fast time τ so that w j i = w j i (τ ) where τ = t/ . To calculate the approximations by iteration, let w cor,j i be the corrections so that w j+1 To construct w cor,j i from the the approximation w j i , we will use a linear variational system around the 0th order approximation, not the linearized system around w j i . The notations on the super and sub scripts of w also apply to its components w = (u, v, y, z) throughout the paper. Finally by the periodicity, w j i = w j i+4 , so we can limit our attention to i = 1, 2, 3, 4.

Constructing the 0th approximation for periodic orbits
In this section, we construct formal periodic solutions for the singular system = 0. Then using the formal solutions we formulate an approximation of the periodic traveling wave for > 0 and small. Let U = (u, v), Y = (y, z) and w = (U, Y ). The singular periodic solution w 0 = (U 0 , Y 0 ) is the concatenation of 4 pieces of orbits w 0 i , 1 ≤ i ≤ 4. The fast orbits w 0 i , i = 1, 3 satisfy (1.5) and are defined on singular layers where the solutions are written as w 0 . The slow orbits w 0 i , i = 2, 4 satisfy (1.6) and are defined on regular layers where the solutions are written as We extend the definition of orbits to i ∈ Z by w 0 i+4 = w 0 i for convenience. The following asymptotic matching conditions must be satisfied: (ii) If (Y 0 2 , Y 0 4 ) are a pair of slow orbits on (S + , S − ) and form a closed loop when projected to the (y, z) plane. Y 0 2 connects y = y 0 1 to y = y 0 3 and Y 0 4 connects y = y 0 3 back to y = y 0 1 .
Although (U i , Y i ) are 2D vector-valued functions and S ± are 3D manifolds in the 4D space. To avoid complicated notation, the projections of S ± to (u, v) components or (y, z) components will be denoted by the same notation S ± if no confusion should arise.

Singular periodic solutions in the internal layers
With (y, s) as parameters, we look for U 0 i , i = 1, 3 that satisfies the first two equations of (1.6): which can be written as the second order ODE: The equilibrium points of (3.1) satisfy (h(u) = y, v = 0) and form a smooth manifold. For convenience, we shall denote the manifold S, S ± , S 0 and its projection to the (u, y) cooredinates by the same notation. By this convention, we say the equilibrium point of (3.1) or the graph y = h(u) is in S and belong to the branch S − , S + or S 0 if u ≤ −u m , u ≥ u m or |u| < u m . Using the characteristic equation Dr 2 + sr − h (u) = 0, we conclude that the equilibrium points on S − and S + are hyperbolic saddle points, while the equilibrium points on S 0 are stable (or unstable) if s > 0 (or s < 0).
To further study the relation between y 0 and the wave speed s 0 , we introduce a linear version of Melnikov's method.
∼ Melnikov's method for bounded solutions on the infinite domain. Let (u(τ ), v(τ )) T be a heteroclinic orbit of (3.1) connecting two saddle points from S − → S + or S + → S − , and let Since (u(τ ), v(τ )) T → S ± as τ → −∞, or ∞, due to the fact that A(±∞) has one positive and one negative eigenvalues. from Lemma 2.1, U = A(τ )U has exponential dichotomy on R − or R + respectively. But the linear system has no exponential dichotomy on R due to the non-transversal intersection of the ranges RP u (0−) and RP s (0+) at τ = 0. It is clear that, up to a constant factor, the derivative of the heteroclinic solution (u , v ) T is the unique bounded solution to the 2x2 system U = A(τ )U . Therefore, up to a constant factor, Ψ = e sτ /D (−v , u ) T is the unique bounded solution to the adjoint equation Ψ + A(τ ) T Ψ = 0.
Let C k (−∞, ∞) be the space of C k and uniformly bounded functions on R. We consider the following nonhomogeneous equation For any H ∈ C(−∞, ∞), the linear system (3.5) has a unique solu- If (3.6) does not hold then let the left hand side be g. There exists a unique piecewise Moreover, we have the following estimates for (g, U ): ∼ Existence of internal layers and the relation between parameters.
We go back to consider (3.1). For the cubic nonlinearity h stated before, it is known that given y = y 0 ∈ (−y m , y m ), there exists a unique s 0 such that (3.5) has a heteroclinic orbit connecting S − → S + or S + → S − . In fact, for the "typical system" u + su + u(u − a)(1 − u) = 0, 0 < a < 1, the exact solution corresponding to the unique wave speed has been obtained by Huxley [24]. In the following lemma we show the relation between s 0 and y 0 for the existence of an increasing heteroclinic orbit connecting saddle points u − (y 0 ) to u + (y 0 ). Similar results on the the decreasing heteroclinic orbit, from u + (y 0 ) to u − (y 0 ), can be obtained by symmetry. Before stating the lemma, consider the limiting case y 0 = −y m where u − (y 0 ) is still a saddle point on S − , while u 0 (y 0 ) = u + (y 0 ) becomes the local minimum point (turning point) of S where S 0 meets S + . It is known that there exists a minimum wave speed s m such that for any s ≥ s m , the heteroclinic orbit connecting the saddle point (u − (y 0 ), y 0 ) on S − to the turning point (u, y) = (u m , −y m ) on S exists, see [2].
there exists a unique s = s * (y 0 ) such that system (3.1) has a unique increasing heteroclinic solution (u 0 (τ ), v 0 (τ )) connecting (u − (y 0 ), 0) to (u + (y 0 ), 0). Moreover, ds * (y0) Proof. Part (1): Assume that (u 0 , v 0 ) is a heteroclinic solution of (3.1) corresponding to parameters (s 0 , y 0 ). Let (u 0 +∆u, v 0 +∆v) be a nearby solution with the parameters (s 0 + ∆s, y 0 + ∆y). The linear variational system for ∆U = (∆u, ∆v) T is The linear homogeneous part of the above has a unique bounded solution To have a C 1 solution near (u 0 , v 0 ) we must have g(∆s, ∆y) = 0. Observe that Therefore equation g(∆s, ∆y) = 0 can be solved for ∆s = s * (∆y) with Part (2): When y 0 = −y m , Billingham and Needham [2] showed that there is a minimum wave speed s m such that a heteroclinic connection from the saddle point (u − (y 0 ), 0) to the turning point (u 0 (y 0 ), 0) exists iff s ≥ s m . When s = s m , the connection is from the unstable manifold W u (u − (y 0 )) of the saddle point to the stable manifold W s (u 0 (y 0 )) of the turning point, while for s > s m the connection is tangent to the weakly stable part of the center manifold of the turning point.
Increasing y 0 from −y m , the function u = h −1 (y 0 ) at (u m , −y m ) splits into two branches -one is u 0 (y 0 ) < u m , and the other one is u + (y 0 ) > u m . If we follow the branch formed by (u + (y 0 ), 0) on S + , then the connection of W u (u − (y 0 )) to W s (u + (y 0 )) persists and the wave speed s * (y 0 ) is a continuous decreasing function of y 0 . This shows that s m = sup{s * (y 0 ) : |y 0 | < y m }. As y 0 → ∓y m , the wave speed s * (y 0 ) monotonically approaches ±s m .
To make the picture more complete, if we follow the other branch where the turning point becomes a stable node (u 0 (y 0 ), 0) on S 0 . The saddle to node connection exists if the wave speed s > s * (y 0 ), which is a well-known fact in ODEs.

Singular periodic solution in the slow layers.
Recall that (s 0 , −y 0 ) are the parameters for U 0 1 and (s 0 , y 0 ) are the parameters for U 0 3 . We look for w 0 2 on S + and w 0 4 on S − so thatY 0 2 connects y = −y 0 to y = y 0 , and Y 0 4 connects y = y 0 back to y = −y 0 . Observe that the slow system and its limit, (1.3) and (1.5), obey the odd symmetry about the origin, so we will assume that the solutions w 0 2 and w 0 4 are symmetric about the origin. Plug u = h −1 ± (y) into the last two equations of (1.5), we have: Several orbits of (3.9) are plotted in Figure 3 with the z-axis pointing to the right and the y-axis pointing upward. The left half of each smooth arc represents Y 0 2 on S + plotted with the nonlinear function h −1 + (y). The right half of each smooth arc represents Y 0 4 on S − plotted with the nonlinear function h −1 − (y). The orbit Y 0 2 (or Y 0 4 ) rotates around the equilibrium points (z, y) = (−c, 0) (or (z, y) = (c, 0)) in clock wise direction. They are symmetric to each other and form a continuous closed loop, but the derivatives are discontinuous at the junction points.
On S + , the z-nullcline is defined by y = 0, and y-nullcline is defined by z = y − h −1 + (y). The quilibrium point (z, y) = (−c, 0) separates the z-nulcline into the left and right half lines, and separates the y nullcline into the upper lower half lines. The smallest loop has junction points on y 0 = 0. To save space, the z-axis is horizontal and the y axis is vertical.
Remark 3.1. We comment that to avoid computing h −1 ± (y) numerically, we first compute orbits on S using the (z, u) coordinates: then map the orbits to the (y, z) plane. The y-nullclines, z = y − h −1 ± (y) on S ± , are also obtained by computing z = h(u) − u first then mapped to the (z, y) plane.
In the rest of this subsection, we study the existence of Y 0 2 , Y 0 4 analytically. The following hypothesis ensures that the divergence of the vector field of (3.9) is negative in the region of interest: (H1): For |y| ≤ y m , we assume that (d/dy)h −1 ± (y) > 1. For the cubic polynomial h(u) = mu(u 2 − c 2 ), condition (H1) is satisfied if m < 1/(3c 2 ). In particular, it is satisfied by the test function H(u). It is easy to see Y A and Y B obey the odd symmetry. In Fig. 4, (A, Y A ) and (B, Y B ) are plotted with 0 >ȳ ≥ −y m . The y coordinate y A (t) is locally decreasing (or increasing) to the right (or left) of the y-nullcline. So Y A (t) will hit the znullcline in backward time at some point C and hit the z-nullcline, y-nullcline and z-nullcline consecutively in forward time at (K, A , D). Similarly Y B (t) will hit the z-nullcline and y-nullcline consecutively at (E, B, L, B , F ).
Using the divergence theorem and (H1), we can show that 0 < A y < |A y |, and 0 > B y > −|B y |. Therefore, (Y A (t), Y B (t)) stay within the region |y| ≤ y m when they rotate one cycle around (−c, 0) and (c, 0) as in Figure 4. By the same method we can show −c < D z < C z and E z < F z < c.  Moreover, if between the intersecting points, Proof. In Figure 4, y-nullclines y − h −1 ± (y) − z = 0 are plotted as part of the mirror image of an N shaped curve. The left branch of the curve is the y-nullcline on S + and the right branch of the curve is the y-nullcline on S − .
Recall that if (z(t), y(t)) is a solution of (3.9) on S + where h −1 = h −1 + , then (−z(t), −y(t)) is a solution of (3.9) on S − where h −1 = h −1 − . From the symmetry of the solutions Y A and Y B and the property (P1), E z = −C z < D z in Figure 4. From (P2), E z = −C z ≥ −c > A z . From (P3), D z ≤ −A z = B z . We then have the following ordering on the z-axis: Consider the y coordinates of the curves. We see that E is under the curve A D, and D is under the curve EB. Thus A D and EB must intersect at some point G, of which the z coordinate is between E and D.
At the point G we have dz A (t)/dt = dz B (t)/dt > 0. By A z < G z < D z we have dy A (t)/dt < 0 , and by E z < G z < B z we have dy B (t)/dt > 0. The intersection of Y A and Y B at G is transversal.
Although (Y A , Y B ) are solutions to (3.9), the phase conditions are not satisfied. Let d > 0 be such that Y A (d) = K, Y B (d) = L. The construction of the singular slow orbits are complete if we let Y 0 The construction of (Y 0 2 , Y 0 4 ) can be viewed as the fixed points of Poincare mappings defined below. The two lines y = ±y 0 are transverse to the flow of (3.9). So we can define the Poincare mapping Π 2 that maps points on the line Σ(−y 0 ) := {y = −y 0 } to the line Σ(y 0 ) := {y = y 0 }, and the Poincare mapping Π 4 : Σ(y 0 ) → Σ(−y 0 ). The two mappings are symmetric about the origin. We may say that Y 0 2 (a 2 ) is a fix point of the the mapping Π 4 • Π 2 , and Y 0 4 (a 4 ) is a fixed point of the mapping Π 4 • Π 2 .
In the future we need the following assumption.
(H2): The Poincare mappings Π 2 and Π 4 are exponentially stable in a neighborhood of Y 0 2 (a 2 ) and Y 0 4 (a 4 ). Assumption (H2) is based on the numerical simulation of our system. It seems to be related to (H1) but we do have a rigorous proof of that. Definition 3.2. We denote −0.5 by for the convenience. We sometimes say w 0 i , i = 1, 2, 3, 4, are defined on [α i , β i ] with the understanding that α i = a i , β i = b i , i = 2, 4, are independent of while α i = − , β i = , i = 1, 3, are unbounded as → 0.
Here w 0 i (t), i = 2, 4 and w 0 i (τ ), i = 1, 3 are formal solutions to the singular limit systems. If = 0, (u 0 , v 0 , y 0 , z 0 ) does not satisfy (1.3), or (1.4). The formal approximation has residual errors in singular and regular layers. It is straightforward to check that in the slow layers, the residual errors R 0 i (t), i = 2, 4 satisfy: (3.10) In the fast layers, the residual errors R 0 i (τ ), i = 1, 3 satisfy: (3.11) When = 0, due to the matching condition, there is no jump error between singular and regular layers. When = 0, the approximation has jump errors at the junction points.
We use J i,i+1 to denote the jump errors of an approximations for > 0. Since the internal layers approach the limits exponentially as τ → ±∞, for the 0th approximation, it is easy to see that there exists γ > 0 such that

Existence of Periodic Solutions for small > 0.
Let us rewrite (1.3) as General conditions for the existence of exact solutions of boundary value problems near the formal approximations have been obtained in [19] where the boundary conditions include isolated boundary conditions and the periodic boundary conditions. For the periodic solutions, these conditions are (based on [19]): (C1): We assume the normal hyperbolicity at the approximation on S ± , (C2): The breaking of the heteroclinic solutions due to the change of Y is nontrivial. It is represented by the following Melnikov integral: After a rescaling of Ψ i , ∆ i = (1, 0) for our system. It reflects the fact that the change of the gap between U 0 i (0+) and U 0 i (0−) depends on the parameter y, not z.
It is straightforward to check that Conditions (C1)-(C3) are satisfied. Condition (C4) is from the hypothesis (H2) which based on the numerical results. Therefore we obtain the existence of a real periodic solution near the singular periodic solution from Theorem 7.1 in [19].
However, the paper [19] is not dedicated to periodic solutions. The focus of that paper was on singular perturbation problems with isolated boundary conditions and asymptotic expansion of the solutions to any order of , which are irrelevant to this paper. Thus we will give an outline of the proof for the existence of a real periodic solution even the result is not new.

Outline of the proofs
We shall construct a sequence of approximations {w j } ∞ j=0 , using the singular periodic orbit as the initial approximation w 0 . We present an iteration method to obtain the j + 1th approximation from the jth approximation in this section. Each approximation should reduce the errors by a factor 0 < c < 1, so the sequence {w j } ∞ j=0 converges to the exact periodic solution. For > 0, the process of obtaining w j+1 from w j can be viewed as an operator that takes any approximation w ap to a better approximationw ap . Since our goal is a better approximation, we have the freedom of dropping higher order terms which can make the iteration much easier to carry out. Let w ex be the exact periodic solution near an arbitrarily approximation w ap . The correction term w = w ex − w ap in I i will be denoted by w cor i . We will drop the superscript cor if no confusion should arise. This goes to all the components (u, v, y, z) of w.
We use fast time τ for corrections in fast layers and use slow time t for corrections in slow layers, so: where δa l , δb l are correction terms to boundaries of the domain I i , i = 2, 4. Let the (U, Y ) components of the residual errors of the approximation w ap in I i be denoted by (p i (t), q i (t)), i = 2, 4 and (p i (τ ), q i (τ )), i = 1, 3. For the initial approximatoin w 0 which are of O( ) as from (3.10) and (3.11). The approximation and the residual error satisfy the following equations: Thus the exact correction terms satisfy: By higher order term we mean a function of the form o( + |w|) and its Lipschitz number with respect to the unknown function w i is o(1). Let (p i , q i ) = −(p i , q i ). The linear variational system for the approximate correction terms can be written as: The coefficients are calculated by linearizing around (U 0 i , Y 0 i ). We use the slow time t for i = 2, 4 and the fast time τ for i = 1, 3 so that Similar rules apply to F Y , G U and G Y . We will drop the h.o.t. in (5.1) since we only look for approximate correction terms. Then we use the smallness of to simplify the system to a triangular form which can easily be solved.
(i) For i = 1, 3: Since |Y i | L ∞ ≤ C|Y i | L 1 , and the L 1 (− , ) norm of the two terms, can be considered as higher order terms. We cannot drop q i even it is small. The simplified fast system is: (ii) For i = 2, 4: From (5.1), the vectors near the slow manifold S decompose into two components. The first is tangent to the slow manifold, and the second is parallel to the U -axis, each of them is invariant under the singular homogeneous flow of (5.1): We describe the decomposition in the following figure.  0) is on the U -plane where system U = F U U has an exponential dichotomy for τ ≥ 0 or τ ≤ 0.

Jump conditions for the intervals with varying end points
Recall that I i = [α i , β i ] for all i. For i = 1, 3, we use the fast time τ , then α i = − , β i = . For i = 2, 4, we use the slow time t, then α i = a i , β i = b i . Define jump between the adjacent intervals of the approximations by: To correct the jump errors, we allow the boundaries of I 2 , I 4 to change: The total change of boundary values on I i , i = 2, 4, partially comes from the correction terms (Y i , U i ) and partially comes from the change of boundary points. Thus, for i = 2, 4, if c = a i or b i and δc is the variation of that point, then For i = 1, 3, the domains are fixed so At the two end points of I 3 = (α 3 , β 3 ), Let J i,i+1 = −J i,i+1 be the jump for the corrrectioni terms and let J i,i+1 Y be its U and Y components. Subtracting the two lines and dropping the higher order terms, we have Here we replacedẎ ap i , i = 2.4 byẎ 0 i with errors that are of higher order terms. Similar results hold for the end points of I 1 . Thus the jump conditions for correction functions Y i , i = 1, 2, 3, 4 and variation of junction times δa i , δb i , i = 2, 4, are : where we have dropped the higher order terms and used the convention Similarly, the linear variational systems for the correction functions U i , i = 1, 2, 3, 4 are Define the jump for the approximation solution between two slow intervals as: Then the corrections to the jumps satisfy To summarize, the correction terms w i , i = 1, 2, 3, 4 must satisfy the following system: where A i (τ ) and A i (t) are 4x4 matrices, f i (τ ) and f i (t) are the nonhomogeneous terms of the linear system. The vector d i ⊥ẇ i ap (0) determines the phase condition of w i . Since solutions of linear the system are determined by the jump conditions, we call the system a linear non-homogeneous jump value problem.
From (5.2), dropping the h.o.t, the precise form of (5.8) is: Naturally we can solve the 2nd equation then plug into the first: (5.14) In the rest of the seciton we look for solutions of (5.13), (5.14) with jump conditions (5.10), (5.11).

(5.18)
We now consider the equations for U i , i = 1, 3 and V i , i = 2, 4.
For i = 2, 4, the system V i (t) = F U (t)V i (t) has exponential dichotomy on [a i , b i ] with the exponents being O(1/ ). This can be shown by changing the time t to τ = t/ so the system becomes V i (τ ) = F U ( τ )V i (τ ), i = 2, 4. For each fixed τ ∈ [a i / , b i / ], the system has exponential dichotomy. The coefficient is slow varying since (d/dτ )F U ( τ ) = O( ). Thus the rescaled system has exponential dichotomy in [a i / , b i / ] as from Lemma 2.1. After scaling back τ → t, the system V i (t) = F u (t)V i (t), i = 2, 4 has an exponential dichotomy and the exponents become O(1/ ).
At the junction point between I i and I i+1 , the vectors w i ap , w i+1 ap are near S, so we split the jump into (V, Y ) components. From (5.10), (5.11), the jump conditions for V i are For c = a i or b i , from the definition of V i , and Remark 5.1. To simplify the notation, we let V i = U i in internal layers i = 1, 3.
Dropping the higher order term, the jump condition on V i , i = 1, 2, 3, 4 are We shall solve (5.19) with the jump conditions (5.20) in three steps.
Step 1: LetV i , i = 1, 2, 3, 4 be the solutions of (5.19) with the boundary con- (5.19) has exponential dichotomies on I i . If the right hand side of (5.19) is denoted by f i (t), the solution exists, and are bounded by C|f i | ≤ C(|p i | + |q i |).
For i = 1, 3 system (5.19) has exponential dichotomies on the left and the right half of the interval I i . From Lemma 5.1, there exists a uniqueȳ 0 i such that if Y i (0) = (ȳ 0 i ,z 0 i ) T , then the solutionV i exist, and is bounded by (|p i | + |q i |) as well.
Step 2: The functionsṼ i = V i −V i , i = 1, 2, 3, 4 are solutions of (5.19) with p i , q i = 0. The jump conditions are modified tõ (5.21) Let Q i,i+1 be the projection to the space RP s (α i+1 ) with the null space N Q i,i+1 = RP u (β i ). Then the jump conditions are approximated by the boundary conditions: We have finally obtained the boundary conditions forṼ i , i = 1, 2, 3, 4, The system forṼ i with boundary conditions is homogeneous: The system forṼ i , i = 2, 4 has an exponential dichotomy in the domain I i . Solutioñ V i (t) = T (t, a i )φ i s + T (t, b i )φ i u that satisfies the boundary conditions uniquely exists and is bounded by |φ i s | + |φ i u |.
The system forṼ i , i = 1, 3 has exponential dichotomies on [− , 0] and [0, ] respectively. There exists a uniqueỹ 0 i such that if Y i (0) = (ỹ 0 i ,z 0 i ), then the solution of the boundary value problemṼ i exists and is bounded by |φ i s | + |φ i u |. One nice fact is that the uncontrolled boundary terms are exponentially smaller than the specified boundary terms |φ i s | + |φ i u |. Although the jump conditions (5.21) are not be satisfied by the solutions of the boundary value problemṼ i , but the errors are smaller than K i,i+1 V by a factor c < 1. Using an iteration scheme, the solutions with exact jumps as in (5.21) can be obtained.
Step 3: If we let V i =V i +Ṽ i , then equations (5.19) and jump conditions (5.20) are satisfied. The parameter Y i (0) = (y 0 i , z 0 i ) T where y 0 i =ȳ 0 i +ỹ 0 i has been determined but z 0 i is arbitrarily. In summary we have obtained corrections V i , i = 1, 2, 3, 4 to any given approximation. Now suppose the given approximation is V j i , then the next approximation But for U j+1 i , i = 2, 4, we need to know Y j+1 i , i = 2, 4, which will be obtained in the next subsection.

Poincare mappings for linear systems and construction of
(Y 2 , Y 4 ).
We have constructed (Y 1 , Y 3 ) by the first equation of (5.13). Let The values of y 0 i , y ± i have been completely determined. For the z coordinate, only the value of the gap z + i − z − i has been determined. The purpose of this subsection is to construct Y j+1 i , i = 2, 4 based on Y j i , i = 2, 4. We will use a linear version of Poincare mappings to find Y 2 , Y 4 . As always, linear Poincare mappings are tricky to define. To illustrate the intuitive idea, we first show how the next iteration Y j+1 i should look like in Fig. 6. The linear variational system for Y i iṡ ds is a solution to (5.24). By subtractingȲ i from (5.24), we only have to look forỸ i = Y i −Ȳ i that satisfies a homogeneous system. We also update the boundary conditions and jump conditions related toỸ i to include the contributions fromȲ i (t). The equation forỸ i , i = 2, 4 must satisfy the following system (we rewriteỸ i as Y i to simplify the notations):  Figure 6. The size of jumps BC and DA are given on the lines y = y 0 + y 4 −, y = y 0 + y + 2 and y = −y 0 + y − 2 , y = −y 0 + y + 4 .
where the boundary conditions in y and the jump conditions in z for i = 2, 4 are: The line {(y, z) : y =ȳ, z ∈ R} is called a y section and denoted by Σ(ȳ). Consider the three y sections: Define the linear affine map (linear Poincare map) from t = 0 to t = b i , Σ(0) → Σ(y + i ) as follows. Given (0, z i ) T ∈ Σ(0), there exists a unique δb i such that Also define the backwards linear affine map from t = 0 to t = a i < 0, Σ(0) → Σ(y − i ) as S(a i , 0)(0, z i ) T + δa i ·Ẏ 0 i (a i ) ∈ Σ(y − i ), i = 2, 4.
In the previous subsection, the solutions Y i , i = 1, 3 has been constructed by integrals with one degree of freedom in the choice of z i (0), i = 1, 3. It is now fully determined too. If Y ap i , i = 2, 4 is the jth approximation Y j i in the interval I i , then the (j + 1)th approximation Y j+1 i (t) = Y ap i + Y i (t), i = 1, 2, 3, 4 has been constructed. Moreover, using Y i , i = 2, 4, we obtain U j+1 i as Before concluding that w j → w ex , as j → ∞, in some suitable topology, we point out that the sequence of functions {w j+1 (t), t ∈ R} ∞ j=0 do not satisfy the Cauchy criterion in the super norm of t since the period of each approximation is different. To complete the proof we can show that orbitally the sequence is convergent, see [19]. Alternatively, using the local coordinates we have a natural way to show in what sense the approximations approach the exact solution.
First, based on our construction, if we compare w j+1 i to w j i , the following residual errors and jump errors decay by a factor 0 < c < 1 The correction terms w i and change of boundary points δa i , δb i at each iteration are controlled by the error terms and jump terms as in (5.8) to (5.11). So sup{(|w i | + |δa i | + |δb i |) : i ∈ Z} also decays by a factor 0 < c < 1 when j increases to j + 1.
Then the existence of an exact periodic solution w ex which is the concatenation of w ex i , i ∈ Z can be proved in 5 easy steps. Details will not be presented.
(1) The domain [a j i , b j i ] approaches the limit [a ex i , b ex i ] as j → ∞.
(2) The function w ex i is well defined on [a ex i , b ex i ]. (3) The function w ex i satisfies the nonlinear system with no residual error. (4) The sequence w ex i , i ∈ Z has no jump between two successive intervals, i.e., w ex i (b ex i ) = w ex i+1 (a ex i+1 ). (5) The function w ex is periodic in t due to the fact w ex i = w ex i+4 .