Bifurcations and chaos in two-cell cellular neural networks with periodic inputs

(1)

World Scientiﬁc Publishing Company

BIFURCATIONS AND CHAOS IN TWO-CELL

CELLULAR NEURAL NETWORKS

WITH PERIODIC INPUTS

SONG-SUN LIN∗

Department of Applied Mathematics, National Chiao-Tung University, Hsin-Chu 30050, Taiwan

WEN-WEI LIN†

Department of Mathematics, National Tsing-Hua University, Hsin-Chu 30050, Taiwan

TING-HUI YANG∗

Department of Applied Mathematics, National Chiao-Tung University, Hsin-Chu 30050, Taiwan

ReceivedMay 30, 2003; RevisedSeptember 8, 2003

This study investigates bifurcations and chaos in two-cell Cellular Neural Networks (CNN) with periodic inputs. Without the inputs, the time periodic solutions are obtained for template

A = [r, p, s] with p > 1, r > p− 1 and −s > p − 1. The number of periodic solutions can

be proven to be no more than two in exterior regions. The input is b sin 2πt/T with period

T > 0 andamplitude b > 0. The typical trajectories Γ(b, T , A) and their ω-limit set ω(b, T , A)

vary with b, T and A are also considered. The asymptotic limit cycles Λ_∞(T , A) with period

T of Γ(b, T , A) are obtainedas b → ∞. When T₀ ≤ T₀∗ (given in (67)), Λ_∞ and −Λ∞ can be separated. The onset of chaos can be induced by crises of ω(b, T , A) and −ω(b, T, A) for suitable T and b. The ratioA(b) = |aT(b)|/|a₁(b)|, of largest amplitude a₁(b) except for T -mode andamplitude of the T -mode of the Fast Fourier Transform (FFT) of Γ(b, T , A), can be usedto compare the strength of sustainedperiodic cycle Λ₀(A) andthe inputs. WhenA(b) 1, Λ₀(A) dominates and the attractor ω(b, T , A) is either a quasi-periodic or a periodic. Moreover, the range b of the window of periodic cycles constitutes a devil’s staircase. When A(b) ∼ 1, ﬁnitely many chaotic regions and window regions exist and interweave with each other. In each window, the basic periodic cycle can be identiﬁed. A sequence of period-doubling is observed to the left of the basic periodic cycle and a quasi-periodic region is observed to the right of it. For large b, the input dominates, ω(b, T , A) becomes simpler, from quasi-periodic to periodic as b increases.

Keywords: Cellular neural networks; CNN; chaos; crises; fractal; Lady’s shoe; Lyapunov

exponent.

∗_{Work partially supported by the NSC under Grant No. 89-2115-M-008-029, the Lee and MTI Center for Networking Research}

and the National Center for Theoretical Sciences Mathematics Division, R.O.C.

†_{Work partially supported by the NSC under Grant No. 91-2115-M-007-006 and the National Center for Theoretical Sciences}

Mathematics Division, R.O.C.

3179

Int. J. Bifurcation Chaos 2004.14:3179-3204. Downloaded from www.worldscientific.com

(2)

1. Introduction

Following the introduction by Chua and Yang [1988a, 1988b], Cellular Neural Networks (CNN) have been extensively studied, see [Chua, 1998; Manganaro et al., 1999; Chua & Roska, 2002] and references therein. Two of their applications are in image processing andpattern recognition. An important class relatedto applications are steady-state solutions including mosaic solutions and defec-tion soludefec-tions [Chua, 1998; Manganaro et al., 1999; Hsu et al., 2000; Juang & Lin, 1997, 2000]. The complexity of steady-state solutions have recently been extensively studied [Ban et al., 2001a; Ban et al., 2002; Ban et al., 2001b; Hsu et al., 2000; Hsu & Lin, 2001; Hsu & Yang, 2002; Juang & Lin, 2000; Juang & Lin, 1997; Lin & Shih, 1999; Lin & Yang, 2000; Lin & Yang, 2002]. Furthermore, with-out the input terms, the theory of complete stabil-ity for CNN with symmetric feedback template have been proven in [Lin & Shih, 1999; Shih, 2001; Wu & Chua, 1997]. However, when the feedback template is antisymmetric, the time dependent periodic solu-tions have been obtainedby Thiran [1997].

Zou andNossek [1991] discovereda chaotic

attractor in a two-cell CNN with an antisymmetric feedback template and a periodic input. Motivated by [Zou & Nossek, 1991], this study addresses the bifurcations andchaos of a two-cell CNN with periodic inputs in a general situation. Indeed,

˙

x₁ =−x₁+ py₁+ sy₂+ bu(t) ,

˙

x₂ =−x₂+ ry₁+ py₂, (1) is studied with the output function

y = f (x) = 1

2(|x + 1| − |x − 1|) , (2) where the feedback template A = [r, p, s], satisﬁes

p > 1, p− 1 < r and p − 1 < −s . (3) The input function (or forcing function), as in [Zou & Nossek, 1991], is

u(t) = sin2π

T t , (4)

with period T > 0 andamplitude b > 0.

The bifurcations of (1) involve ﬁve parameters: r, p, s, T and b. The strategy employedis to begin with b = 0 anda template A = [r, p, s], which sat-isﬁes (3). Without input, we are mainly concerned

with the existence anduniqueness of limit cycle Λ₀.

Λ₀ will interact with inputs bu(t) and may cause

complicateddynamics later. Then a suitable range

of T and b is identiﬁed to ensure that (1) have chaotic attractors.

The template A governs the basic dynamics of (1). When b = 0 and(3) holds, (semi-)stable limit cycles always exist. Indeed, all trajectories, except

the origin, will tendto the limit cycles as t → ∞.

Our numerical experience indicates that a unique limit cycle always applies. See Sec. 4 for details.

The impact of an input bu(t) on its period T andamplitude b are studied. Consider (1) with initial conditions

x₁(0) = ξ₁ and x₂(0) = ξ₂. (5)

The solution of (1) and(5) is denotedby

(x₁(t, ξ₁, ξ₂; b, T, A), x₂(t, ξ₁, ξ₂; b, T, A)) . (6)

The ω-limit set of (6) is denoted by

ω(ξ₁, ξ₂; b, T, A) , (7)

and the nonwandering set of (1) is denoted by

Ω(b, T, A) =

(ξ1, ξ2)∈R2

ω(ξ₁, ξ₂; b, T, A) . (8)

Since the input is T -periodic, for a ﬁxed parameter A, T and b, a two-dimensional Poincar´e map of (1) can be deﬁned as

F (ξ₁, ξ₂) = (x₁(T, ξ₁, ξ₂), x₂(T, ξ₁, ξ₂)) . (9)

Now, the study of the bifurcations problem of (1) is equivalent to the stud y of how Ω(b, T, A) changes when b, T and A vary. To simplify the prob-lem, rather than studying Ω(b, T, A), this paper is concernedmainly with how “typical” trajectories vary with b, T and A. In particular, when b > 0,

the trajectory Γ_b ≡ Γ(b, T, A) of (6) and ω-limit

set ω_b ≡ ω(b, T, A) of (7) with the initial condition

at the origin O = (0, 0) are considered. The ω-limit

set of the Poincar´e map is denoted by ˆω(b, T, A). To

show Ω(b, T, A) is a chaotic attractor, the following conditions must be proven to hold.

(i) Γ(b, T, A) has a positive Lyapunov exponent,

(ii) ˆω(b, T, A) is fractal,

(iii) FFT (Fast Fourier Transform) of Γ(b, T, A) has a broad-band.

An eﬀective approach of studying eﬀects of the input period T is to examine the asymptotic limit

cycle Λ_∞(T, A) by letting b → ∞. When T ≤ T₀∗

(deﬁned in (67)) then Λ_∞ and −Λ_∞ can be

sep-arated. Therefore ω(b, T, A) and −ω(b, T, A) are

separatedfor large b but collide when b is small. If b becomes even smaller, a chaotic attractor may

(3)

develop. Indeed, the onset of chaos induced by crises

of ω(b, T, A) and −ω(b, T, A) were observedfor

suitable T and b. For details, please see Fig. 13 and Sec. 6.

After an interesting range of T is identiﬁed, the eﬀect of b can be examined. Intuitively, the

unperturbedlimit cycle will dominate when b is

small. Indeed, FFT of Γ(b, T, A) is considered when b > 0 andis relatively small. Let T_b be the

periodwith the largest amplitude a₁(b) of FFT on

x₁(t, 0, 0; b, T, A) except for T -mode, and a_T =

a_T(b, T, A) be the amplitude of the period T mode.

The ratio

A(b) ≡aT(b) a₁(b)

(10)

represents the relative strength of the T -mode with

respect to the T_b-mode as b varies. Equation (1) is

called T_b dominant ifA(b) 1, the Tb and T modes

are comparable if A(b) 1 but T is dominant if

A(b) 1.

When T_b is dominant, ω_b is foundto be either

quasi-periodic or periodic. The periodic windows

typically form a devil’s staircase when b∈ (b∗, b∗₀),

where 0 < b_∗ < b∗₀ depends on A and T . For

exam-ple, Figs. 3 and11 present the ZN-case in [Zou &

Nossek, 1991], that is A = [1.2, 2, −1.2] and T = 4.

When T_b and T modes are comparable, the

Lyapunov exponents of Γ(b, T, A) and ω-limit set ˆ

ω(b, T, A) of Poincar´e map are computed. In many interesting cases, including the ZN-case, ﬁnitely many chaotic andwindow regions interweave with each other. In chaotic regions, the largest Lyapunov

exponent is positive andˆω(b, T, A) is fractal, as in

Fig. 15. ˆω(b, T, A) looks like a lady’s shoe as in

the ZN-case andcontains a horseshoe as in asym-metric templates case, as shown in Figs. 20 and 23. In each window, the basic periodic cycle can be identiﬁed, i.e. the periodic cycle with the mini-mum period. In the window, a sequence of period-doubling is observed to the left of the basic periodic cycle. A quasi-periodic region is to the right of the basic periodic cycle.

For b large, the T -mode dominates, the attrac-tors ω(b, T, A) gets simpler, from quasi-periodic to periodic as b increases. See Sec. 6.

The rest of this paper is organizedas follows. Section 2 introduces some properties of solutions of (1) which will be useful later. Section 3 introduces a program for studying bifurcations and chaos since many parameters are involved. The limit cycle of (1) is ﬁrst studied without input. Then, methods

are developed to identify possible ranges of T and b to ensure the occurrence of interesting bifurcation andthe existence of chaotic attractors. Section 4 addresses the existence and uniqueness of the limit cycle of (1) when b = 0 and(3) holds. Section 5 uses the FFT of Γ(b, T, A) to study the

bifurca-tions when b is relative small, i.e. when T_b

domi-nates. Section 6 studies the asymptotic limit cycle

when b→ ∞. Section 7 studies chaos when T_band T

modes are comparable. Section 8 introduces our nu-merical methods. Section 9 brieﬂy discusses results andoﬀers suggestions for future study.

2. Preliminaries

This section provides some preliminary results of (1). Given an initial condition

(x₁(0), x₂(0)) = (ξ₁, ξ₂) , (11)

the solution of (1) with (11) is denoted by

(x₁(t; ξ₁, ξ₂), x₂(t; ξ₁, ξ₂)). We ﬁrst state some

symmetric properties of the solutions of (1).

Theorem 2.1.

(i) If (x₁(t), x₂(t)) is a solution of (1), then

−x1 t +T 2 ,−x₂ t +T 2 (12) is also a solution of (1). In particular, if b = 0, then (−x₁(t), −x₂(t)) is also a solution.

(ii) If b = 0 and (x₁(t), x₂(t)) is a periodic

so-lution of (1), then its period is mT for some positive integer m.

(iii) When b = 0 and A is antisymmetric, i.e. s = −r, if (x1(t), x2(t)) is a solution of (1), then (x₂(t), −x₁(t)) (13) is also a solution. Proof (i) Since f (−x) = −f(x) and u t + T 2 =−u(t) , (14) the function given in (12) is clearly also a solution.

(ii) Assume that (x₁(t), x₂(t)) is a periodic

solu-tion with period ˜T > 0; then x₁(t + ˜T ) = x₁(t)

and x₂(t+ ˜T ) = x₂(t) imply sin(2π/T )(t+ ˜T ) =

sin(2π/T )t for all t. Hence, ˜T = mT for some

positive integer m.

(4)

(iii) When b = 0 and s = −r, let (v₁(t), v₂(t)) =

(x₂(t),−x₁(t)). Then, (v₁, v₂) satisﬁes ˙v₁ =

−v1 + pw₁ + sw₂ and˙v₂ =−v₂+ rw₁+ pw₂,

where w_j = f (v_j) and j = 1, 2. Hence, (13) is

also a solution. The proof is complete.

A set S ⊆ R2 is calledsymmetric with respect

to O = (0, 0), if

−S = S , (15)

where −S = {(−x₁,−x₂) ∈ R2|(x₁, x₂) ∈ S}.

Otherwise, S is calledasymmetric. In

partic-ular, a trajectory (x₁(t; ξ₁, ξ₂), x₂(t; ξ₂, ξ₂)) of

(1) is calledsymmetric if the set Γ(ξ₁, ξ₂) =

{(x1(t), x2(t))|t ∈ R1and(x1(0), x2(0)) = (ξ1, ξ2)}

is symmetric.

The isoclines are useful in studying (1). The

x₂-isocline ˙x₂ = 0 is independent of time, i.e.

h(x₁, x₂)≡ −x₂+ ry₁+ py₂= 0 . (16)

The x₁-isoclines are time periodic with period T if

b > 0, i.e.

g(x₁, x₂)≡ −x₁+ py₁+ sy₂ =−b sin2π

T t , (17) See Figs. 1 and2.

Moving isoclines are ﬁrst usedto discuss the

possible trajectories of (1). When (3) holds and b = 0, the origin O = (0, 0) can be easily veriﬁedto be the only steady-state solution of (1). Furthermore, O is an unstable spiral with eigenvalues λ = (p− 1)± i√−rs. Figure 1 presents vector ﬁelds of (1) when b = 0. In this case, apart from O, all trajec-tories move counterclockwise around O andtendto a limit cycle. See Theorem 4.1 for details. However,

when b > 0, the periodically moving x₁-isocline

g(x₁, x₂) = −bu(t) oscillates horizontally. At a

given instant ¯t, g(x₁, x₂) = −bu(¯t) may

inter-sect h(x₁, x₂) = 0 at point (¯x₁, ¯x₂). In that case,

(¯x₁, ¯x₂) can be regarded as a “temporary or

instan-taneous steady-state”. The trajectories which are

near (¯x₁, ¯x₂) at time ¯t may circle around(¯x₁, ¯x₂)

thereafter. This basic mechanism can generate com-plicatedtrajectories as easily observedfrom the numerical simulations. See Figs. 2 and13(d). The following sections will describe global trajectories.

With reference to a dynamical system of (1), the asymptotic behavior of trajectories as t tends to inﬁnite is of most interest. Therefore, the ω-limit

x₂ x₁ x₂= 1 x₂=−1 x₁= 1 x₁=−1 C+ C− O h(x₁, x₂) = 0 g(x₁, x₂) = 0

Fig. 1. Isoclines and vector ﬁelds of system (1) whenb = 0.

set for each trajectory must be studied. The ω-limit set of (1) and(11) is deﬁnedby

ω(ξ₁, ξ₂) ={(¯x₁, ¯x₂)∈ R2|

∃tk→ ∞ such that xi(tk; ξ1, ξ2)

→ ¯xi, i = 1, 2} . (18)

The nonwandering set Ω of (1) is deﬁned by

Ω =

(ξ1, ξ2)∈R2

ω(ξ₁, ξ₂) . (19)

Note that ω andΩ dependon the template A, T and b. To simplify the notation, the dependency is

omittedif it does not cause confusion. However,

ω_b(ξ₁, ξ₂) or ω(ξ₁, ξ₂; b) and Ω_b may be usedto

emphasize the dependency on b.

The main goal of this paper is to analyze

ω_b(ξ₁, ξ₂) and Ω_b, andto study their bifurcations

as parameters A, T and b vary. The following re-sults can be derivedfrom these isoclines andtheir associatedvector ﬁelds in phase-plane as shown in Figs. 1 and2.

Theorem 2.2. Assume (3) and b ≥ 0. The non-wandering set Ω_b⊆ [−b − p + s, b + p − s] × [−(p + r), p+r]. Furthermore Ω_b is symmetric and attracts all trajectories as t→ ∞.

Proof. For each b ≥ 0, Ω_b ⊆ [−b − p + s, b + p − s]× [−(p + r), p + r] can be easily veriﬁedfrom the associatedvector ﬁelds in phase-plane. See Figs. 1

and2. Clearly, Ω_b attracts all trajectories as t→ ∞.

(5)

By Theorem 2.1(i), Ω_b is symmetric. The proof is

complete.

Remark 2.3. b > 0 may have an asymmetric pe-riodic orbit Λ. In that case, an other asymmetric

periodic orbit −Λ exists, and −Λ ∪ Λ ⊆ Ω_b.

Since the inputs are periodic with period T ,

introducing a two-dimensional Poincar´e map F :

R2_{→ R}2 _by

F (ξ₁, ξ₂) = (x₁(T ; ξ₁, ξ₂), x₂(T ; ξ₁, ξ₂)) (20)

is natural. Clearly the periodic orbits of (1) with period mT are the periodic points of F with period m, andvice-versa.

The ω-limit set ˆω_b(ξ₁, ξ₂) andthe

nonwander-ing set ˆΩ_b of Poincar´e map F can also be studied.

Indeed, ˆ ω_b(ξ₁, ξ₂) ={(η₁, η₂)∈ R2| ∃nk→ ∞ such that Fnk(ξ1, ξ2) → (η1, η₂) as n_k→ ∞} , (21) and ˆ Ω_b = (ξ1, ξ2)∈R2 ˆ ω_b(ξ₁, ξ₂) . (22) Clearly, ˆω_b(ξ₁, ξ₂)⊂ ω_b(ξ₁, ξ₂) and ˆΩ_b ⊆ Ω_b.

Now, the Lyapunov exponents of (1) can be

studied using its Poincar´e map F . Recall that the

Lyapunov exponents of a smooth map F on Rm →

Rm _{are deﬁned as follows [Alligood et al., 1997,}

pp. 194–195].

Definition 2.4. For a smooth map f on Rm, let

J_n = Dfn(v₀), andfor k = 1, . . . , m, let Rn_k be

the length of kth longest orthogonal axis of the

ellipsoid J_nU for an orbit with initial point v₀.

Then Rn_k measures the contraction or expansion

near the orbit of v₀ during the ﬁrst n iterations.

The kth Lyapunov exponent of v₀ is deﬁned by

α_k= lim_n→∞log((Rn_k)1/n), if the limit exist.

In this paper, the system (1) is calledchaotic if the following conditions hold:

(i) the largest Lyapunov exponent of Ωˆ_b is

positive,

(ii) ˆΩ_b is fractal,

(iii) some typical trajectories of (1) have broad-bands under FFT.

Proving that a typical trajectory, say Γ_b(0, 0)

satisﬁes (i) and(ii) suﬃces to verify conditions

(i) and(ii). The following sections present the rele-vant details.

3. Programs for Studying Bifurcations and Chaos

The rest of this paper addresses the bifurcations and chaos of (1) as the parameters A = [r, p, s], T and b vary. The following programs are appliedto study thoroughly a complex andinteresting phenomenon over a range of parameters, since the problems involve ﬁve parameters.

g(x₁, x₂) =−bu(T₄) g(x₁, x₂) =−bu(¯t) g(x₁, x₂) =−bu(3T₄ ) x₂ x₁ x₂= 1 x2=−1 x1= 1 x1=−1 C+ C− O

Fig. 2. Isoclines and vector ﬁelds forb > 0.

(6)

Fig. 3. FFT of the largest 20 modes for the ZN-case:A = [1.2, 2, −1.2] and T = 4.

(I) Take b = 0 andstudy how the sustainedlimit

cycles Λ₀(A) vary with the template A =

[r, p, s]. In particular, examine how the

pe-riod T₀(A) of Γ₀(A) varies with A.

(II) Fix A. Findpossible range of input periods T such that (1) exhibit chaotic behavior for suitable b > 0. In particular, try to ﬁndthe

relation between T and T₀(A) such that (1)

have complex trajectories for some b > 0. (III) Fix A and T obtainedin (I) and(II), try

to identify critical numbers of b, say, b∗₀ <

b∗₁ < · · · < b∗_k, which represent various types of trajectories of (1) andmay cause distinct

bifurcations when b∗_j is crossed.

With reference to program (I), Sec. 4 discusses the existence anduniqueness of limit cycles. To ex-plain how programs (II) and(III) are implemented, a series of numerical experiments with varying b > 0 are presented, as follows.

For ﬁxed A and T , d enote by Γ_b the

forward-trajectory of (1) with the initial condition at the

origin O, and ω_b the corresponding ω-limit set of

Γ_b. Λ₀ is the (inner) limit cycle for b = 0 and is

obtainedfrom Theorem 4.1. Apply FFT to the x₁

-component of Γ_b, i.e. x₁(t; 0, 0), t > 0. Pick up the

ﬁrst N frequencies of these data, i.e. let{akeiωkt}N

k=1

satisfy

|a1| ≥ |a2| ≥ · · · ≥ |aN| ≥ aω| , (23)

for other frequency ω, where a_k = a_k(b) and ω_k =

ω_k(b), denote

τ_k(b) = 2π

ω_k(b), (24)

the periodof the kth mode. For simplicity, denote

T_b = τ₁(b) , (25)

which corresponds to the largest amplitude except for T -mode. It is not diﬃcult to verify

lim

b→0+Tb = T0. (26)

The normalizedcurves

R_k(b) = τk(b)

T_b (27)

of τ_k(b), and1≤ k ≤ N, are very useful for ﬁnding

periodic orbits. To be more speciﬁc, in the ZN-case,

R_k(b) with 1≤ k ≤ 20 and b ∈ [0, 4] are as in Fig. 3.

Figure 3 can be explainedas follows.

(1) The amplitude of the T = 4 mode (represented by a redthick line in Fig. 3) grows steadily as b increases in (0, 3.826). It is comparable to T_b when b is close to 4, near the onset of chaos.

(7)

(a)b∗₁= 3.98 (b)b∗₂= 4.284

(c)b∗₃= 4.365 (d)b∗= 4.2697

Fig. 4. Critical trajectories ofb∗₁,b∗₂,b∗₃ andb∗, whenA = [1.2, 2, −1.2] and T = 4.

(2) Curve number 2 decreases and curve number

3 increases andmerges into T_b/2, giving rise to

4T periodic cycles. The 4T cycle will survive for quite a large range of parameters in (0.43, 0.66). Curves merging is very common andinduces a periodcycle.

(3) The T_b/3 mode maintains the largest

parame-ters in (0, 3.826) andgives rise to a 3T periodic cycle in (1.2, 3.826).

(4) The dotted regions and window regions

(steppedregions) interweave with each other. Steppedregions represent periodic cycles and dotted regions represent quasi-periodic orbits. Section 5 will analyze the bifurcations before the onset of chaos.

In the ZN-case, when b≥ 3.826, the strength of

the T -mode is comparable with or larger than the

strength of the T_b-mode. In the following, a

heuris-tic argument is usedto derive relations among b, T

and T₀ when T_b and T are comparable.

Let

γ(t) = Λ₀(t + t₀) , (28)

Λ₀(t) be the limit cycle of (1) with b = 0. The ﬁrst

equation of (1) is modeled as dx

dt = γ(t) + bu(t) . (29)

Now, γ(t) is a periodic function with period T₀ and

u(t) is a periodfunction with periodT . The two time scales of the functions γ and u can be

nor-malizedto a single time scale τ ∈ [0, 1] by setting

t = T₀τ for γ and t = T τ for u. Hence,

x(t) = x(0) + T₀Γ(τ ) + bT U (τ ), τ ∈ [0, 1] , (30) where Γ(τ ) = _T₀_τ 0 γ(s)ds and U (τ ) = _{T τ} 0 u(s)ds

are normalizedperiodic functions with periodone.

From (30), x(t) can vary maximally if T₀ and bT

(8)

(a)b = 8.44 (b)b = 8.9

Fig. 5. Critical trajectories can circle aroundC+ orC−many times whenA = [1.2, 2, −1.2] and T = 2.

have the same order of magnitudes and an

appro-priate time shift t₀ occurs in (28). Therefore, for a

ﬁxedtemplate A = [r, p, s], let T₀ = T₀(A) be the

periodof the sustainedlimit cycle Λ₀(A) (without

input), deﬁne

b∗(T ) = c0T0(A)

T , (31)

where c₀ ∼ 1 is a constant that depends on A

and T .

From our experience, for a given A and T ,

c₀ = 1 in (31) is a goodestimate for the position at

which to start the search for interesting ranges of

b. c₀(A, T ) may decrease as T increases. In the

ZN-case andmany other templates, (31) workedvery well. See Figs. 4, 15 and19.

The program (II) will be supportedin Sec. 6,

where asymptotic limit cycles Λ_∞ are studied as

b→ ∞.

When the largest Lyapunov exponent of the

Poincar´e-map is close to zero andthen becomes

pos-itive, (1) enters a chaotic region. In the ZN-case,

eight chaotic regions C_k, 1 ≤ k ≤ 8 can be id

en-tiﬁed, followed by successive window regions W_k,

1≤ k ≤ 8. See Fig. 15. Bifurcations, typically back-ward period-doubling in window regions were dis-covered. See Fig. 16. Section 7 provides details. In

the chaotic regions, three critical trajectories at b∗₁,

b∗₂, and b∗₃ can be id entiﬁed which are important to

implement program (III). See Figs. 4(a)–4(c). The

parameter b∗₁ is relatedto the onset of chaos while

b∗₂ and b∗₃ are relatedto fully-developedchaos. They

occur when T is near four in the ZN-case, but b∗₂and

b∗₃ may change considerably when T is relatively

small, say T ≤ 2. In that case, b∗₂ and b∗₃ d o not

exist. Instead, the trajectories may circle around

C+ and C− many times. See Fig. 5 andSec. 7 for

details.

When b is relatively large, in the ZN-case b ≥

4.432, the T -mode dominates, i.e. the sum of the strength of all other modes is less than a few per-cents of T -mode. The chaotic regions disappear and

are followedby quasi-periodic regions andthen,

eventually, a periodic region. Now, Λ_b is either a

symmetric or an asymmetric periodic cycle depend-ing on T. Section 6 settles this issue by considerdepend-ing

the asymptotic limit cycle Λ_∞ as b→ ∞.

4. Limit Cycles

This section addresses the existence and multiplic-ity of limit cycles of (1) when b = 0 and (3) holds.

The existence of limit cycles can be easily proven.

Theorem 4.1. Assume (3) and b = 0, limit cycles exist. Moreover, apart from O = (0, 0), all trajecto-ries will tend to one of the limit cycles as t→ ∞. Proof. Under the assumptions (3), the origin O = (0, 0) can be easily veriﬁedto be the only steady state of (1); moreover, O is an unstable spiral. Indeed, the associated eigenvalues at O are given by

λ_± = p− 1 ± i√−rs .

By Theorem 2.2 andthe Poincar´e–Bendixson

Theorem, a limit cycle exists. Apart from the ori-gins, all trajectories will tendto one of the limit

cycles as t→ ∞. The proof is complete.

(9)

I M1 M2 M3 M4 T1 T2 T3 T4 x₂ = 1 x₂ =−1 x₁ = 1 x₁ =−1 (α,-1) (α₁,1) (1,β₂) (-1,β₃) (α₄,1)

Fig. 6. A typical orbit of (1) with initial condition (α, −1), and 1≤ α ≤ p − s.

Since the nonlinear output function is piecewise linear in (2), the phase-plane can be divided into nine regions which are the mosaic (saturated)

re-gionMj, the transitional (partial saturated) region

Tj andthe interior (not saturated) regionI, j = 1,

2, 3, 4. See Fig. 6.

It is easy to see that the periodic orbit does not

lie entirely in the interior regionI. Therefore,

peri-odic orbits have to intersect the exterior region E,

here

E = R2_{− I =} 4

k=1

(T_k∪ M_k) . (32)

A periodic orbit Λ is called an exterior periodic

cycle (exterior cycle) if Λ⊆ E, otherwise Λ is called

a nonexterior periodic cycle, i.e. ΛI = ∅.

Now, the multiplicity of the exterior periodic cycles can be proven as follows.

Theorem 4.2. Assume (3) and b = 0. No more than two limit cycles are present in the exterior regionE.

Proof. Periodic solutions as in [Thiran, 1997] are constructedto show that no more than two

peri-odic orbits exist in exterior region E.

Now starting at the point (α,−1) at t = 0,

where 1 ≤ α ≤ p − s, the trajectory Γ_α in T₁ is

followed; it intersects x₂ = 1 at the point (α₁, 1) on

t = t₁, 1 < α₁, entersM₁; then intersects x₁ = 1 at

(1, β₂) on t = t₂, entersT₂; then intersects x₁ =−1

at the point (−1, β₃) on t = t₃, andﬁnally enters

M2 andintersects x2 = 1 at the point (α4, 1) on

t = t₄, i.e. (x₁(0), x₂(0)) = (α, −1) , (x₁(t₁), x₂(t₁)) = (α₁, 1) , (x₁(t₂), x₂(t₂)) = (1, β₂) , (x₁(t₃), x₂(t₃)) = (−1, β₃) , (x₁(t₄), x₂(t₄)) = (α₄, 1) . (33)

See Fig. 6. Since b = 0, (1) is an autonomous equation. The periodic orbit cannot intersect

it-self. Therefore, by Theorem 2.1(i), Γ_α is a periodic

(closed) orbit if and only if

α₄ =−α . (34)

α₁, β₂, β₃, α₄ and t₁, t₂, t₃, t₄ must be

com-putedin terms of α. The following expressions can be straight-forwardly obtained. The details are omittedhere. Denote by ξ = r + 1− p, η = r + p − 1 , γ = 1− s − p, δ = p − s − 1 , and q = 1/(p− 1) . (35) Then, α₁ = p +s(1− r) p + α− p +s(r + 1) p ξ η _q , (36) β₂ = p + r− ηγ α₁− p − s, (37) β₃ = p−r(s + 1) p + β₂− p − r(1− s) p _γ δ _q , (38) α₄ = δξ β₃+ r− p+ s− p . (39)

α₄is written as a function of α to show that (34) has

at most two solutions for α ∈ [1, p − s]. Indeed, in

the following, k_i, i = 1, . . . , 17, are constants that

depend on p, r, s, but are independent of α, α₁= k₁α + k₂, β₂= k4 k₁α + k₃ + k5, β₃= k₆β₂+ k₇, α₄= k9 β₃+ k₈ + k10,

(10)

and

α₄ = k9α + k13

k₁₁α + k₁₂ + k14. (40) Substituting (40) into (34) yields, a quadratic equation for α, i.e.

k₁₅α2+ k₁₆α + k₁₇= 0 . (41) Therefore, (34) has at most two solutions in

[1, p− s]. The proof is complete.

The uniqueness of limit cycle in the

exte-rior region E can be proven by making further

assumptions:

Theorem 4.3. Assume (3) and b = 0. If 1 < p≤ 2 then (1) has at most one limit cycle in the exterior region E.

Proof. Note that ∂ ∂x₁(−x1+ py1+ sy2) + ∂ ∂x₂(−x2+ ry1+ py2) = p− 2 if (x₁, x₂)∈ T_i, −2 if (x₁, x₂)∈ M_i,

1 ≤ i ≤ 4. The sign is nonpositive if p ≤ 2. The

Dulac criteria rule out the secondclosedorbit inE.

The proof is complete.

The existence andnonexistence of periodic

cycles in the exterior region E can also be proven

by making additional assumptions.

Theorem 4.4. Assume (3) and b = 0. Let ξ, η, γ, δ and q be given by (35).

(i) There is a periodic cycle in the exterior region E if the following conditions are satisﬁed.

(E₁) p− 1 +s(1− r) p + 1− p + s(r + 1) p ξ η _q ≥ 0 , (42) p− 1 −r(1 + s) p + 1− p − r(1− s) p _γ δ _q ≥ 0 . (43) In particular, if A is antisymmetric, i.e. −s = r, (E₁) and (E₂) are equivalent to

(E) p(p− 1) + r(r − 1) − [p(p − 1) + r(r + 1)] ξ η _q ≥ 0 . (44)

(ii) There is no periodic orbit in the exterior region E if one of the following conditions holds.

(N₁) p− 1 +s(1− r) p + sξ p ξ η _q < 0 , (45) or (N₂) p− 1 −r(1 + s) p − rγ p _γ δ _q < 0 . (46)

In that case, all periodic cycles necessarily in-tersect the interior region I.

Proof. The existence results are ﬁrst proved. It is easy to verify that if Λ is an exterior

peri-odic cycle then Λ{(x₁,−1)|x₁ ∈ [1, p − s]} = ∅.

From (36) and(42),

α₁(α)≥ 1 for all α ∈ [1, p − s] . (47) Similarly, from (38) and(43),

β₃(β₂)≥ 1 for all β₂ ∈ [1, p + r] . (48)

Therefore, x₁(α, 1) maps [1, p−s] into [−p+s, −1].

It implies −x₁(α, 1) maps [1, p− s] into itself and

then has a ﬁxedpoint in [1, p− s]. Hence, (34) has

at least one solution in [1, p− s]. This proves that

exterior periodic cycle exists.

Clearly, (E₁) and (E₂) is equivalent to (E) when

s =−r.

Finally, from (36) and(45),

α₁(α) < 1 for all α∈ [1, p − s] , (49)

andfrom (38) and(46),

β₃(β₂) < 1 for all β₂ ∈ [1, p + r] . (50)

Hence, there no exterior periodic cycle exists. The

proof is complete.

Notably, (45) and(46) can be replacedby stronger conditions that can be veriﬁed easily as follows. 0 < p− 1 < r < 1 and − s ≥ p(p− 1) 1− r , (51) and 0 < p− 1 < −s < 1 and r ≥ p(p− 1) 1 + s . (52)

Figure 7 shows some typical exterior periodic cycles andnonexterior periodic cycles.

Although, the expressions (36)–(39) are ex-plicit, analytically showing that (41) has a unique

solution in exterior regionE remains quite diﬃcult.

(11)

(a)A = [1.2, 2, −1.2], b = 0 (b)A = [0.5, 1.2, −0.5], b = 0

(c)A = [0.5, 1.217, −1], b = 0 (d)A = [0.5, 1.217, −2], b = 0 Fig. 7. Typical limit cycles of (1) withb = 0, (a) exterior cycle, (b)–(d) nonexterior cycles.

However, numerical results indicate that a unique

limit cycle exists in the entireR2 over a quite large

parameter range.

Numerical computation also reveals, the period

function T₀(r, p, s) of the limit cycles for template

A = [r, p, s] that satisﬁes the condition (3) is

de-creasing with respect to r and −s, andincreasing

with respect to p. See Figs. 8 and9. An analytic study of this monotonicity has some progress.

5. Bifurcations Precede Chaos

This section considers the bifurcations before chaos when the amplitude b of the input is relatively small. Section 3 explains the methods used.

Given a template A which satisﬁes (3), assume

that (1) has a unique limit cycle Λ₀(A) with period

T₀ = T₀(A). For each T ∈ (0, T₀) and b > 0, let

T_b = T (b, T, A) be the periodof the largest

am-plitude a₁(b, T, A) of the FFT except for T -mode

appliedto x₁(t, 0, 0; b, T, A) and let a_T(b, T, A) be

the amplitude of T -mode. Let b∗₀> 0 such that

|a1(b, T, A)| > |aT(b, T, A)| (53)

holds in (0, b∗₀), b∗₀ can be assumedto be the least

upper boundof ˜b such that (53) holds in (0, ˜b).

Consider the ﬁrst N highest modes and plot the curves

R_k(b)≡ τk(b)

T_b (54)

for b∈ (0, b∗₀) and k = 1, . . . , N , where τ_k(b) is the

periodof the kth-largest amplitude of the FFT and

T_b = τ₁(b). See Fig. 3. The R_k(b) curve is well

de-ﬁnedlocally in the window regions andcan merge with its neighbor curves. After they merge, they are

considered to be one curve. See curves2 and3 in

Fig. 3.

When a periodic window appears on the open

interval B ⊂ (0, b∗₀), the curve R_k(b) will be

a horizontal line, or approximately one, on B.

(12)

Fig. 8. The period functionT₀(r, p, −r) with b = 0.

Fig. 9. The period functionT₀(r, 2, s) with b = 0.

Furthermore, on B

T_b = m

nT (55)

for some positive integers m and n and(m, n) = 1, i.e. m and n are relative prime. Therefore,

B_m,n = b∈ (0, b∗₀)|Tb = m nT (56) is deﬁned. Denote by [T₀/T ] = m∗, (57)

where [x] is the largest integer which is equal to or smaller than x.

The solutions of (1) can be written explicitly

on each of the nine regions Mj, Tj and I.

There-fore, the exact solutions of periodic orbits in B_m,n

can be rigorously checkedusing a computer

pro-vided n is not too large. The periodic cycles in B_m,n

are of period mT andcircle aroundthe origin O n-times (n-copies). This explanation partially proves the following results.

Conjecture 5.1. Assume (53) holds. Then

(i) B_m∗_,1 = ∅, i.e. a stable m∗T periodic cycle of

(1) exists.

(ii) If (1) has another stable limit cycle with m_∗T

period in B_m_∗_,n_∗ ⊂ (0, b∗₀) and m_∗/n_∗ < m∗, then ∪_m_∗_/n_∗_≤m/n≤m∗B_m,n is open and dense

in (ˆb₁, ˆb₂), where ˆb₁ = inf B_m∗_,1 and ˆb₂ =

sup B_m_∗_,n_∗, i.e. T_b is a function of b is a devil’s

staircase in b. See Fig. 10.

(a)T = 4 (b)T = 2

Fig. 10. Devil’s staircase like period functionT_b,A = [1.2, 2, −1.2], (a) T = 4 and (b) T = 2.

(13)

Figure 11 plots the typical periodic orbits in

B_m,n in the ZN-case andFig. 12 plots the typical

quasi-periodic orbits in the ZN-case.

6. Asymptotic Limit Cycles for Large Inputs

This section addresses asymptotic periodic cycles for various T when b is large. Whether the ω-limit

sets ω_b and −ω_b can be separatedfrom each other

by the x₁-axis such that one lies in the upper half

of phase-plane andthe other lies in the lower half of phase-plane is the main concern. The answer is aﬃr-mative when T is relatively small. See Theorem 6.1. In any case, the system can always support a limit-ing cycle even for large T .

For a given template A = [r, p, s], w₁ = x1

b (58)

(a)b = 0.5, T_b= 4T (b)b = 0.78, T_b=15₄T

(c)b = 0.9, T_b= 11₃T (d)b = 1, T_b=7₂T

(e)b = 1.1, T_b= 10₃T (f)b = 1.2, T_b= 3T

Fig. 11. Some typical orbits inB_m,nprior to chaotic regions,A = [1.2, 2, −1.2] and T = 4.

(14)

(a)b = 0.37 (b)b = 0.37

(c)b = 1.14 (d)b = 1.14

Fig. 12. Some typical quasi-periodic orbit (a) and (c) and their ω-limit set ˆω of the Poincar´e map (b) and (d), A = [1.2, 2, −1.2] and T = 4.

is written as b→ ∞, andthe limiting equation for

w₁ is

dw₁

dt =−w1+ u . (59)

The solutions of (59) are

w₁(t) = ce−t+ 1

1 + Ω2(sin Ωt− Ω cos Ωt) , (60)

where

Ω = Ω(T ) = 2π

T and c is a constant. Consequently,

x₁(t, ξ₁, ξ₂; b)∼ b

1 + Ω2(sin Ωt− Ω cos Ωt) (61)

for large b and t. Notably,

−p − r < x2(t, ξ1, ξ2; b) < p + r (62)

always holds for large t. Now, (1) is assumed to

have a asymptotic limit cycle Λ_∞ with period T as

b→ ∞. From (61), Λ_∞will always almost be in the

region |x₁| ≥ 1. In the limit, denoted by w₂(t) for

x₂(t; b), w₂ satisﬁes dw₂ dt =−w2+ p + r if w2 ≥ 1 , (63) dw₂ dt = (p− 1)w2+ r if |w2| ≤ 1 , (64) dw₂ dt =−w2+ r− p if w2 ≤ −1 . (65)

for a total time of T /2 in the region w₁ ≥ 0. Similar

equations holdin region w₁ ≤ 0 for another T/2

time. The separation theorem is statedas follows.

Theorem 6.1. The system (1) can support a limiting limit cycle Λ_∞ with period T provided

(i) in region w₂ ≤ −1 if

T < T₁∗≡ 2 log 2r

r + 1− p, (66)

(15)

(a)b = 10 (b)b = 4.8

(c)b = 4.428 (d)b = 4.336

Fig. 13. Crises induced byω_b and−ω_b whenA = [1.2, 2, −1.2] and T = 4.

(ii) in region w₂ ≤ 0 if T < T₀∗ ≡ 2 log 2r r + 1− p+ 1 p− 1log r r + 1− p . (67) Similarly, −Λ_∞ lies in w₂ ≥ 1 and w₂ ≥ 0, respectively.

Proof

(i) Assume that Λ_∞remains in the region w₂ ≤ −1

for T /2 time. Then general solutions of (65) are w₂(t) = ce−t+ r− p . (68) For 1≤ β ≤ α < p + r, assume w₂(t₀) =−α , and w₂ t₀+T 2 =−β . (69)

Then (68) and(69) imply T = 2 logα + r− p

β + r− p. (70)

Since 1≤ β ≤ α < p + r, (66) follows.

(ii) Assume that Λ_∞ remains in the region w₂ ≤ 0

for T /2 time. Let

0≤ β ≤ 1 ≤ α < p + r , w₂(t₀) =−α , w₂(t₀+ T) =−1 , w₂ t₀+T 2 =−β , (71)

where 0 < T < T /2. Then from (64) and

(65),

w₂(t) = c₁e−t+ r− p for t ∈ [t₀, t₀+ T] ,

(72)

(16)

Fig. 14. ω_b=−ω_b whenA = [1.2, 2, −1.2] and T = 10. and w₂(t) = c₂e(p−1)t − r p− 1 for t∈ t₀+ T, t₀+T 2 . (73) From (72), T = logα + r− p 1 + r− p. (74)

Similarly, from (71) and(73), T 2 − T ₌ 1 p− 1log r− β(p − 1) r + 1− p . (75) From (74) and(75), T = 2 logα + r− p r + 1− p + 1 p− 1log r− β(p − 1) r + 1− p which implies (67).

The proof is complete.

Remark 6.2

(i) The perturbation methodcan be usedto prove

that there exists a limit cycle Λ_b with period

T for large b. Λ_b lies in the region according to Theorem 6.1. The details are omitted here. See Fig. 13.

(ii) For large T , Λ_∞can be proven to exist with

pe-riod T . Furthermore, Λ_∞ spends most of T /2

time near p + r. Therefore, for large T andlarge

b, Λ_b is a symmetric T -periodic cycle like a

rhombus, with two vertices close to (1, p + r) and(−1, −p − r), respectively. See Fig. 14. Remark 6.3. When T satisﬁes either (66) or (67),

Λ_∞ and −Λ_∞are separated. Then Λ_b and −Λ_b are

also separatedwhen b is large. For example, in the

ZN-case i.e. A = [1.2, 2,−1.2],

T₁∗ = 4.9698 (76)

and

T₀∗ = 8.5533 . (77)

Note that T = 4 < T₁∗ has been usedin [Zou

& Nossek, 1991], when b decreases to some

criti-cal number b∗_∞, ω_b and −ω_b cross each other and

cause crises; chaos may occur when b decreases a little further. See Figs. 13(a)–13(d). If T is too large, chaos may not occur, for example, in the ZN-case A = [1.2, 2,−1.2] and T = 10, in that case, ω_b is like rhombus. See Fig. 14.

7. Chaos

This section considers the chaotic phenomena that

occurs when the strengths of Γ_b andthe input bu

are comparable. A specific model of the ZN-case is studied first to elucidate the methods of the study and the chaotic behavior. Section 7.1 addresses bi-furcation to chaos by fixing the template A and

the input period T andvarying the amplitude b.

Section 7.2 considers the eﬀect of an input period T by ﬁxing the template A. The asymptotic limit

cycles Λ_∞ studied in Sec. 6 will guide the approach

taken to solving the problem. Section 7.3 addresses the fundamental role of the template A.

7.1. Eﬀects of input amplitude

As statedin Sec. 3, ω_b is a chaotic attractor if the

following three conditions are satisﬁed.

(ii) ˆω(b, T, A) is fractal,

(iii) Γ(b.T, A) has a broad-band in FFT.

The numerical methods for computing

Lyapunov exponents, Poincar´e maps andFFTs

are standardandhave been appliedby many

researchers. Section 8 will detail these methods. In all cases studied, the largest Lyapunov exponents

must exceed0.02 to be consideredto be positive.

When ˆω_b appears as a partial lady’s shoe or as a

horseshoe, it is considered to be fractal. Quantita-tive results concerning fractal dimensions can also

(17)

Fig. 15. Lyapunov exponents diagram for the ZN-case with b in (30) and (37), b∗₁ = 3.98, b∗₂ = 4.284, b∗₃ = 4.365 and b∗= 4.2697.

(a) Type I (b) Type II

Fig. 16. (a) Type I: the orbit Γ_b circles around all three points C+, O and C−. (b) Type II: the orbit Γ_b does not circle around all of three pointsC+,O and C−.

be computedin these cases. Finally, the broad-band of FFT is considered in the classical sense.

The ZN-case, with A = [1.2, 2,−1.2] and T =

4, is ﬁrst considered as a model to help discuss the bifurcations of chaotic phenomena. The Lyapunov exponents were computedfor a long b > 0.

Table 1. Characteristics (1) (2) (3) (4) (5) W0 (3.956, 3.96) 7 a T I W1 (3.992, 4.008) 11 s T_b I W2 (4.052, 4.068) 6 a T I W3 (4.092, 4.104) 10 a T I W4 (4.124, 4.168) 4 a T_b I W5 (4.252, 4.260) 9 s T I W6 (4.368, 4.412) 8 a T I W7 (4.42, 4.431) 4 a T II W8 (4.433,∞) 2 a T I

The largest Lyapunov exponent that is close

to or above zero is recorded in b ∈ (3.8, 4.6). See

Fig. 15. It can be usedto identify eight regions

C_k, 1 ≤ k ≤ 8, which are chaotic since they have

positive Lyapunov exponents. Notably, C₄ is the

region that has been studied by Zou and Nossek

[1991]. Each chaotic region C_k is followedby a

window region W_k and1 ≤ k ≤ 8. The window

W₀ precedes C₁.

In each W_k and0 ≤ k ≤ 8, the basic

peri-odic cycle — the periperi-odic solutions with the small-est period — can be identiﬁed. These windows are ﬁrst comparedin terms of the following

character-istics of the basic periodic cycle in each W_k.

(1) Range of parameters in window, (2) Periodin T units,

(3) Symmetry: “s” for symmetric and“a” for asymmetric cycles,

(18)

Fig. 17. Bifurcations diagram for ZN-case withb in window W4= (4.124, 4.172).

(4) Dominating mode in FFT: T_b for superiority of

T_b and T for input,

(5) Type of attractor: “I” for type I, “II” for type II. See Figs. 16(a) and16(b).

Considering W_k carefully, for example, in W₄,

reveals that at the middle point of (4.124, 4.168) the 4T basic periodic cycle is asymmetric. To its left, a sequence of periodic-doubling occurs; to its right is

a quasi-periodic region. See Fig. 17. Similarly, W₅

includes a symmetric 9T basic periodic cycle, with periodic-doubling to its left and a quasi-periodic region to its right.

To analyze ω_b in C_k, it is considered to be

a chaotic bifurcation from diﬀerent types of

win-dows W_k−1 and W_k. ˆω_b of the Poincar´e map in C_k

will track the change of the basic periodic orbits ˆ

ω_b in W_k−1 and W_k. See Figs. 18(a)–18(n) for the

(a)b = 3.976 ∈ C₁ (b)b = 3.996 ∈ W₁

(c)b = 4.04 ∈ C₂ (d)b = 4.056 ∈ W₂

Fig. 18. Typical Poincar´e section in chaotic regions and basic periodic cycles in windows,A = [1.2, 2, −1.2] and T = 4.

(19)

(e)b = 4.084 ∈ C₃ (f)b = 4.104 ∈ W₃

(g)b = 4.112 ∈ C₄ (h)b = 4.144 ∈ W₄

(i)b = 4.204 ∈ C₅ (j)b = 4.256 ∈ W₅ Fig. 18. (Continued )

(20)

(k)b = 4.264 ∈ C₆ (l)b = 4.412 ∈ W₆

(m)b = 4.416 ∈ C₇ (n)b = 4.42 ∈ W₇ Fig. 18. (Continued )

transitions. Notably, the basic periodic cycles in

windows W_k like those in Fig. 18, which circle

around C+, O and C−, diﬀer greatly from the

pe-riodic cycles in B_m,n as shown in Fig. 11.

Notably, in the chaotic regions, the time be-tween a maximum point anda minimum of point

of x₁(t, 0, 0, b) is approximately T /2. Tracing

the maximum andminimum points of the x₁

-coordinates of ω_b, yields an approximate Poincar´e

T /2-map, ˆ F (ξ₁, ξ₂)≡ x₁ T 2; ξ1, ξ2 , x₂ T 2; ξ1, ξ2 , (78)

which is double the Poincar´e T -map F given in (9).

7.2. Impacts of the input periods

This subsection brieﬂy discusses the eﬀects of the input period. Section 6, for a given T , discusses

the types of asymptotic limit cycles Λ_∞ that arise

as b → ∞. Two diﬀerent types arise according to

separability. Λ_b and −Λ_b is separable if −Λ_b = Λ_b

andseparatedby x₁-axis, which is ensuredwhen T

is relatively small. See Theorem 6.1. Otherwise, Λ_b

and−Λ_b are inseparable.

The numerical evidence indicates that separa-bility is importantly involvedin the occurrence of

chaotic attractors. Apparently, separable Λ_b and

−Λb cause crises as b decreases to a particular

threshold, for example, near b∗₃ in the ZN-case.

In the search for chaotic regions, Theorem 6.1

is ﬁrst appliedto T ≤ T₀∗ as in (67). For those T ,

b near b∗(T ) = T₀(A)/T is ﬁrst tried. Computing

three critical trajectories at b∗₁, b∗₂ and b∗₃ is

im-portant. Normally, b∗₁, b∗₂ and b∗₃ near b∗(T )

when-ever they exist. In the ZN-case, the graphs of b∗,

b∗₁, b∗₂ and b∗₃ are drawn and compared with the

re-gions in which (1) has a positive Lyapunov exponent (≥ 0.02). See Fig. 19. Notably, the chaotic regions

of (1) are markedby◦ in Fig. 19. Since the chaotic

regions andwindows regions interweave each other,

the marker ◦ in Fig. 19 is necessarily isolatedfor

each T .

(21)

Fig. 19. Critical numbersb∗,b∗_k,k = 1, 2, 3 for varying T , A = [1.2, 2, −1.2].

The computedresults are statedas follows.

Theorem 7.1. Given a template A = [r, p, s] that satisﬁes (3), (1) has a chaotic region on

[b∗₁(T ), b∗₃(T )] provided the input period T satisﬁes

(67) and b∗₁(T ), b∗₂(T ) and b∗₃(T ) exist.

Notably, critical trajectories Γ_b(0, 0) that circle

around O 6-times, C+ m-times and C− n-times at

critical numbers b∗_,m,n can be induced. See Figs. 19

and20 for T = 2 or T = 3. In that case, analogue results similar to Theorem 7.1 can be stated. The details are omitted here.

For T ∈ [3.5, 4.5], chaotic phenomena

sim-ilar to T = 4 were observed. For T ≥ 5, no

chaotic regions are found. For T ∈ [2, 3], chaotic

regions exist but ˆω_b are not a lady’s shoe. See

Figs. 20(a)–20(h).

(a)T = 2 and b = 8.88 (b)T = 2 and b = 8.872

(c)T = 3 and b = 5.844 (d)T = 3 and b = 5.828 Fig. 20. Chaotic attractors and basic periodic cycles forA = [1.2, 2, −1.2] with varying T .

(22)

(e)T = 3.5 and b = 4.94 (f)T = 3.5 and b = 4.901

(g)T = 4.5 and b = 3.86 (h)T = 4.5 and b = 3.822 Fig. 20. (Continued )

7.3. Varying templates

The role of template A = [r, p, s] is fundamental. It governs the basic dynamics among the inputs. A more complete study of the eﬀects of the template is required. Some preliminary results are presented here; very interesting results can be obtainedby varying the templates.

Since the ﬁrst criterion that determines

whether ω_bis a chaotic attractor is that ω_bmust has

a positive Lyapunov exponent, let λ₁(b, T, r, p, s)

be the largest Lyapunov exponent of ω(b, T, r, p, s) anddeﬁne

λ∗₁(r, p, s) = sup

b>0, T >0λ1(b, T, r, p, s) . (79) Then, the template A = [r, p, s] can introd uce a chaotic attractor with a suitable input bu(t) if

λ∗₁(r, p, s) > 0. Insteadof considering all b > 0 and

T > 0 in (79), denote

λ∗(r, p, s) = max{λ₁(b, T, r, p, s)| (80)

δT₀∗ ≤ T ≤ T₀∗ and δ₁b∗(T )≤ b ≤ δ₂b∗(T )} ,

where T₀∗ is deﬁned in (67) and b∗(T ) is deﬁned

in (31), δ and δ₁ are small positive numbers, for

example δ = δ₁ = 0.1, and δ₂ = 2. Numerical

evi-dence suggests that λ∗(r, p, s) closely approximates

to λ∗₁(r, p, s).

Antisymmetric A is ﬁrst considered, i.e. s =−r,

andwrite

λ∗(r, p) = λ∗(r, p, −r) . (81)

Taking p = 2 in (81), the graph of λ∗(r, 2) is

plot-tedfor r ∈ (1, 6). See Fig. 21. Notably, in the

ZN-case, i.e. r = 1.2 is not a maximum point for

λ∗(r, 2). Indeed, when most r are in [3, 4.5], a

higher Lyapunov exponent exists, andcan induce more complex chaotic behaviors. See Fig. 23(c).

Some preliminary results are also obtainedfor asymmetric templates. See Fig. 23. For

asymmet-ric templates, λ∗(r, 2, s) are computedover some

ranges. See Fig. 22. The marker indicates for

(23)

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 1.2

r

λ

∗

(r, 2)

Fig. 21. The maximum Lyapunov exponent function λ∗(r, 2) for A = [r, 2, −r].

Fig. 22. The maximum Lyapunov exponent for λ∗(r, 2, s) maker for λ∗> 0 and · for λ∗< 0.

(a)A = [1.2, 2, −1.3], b = 4.16, T = 4.3 (b)A = [1.5, 2, −1.7], b = 4.7, T = 2.62

(c)A = [4.5, 2, −4.5], b = 6.9883, T = 1.31 (d)A = [1.5, 2, −1.75], b = 4.28, T = 2.08 Fig. 23. Some typical chaotic attractors for generalA = [r, 2, s].

(24)

(e)A = [1.5, 2, −4], b = 9.06, T = 0.88 (f)A = [1.05, 2, −3.95], b = 2.72, T = 6.2 Fig. 23. (Continued )

those A = [r, 2, s] for which λ∗(r, 2, s) is positive

and · is nonpositive. Figure 23 plots some typi-cal chaotic attractors for antisymmetric andgen-eral asymmetric templates. Notably, most are not lady’s shoes. The more detailed results concerning bifurcations andchaos will be reportedelsewhere.

8. Numerical Methods

This section describes the several numerical

meth-ods used herein, including the Poincar´e map, the

FFT andthe Lyapunov exponent.

The trajectory of the system (1) must first be generated. Numerically, for a given set of parame-ters, a template A = [r, p, s] that satisfies (3), an amplitude b andperiodT , the system of differen-tial equations is solvedin FORTRAN 90 by call-ing a subroutine, RKF45, uscall-ing the Runge–Kutta– Fehlberg (4, 5) methods described in [Fehlberg,

1968], with step size = 0.05, absolute error 1×10−10

andrelative error 1× 10−8.

Since the ω-limit set ω(b, T, A) is of greatest

concern, 2× 106 steps are taken in the RKF45

in-tegration. The ﬁrst 1× 106 steps were ignored, and

the following numerical methods applied to the

re-maining data; the last 1× 106 points were taken as

the ω-limit set ω(b, T, A).

The ω-limit set ˆω(b, T, A) of Poincar´e T -map is

taken every T /stepsize points from ω(b, T, A). The

relative error of the Poincar´e map can be easily

com-puted. For example, in the ZN-case T = 4 with a step size 0.05, 80 steps must be integratedfor each

point on the Poincar`e map. Therefore, the relative

error 1× 10−8× 80 = 8 × 10−7 is obtainedfor each

successive point of the Poincar´e map.

The Lyapunov exponents are obtainedby

av-eraging eigenvalues of DF (ξ₁, ξ₂) on each point in

ˆ

ω_b. Here, a convergent condition is imposed that the

relative error is less than 1× 10−4. Moreover, the

ﬁrst 1× 106 steps in the numerical integration are

ignoredto accelerate the convergence.

9. Conclusions

Zou andNossek [1991] discovereda chaotic

at-tractor in a two-cell CNN with template A =

[1.2, 2,−1.2] andinput b sin(2π/T ) with T = 4 and

b ∼= 4.08. This work investigates bifurcations and

chaos for a broadrange of templates A = [r, p, s], input period T > 0 andinput amplitude b > 0.

The bifurcations of (1) involve ﬁve parameters; r, p, s, T and b. The strategy usedherein is to begin with b = 0 and A satisfying (3). Section 4 studied the existence andmultiplicity problems of periodic

cycle of (1). The existence of the limit cycle Λ₀(A)

is proven when (3) holds. The existence and multi-plicity of exterior periodic cycles are studied under further assumptions. The uniqueness problem is still open for general A that satisﬁes (3). The numerical evidence suggests that the limit cycle is unique.

The bifurcation problem is studied by examin-ing how “typical” trajectories vary with b, T and A. In particular, the trajectory Γ(b, T, A) and its ω-limit set ω(b, T, A) with initial conditions at the origin O = (0, 0), are considered. The system (1) is considered to be chaotic, if

(ii) The Poincar´e T -map ˆω(b, T, A) of ω(b, T, A)

is fractal and