Part 1, No. 10, October 2001
2001 The Japan Society of Applied Physicsc
Chaotic Behaviors of Bistable Laser Diodes and Its Application in Synchronization of Optical Communication
Weichung WANG1, Tsung-Min HWANG2, Cheng JUANG3, Jong JUANG4, Chin-Yueh LIU5and Wen-Wei LIN5
1Department of Mathematics Education, National Tainan Teachers College, Tainan 700, Taiwan, R.O.C.
2Department of Mathematics, National Taiwan Normal University, Taipei 106, Taiwan, R.O.C.
3Electronics Department, Ming Hsin College, Hsinchu 300, Taiwan, R.O.C.
4Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan, R.O.C.
5Department of Mathematics, National Tsing Hua University, Hsinchu 300, Taiwan, R.O.C.
(Received April 25, 2001; accepted for publication July 10, 2001)
Period-doubling routes to chaos behavior in bistable laser diodes is examined numerically. An external sinusoidal electronic drive is injected to generate chaotic light output. Phase portrait, Poincar´e map, bifurcation diagram, and Lyapunov exponent are then calculated to assert the existence of chaos. Furthermore, according to Pecora and Carroll’s theory, a drive-response system is constructed by using the bistable laser diodes to mask the sinusoidal electronic signal with chaotic light. Synchronization can be achieved for optical simplex and duplex transmissions using bistable laser diodes. The proposed synchronization scheme between two separate chaotic systems provides a key step toward optical chaotic communication.
KEYWORDS: bistable laser diodes, period-doubling routes to chaos, synchronization, drive-response system, optical chaotic communication
1. Introduction
Optical bistability can be achieved by passive type bistable systems1)or active-type bistable laser diodes.2) In some pas- sive components, optical bistability routes to chaos can be achieved.3–5) However, the study of optical bistability routes to chaos in active components remains sparse. Kawaguchi6) studied bistable operation and self-sustained pulse oscillation in a semiconductor laser. He also predicted the occurrence of chaos in the self-pulsating state of the model.
In this work, we further assert existence of chaos and period-doubling routes to chaos in the bistable region of the model suggested by Kawaguchi7, 8) numerically. We used an electronic-controlled technique to produce chaotic light output from the laser diodes. These techniques involved adding a sinusoidal signal to a DC bias current and were also implemented with other types of laser diodes.9, 10) The chaotic behaviors are investigated by computing phase por- trait, Poincar´e map, bifurcation diagram, and Lyapunov ex- ponent.
Furthermore, the chaotic light output provides a possible light source for a synchronization scheme in secure optical communication. Pecora and Carroll11, 12) demonstrated cri- teria that result in synchronization between two chaotic sys- tems. According to their theory, this work addresses three issues of optical chaotic communication while using drive- response system formed by the bistable laser diodes. First, the bistable laser diodes can be used to mask the electronic signal with optical chaotic light. Second, one state variable between the drive-response system should be transmitted, re- ceived, and coupled. Third, all of the conditional Lyapunov exponents for the drive-response system have to be negative.
The ability that two separate chaotic systems can be syn- chronized exposes an application of chaos toward secure com- munications. For example, Lorenz-based chaotic circuits13) and self-pulsating laser diodes14) were used to achieve the goal. The synchronization scheme of the drive-response system using bistable laser diodes also provides a promis- ing model for chaotic communication. While both the self-
pulsating14) and bistable laser diodes based chaotic systems can be synchronized, differences between the systems exist.
The chaos of the self-pulsating laser diodes system occurs due to the conflict between the injected sinusoidal signal fre- quency and the system natural frequency. The quasi-two pe- riodic routes to chaos was thus identified. In this work, phase portrait illustrates the phenomenon that two stable equilib- rium points attractors compete with each other and thus in- duce the chaotic behavior. Period-doubling routes to chaos in the bistable laser diodes system are also spotted.
This article is organized as following. Section 2 describes the three-dimensional bistable laser diodes model. Section 3 demonstrate numerical simulation chaotic behaviors from dif- ferent perspectives. This section also demonstrate the period- doubling routes to chaos. We first inject a sinusoidal elec- tronic signal into the bistable laser diodes. Phase portrait, Poincar´e sections, bifurcation diagram, and Lyapunov expo- nents are all then computed to show various nonlinear dynam- ical effects. Using the bistable laser diodes as chaotic mask- ing devices, §4 considers a drive-response system in optical chaotic simplex and duplex transmissions according to Pec- ora and Carroll’s theory. Synchronization of the simplex and duplex transformations using the bistable laser diodes is also studied in this section. Section 5 concludes the paper finally.
2. Formulation of Bistable Laser Diodes
Kawaguchi described the simplest model for absorptive bistability by the following rate equations:7, 8)
dne1
dt = P1− Bn2e1− g1(ne1)npνg− ne1
τnr1
, (2.1)
dne2
dt = P2− Bn2e2− g2(ne2)npνg− ne2
τnr 2
, (2.2)
and dnp
dt = npνg[γ1g1(ne1) + γ2g2(ne2)]
−np
τp
+ βS PB(γ1n2e1+ γ2n2e2).
(2.3)
5914
Subscripts 1 and 2 in ne, P, g,τnr, andγ refer to parameters in the gain and absorption regions, respectively. The parameter nei is the numbers of injected carriers in the i th region and npis the numbers of photons in the laser resonator. Pi is the pump rate per unit volume and is defined as Pi = Ii/qV , where Ii is the total input current in the i th region and V is the volume. B is the recombination coefficient, gi(nei) is the gain function,νgis the group velocity,βSPis the spontaneous emission coefficient. The parameterτnriis the carrier lifetime in the i th region andτpis the photon lifetime. The parameter γi is a ratio of the length of the region i to the whole cavity length, andγ1+ γ2= 1.
The rate equations (2.1)–(2.3) can be rewritten as following dimensionless form15)
ηd Ne1
d T = Sp1− α1Ne12
− (α2Ne12 + Ne1+ α3)Np− Ne1,
(2.4)
ηd Ne2
d T = Sp2− α1Ne22
− (α2Ne22 + Ne2+ α3)Np− α4Ne2,
(2.5)
d Np
d T = Np[γ1(α2Ne12 + Ne1+ α3) + γ2(α2Ne22 + Ne2+ α3)]
− Np+ ε(γ1Ne12 + γ2Ne22),
(2.6)
where the dimensionless variables are defined by
Ne1= bνgτpne1, Ne2= bνgτpne2, Np = bνgτnr1np, Sp1 = bνgτpτnr1P1, Sp2= bνgτpτnr1P2,
α1= Bτnr1
bνgτp
, α2 =a b
1 bνgτp
, α3= cνgτp, α4= τnr1
τnr 2
.
The new parameters are defined as follows. Let the re-scaled time T be equal to t/τp, the key parameters η be equal to τnr1/τp andδ = 1/η, and the normalized spontaneous emis- sion coefficientε be equal to BβSPτnr1/bνgτp. To carry out numerical simulations, we used the parameter values listed in Table I. These values were also used in.15)The dimensionless forms of the rate equations and these values are used through- out in this paper.
Optical bistability can occur in the bistable laser diodes.
Table I. Parameters used to simulate bistable laser diodes.
Parameter Value
P2 0
γ1 0.5
γ2 0.5
α1 0.532192
α2 0.106838
α3 −4.12352
α4 1
δ 3.73 × 10−4
ε 0.532192 × 10−5
Sp1 28.5
Sp2 0.0
10 20 30 40
Sp1 -5
-1 3
Np (log)
Fig. 1. The bistable operation shown by photon density Npversus pump rate Sp1.
Figure 1 shows bistability diagrams for the case ofγ1 = 0.5.
The lines show computed dimensionless photon intensity Np
plotted against the dimensionless pump rate Sp1. Asγ1= 0.5, the bistable laser diodes model has a bistable region in Sp1 ∈ [25.00, 26.80]. The “bistable region” here denotes a region in which two stable equilibrium points are associated with a single dimensionless pump rate. The bistable phenomenon presents a certain degree of complexity: a single input state is associated with two stable optical output states; and the initial condition of a certain input state affects which of the stable states the system would converge to. We are thus motivated to see if chaos exists when the model laser diodes presents bistability.
3. Chaos in Bistable Laser Diodes
We examine the chaotic behaviors in the model by apply- ing an electronic-controlled external drive in the gain region described by (2.4). The external sinusoidal drive is modulated as
I = Sp1+ mcsin(2π · mf· T ), (3.1) where mcis the dimensionless modulation current, and mfis the dimensionless modulation frequency.
We now demonstrate that suitable chosen parameters can lead to chaos. We observe that the attractors keep chang- ing in the bistable region. That is, two stable equilibrium points attractors compete with each other when the value Sp1+ mcsin(2π · mf· T ) sits in the bistable region. We thus hope that the nonlinearity of the model and such complexity result in chaos. Following sections assert the occurrence of period-doubling routes to chaos by phase portrait, Poincar´e map, bifurcation diagram, and Lyapunov exponents.
3.1 Phase portrait
We first examine the phase portraits of the dynamical sys- tem in Fig. 2. In the figure, the dotted lines represent the com- puted stable equilibrium points showing bistable operation.
The solid lines are the phase portraits projected onto the pump rate (Sp1) space. The y-axis shows photon density. The x-axis denotes the dimensionless pump rate value combined with the external sinusoidal drive, namely Sp1+ mcsin(2π · mf · T ).
Parameters Sp1 = 28.0 and mf = 3.73 × 10−3 are used in
Fig. 2. The period-3, -6, -12 and chaos phase portraits. The dimensionless modulation current mc equals (a) 8.00, (b) 8.30, (c) 8.43 and (d) 8.50, respectively. The dimensionless modulation frequency mf= 3.73 × 10−3 for the four cases.
the tests. In the sub-figures of Fig. 2, mcequals to (a) 8.00, (b) 8.30, (c) 8.43, and (d) 8.50. As shown in the figure, sub- figure (a), (b), (c), (d) presents period-3, period-6, period-12, and chaos, respectively. Note that although all four phase por- traits pass through the bistable region, only part (d) demon- strates chaos.
3.2 Poincar´e map
The asymptotic behavior of the system can be characterized by Poincar´e mapping. Let Poincar´e section be a properly chosen(n − 1)-dimensional hyperplane, such that the hyper- plane is not tangent to the trajectory and does not contain tra- jectory planes. Letφt(x) be the trajectory under observation.
The trajectory will repeatedly pass through from one side to the other. The Poincar´e map can thus describe the discrete intersection points. In other words, let x ∈ and the initial time be t0. The Poincar´e map P(x) : → is defined by P(x) = φt+tx(x, t0). P(x) denotes the point at which the tra- jectoryφt+tx(x, t0) first returns to at time txfor the point x emanating from at time t0. For further information, see the works of Alligood et al.16)and Parker and Chua.17)
Let a finite set {x1, . . . , xk} denote the points at which trajectory φt(x) intersects hyperplane . The system will has a period-k solution. To illustrate this concept, we consider a two-dimensional hyperplane that passes through point (4.98, 0.00, 0.00) and has normal direction [−1.00, 0.00, 0.00]. The hyperplane is then applied to the three-dimensional Ne1–Ne2–Np phase space of the system.
Figure 3 presents the limit sets of the Poincar´e map described by the hyperplane for (a) mc= 8.00 (period-3), (b) mc= 8.30 (period-6), (c) mc = 8.43 (period-12), and (d) mc = 8.50 (chaos). As explained above, the trajectory repeatedly passes through the hyperplane at three separate points in the period- 3 case. The period-6 and 12 cases are similar. In the chaotic case, on the other hand, the intersection points do not form a simple geometrical object. Note that more than 500 points are plotted on each of the diagrams in Fig. 3.
(a)
0 25 50
Np
(b)
(c)
2 2.1 2.2
Ne2 0
25 50
Np
(d)
2 2.1 2.2
Ne2
Fig. 3. The Poincar´e sections. In part (a) mc = 8.00 (period-3), (b) mc = 8.30 (period-6), (c) mc = 8.43 (period-12) and (d) mc = 8.50 (chaos). For four sub-figures, mf= 3.73 × 10−3.
3.3 Bifurcation diagram
We next explore chaotic behavior using a bifurcation diagram of mc. We first find the limit set points of the Poincar´e map on the hyperplane that passes through point (4.98, 0.00, 0.00) and has normal direction [−1.00, 0.00, 0.00]. The Poincar´e map points are then pro- jected onto the Npaxis for each mc. Figure 4 shows the bifur- cation diagram with photon intensity Np plotted against the dimensionless modulation current mc. The parameter mcis increased from 8.00 to 9.00 at a parameter resolution 0.005.
Examining the bifurcation diagram, we see that the period- 3 closed orbit is indicated by three points on the hyperplane for 8.00 ≤ mc ≤ 8.24. As mcpasses through 8.24, a sta- ble period-6 closed orbit is spawned. As mcis increased fur- ther, a period-12 and a period-24 closed orbit are spawned when mc= 8.395 and mc= 8.435, respectively. The period- doubling process continues until mc = 8.45, at which the system becomes chaotic.
Fig. 4. Bifurcation diagram of a Poincar´e map with Npplotted against di- mensionless modulation current mc. The dimensionless modulation fre- quency mf= 3.73 × 10−3.
3.4 Lyapunov exponents
We further investigate the stability of periodic solutions of the rate equations by computing Lyapunov exponents. Lya- punov exponents estimate rates of change when nearby or- bits converge or diverge. They act as indicators determin- ing the stability of equilibrium points, periodic solutions, quasi-periodic solutions, and chaotic behavior. Lyapunov ex- ponents are a generalization of eigenvalues at equilibrium points and of characteristic multipliers. For an n-dimensional continuous-time system, Lyapunov exponents can be defined as follows. Let t(x0) denote the transition matrix with
0(x0) = In, where In is the n-dimensional identity matrix.
Define m1(t), . . . , mn(t) as the eigenvalues of t(x0). The i th Lyapunov exponentsλiof x0are then defined as
λi = lim
t→∞
1
t ln|mi(t)|, (3.2) whenever the limits exist for i = 1, . . . , n. Although there are as many Lyapunov exponents as there are dimensions, the largest nonzero Lyapunov exponent is usually the most im- portant one. If it is positive, the distance between neighbor- ing orbits grows exponentially over time. In other words, the system is sensitively dependent on initial conditions. If the largest Lyapunov exponent is negative, two nearby orbits will converge over time. We adopt the practical algorithm pre- sented by Parker and Chua17)to compute the Lyapunov expo- nents. The algorithm takes a great many iterations to compute the Lyapunov exponents with reasonable accuracy.
Figure 5 shows the largest nonzero Lyapunov exponents, λ1, plotted against mcfor mc∈ [8.00, 9.00]. When mclies in the intervals[8.00, 8.44], the corresponding λ1’s are negative.
These negative Lyapunov exponents indicate the presence of a periodic attractor in the corresponding systems. The Lya- punov exponents becomes positive when mc= 8.45 and this is the point that chaos starts. However, periodic windows do exist in the region[8.45, 9.00], where the corresponding Lya- punov exponents are negative. Comparing Figs. 4 and 5, we see the results are agree with each other nicely.
8 8.5 9
mc
2 103 0 2 103
Lyapunov exponent
Fig. 5. The largest nonzero Lyapunov exponent plotted against the dimen- sionless modulation current mc. The dimensionless modulation frequency mf= 3.73 × 10−3.
4. Synchronization of the Drive-Response System Consider a drive-response system using bistable laser diodes as chaotic masking devices, a drive system described by a three-dimensional rate equation is given by
ηd Ne1
d T = [Sp1+ mcsin(2π · mf· T ) + δ( ˆNp− Np)]
− α1Ne12 − (α2Ne12 + Ne1+ α3)Np− Ne1, (4.1)
ηd Ne2
d T = Sp2− α1Ne22 − (α2Ne22 + Ne2+ α3)Np− α4Ne2, (4.2) d Np
d T = Np[γ1(α2Ne12 + Ne1+ α3) + γ2(α2Ne22 + Ne2+ α3)]
− Np+ ε(γ1Ne12 + γ2Ne22).
(4.3)
The corresponding response system is given by ηd ˆNe1
d T = [Sp1+ mcsin(2π · mf· T ) + ˆδ(Np− ˆNp)]
− α1 ˆNe12 − (α2 ˆNe12 + ˆNe1+ α3) ˆNp− ˆNe1, (4.4)
ηd ˆNe2
d T = Sp2− α1 ˆNe22 − (α2 ˆNe22 + ˆNe2+ α3) ˆNp− α4 ˆNe2, (4.5) d ˆNp
d T = ˆNp[γ1(α2 ˆNe12 + ˆNe1+ α3) + γ2(α2 ˆNe22 + ˆNe2+ α3)]
− ˆNp+ ε(γ1 ˆNe12 + γ2 ˆNe22).
(4.6)
Note thatδ and ˆδ are the coupling coefficients.
The drive-response system described by (4.1)–(4.6) can express simplex (one-way) and duplex (two-way) transmis- sion. Figure 6 illustrates a possible simplex transmission schematic diagram, where δ = 0 and ˆδ = 0. To imple- ment a chaotic simplex transmission system, the drive sys- tem sends Np(t) to the response system. From the monitor photodetector of the identical bistable laser diode in the re- sponse system, ˆNp(t) is subtracted from Np(t). The coupling coefficient ˆδ corresponds to the tunable gain stage. By mixing with the signal current again, the total injection drive becomes Sp1+ mcsin(2π · mf· T ) + ˆδ(Np− ˆNp). A chaotic simplex
Signal Gain
Signal Sync.
LD
LD PD PD
(Drive) PD
(Response)
Fig. 6. An optical chaotic simplex transmission schematic diagram.
transmission thus can be implemented. Note that the coupling term is chosen to be ˆδ(Np− ˆNp), so the attractor of the coupled system belongs to the 45◦line in the Np− ˆNp space. On the other hand, ifδ = 0 and ˆδ = 0, the system becomes duplex transmission. In the case, the drive system sends a chaotic state variable Np(t) to the response system and, in return, the response system sends its chaotic state variable ˆNp(t) to the drive system.
4.1 Synchronization of simplex transmission
We verify the synchronization in simplex transmission by computing the solutions of the drive-response system and the conditional Lyapunov exponents. For all the following exper- iments, mc= 9.0 and mf= 3.73 × 10−3.
One way to verify the synchronization is to directly solve the drive-response system defined by (4.1) to (4.6). Figure 7 and 8 respectively show Npversus ˆNpfor a weak and a strong coupling for simplex transmission. In both of the weak and strong coupling cases, the initial values of Ne1, Ne2, and Np are all equal to 0.0 for the drive system; in contrast, the initial values of ˆNe1, ˆNe2, and ˆNpare all equal to 1.0 × 10−3for the response system. For the weak coupling where ˆδ = 1.0×10−4 andδ = 0, Fig. 7 shows trajectories of the drive and re- sponse systems quickly become uncorrelated, even both sys- tems have closely correlated initial conditions. In the case that the strong coupling is tuned to ˆδ = 1.0×10−1andδ = 0, syn- chronization occurs as shown in Fig. 8. Figure 9 further plots ˆNp/Npagainst ˆδ for δ = 0. The values of ˆNp/Npequal to one for ˆδ ≥ 4.2 × 10−3, which implies that the drive and response system are synchronized. In contrast, when ˆδ < 4.2 × 10−3, the values of ˆNp/Npvary and the two systems are not neces- sary synchronized.
Conditional Lyapunov exponents can also indicate syn- chronization. Pecora and Carroll11, 12) gave a necessary and sufficient condition for the synchronization: all the condi- tional Lyapunov exponents (CLE’s) associated with the vari- ational equations must be negative. To compute the CLE’s, we first let ENe1 = Ne1 − ˆNe1, ENe2 = Ne2 − ˆNe2 and
Fig. 7. ˆNp versus Np for ˆδ = 1.0 × 10−4 and δ = 0 in the simplex drive-response system.
=1 10
0 20 40 60
Np 0
20 40 60
Np
Fig. 8. ˆNp versus Np for ˆδ = 1.0 × 10−1 and δ = 0 in the simplex drive-response system.
0 0.005 0.01
0 4
Np/Np (log)
Fig. 9. The value of ˆNp/Np against ˆδ for δ = 0 in the simplex drive-response system.
EP = Np− ˆNp. Then construct the difference system
d E(t)
dt = AE(t). (4.7)
The real part of the three eigenvalues of A−1are the CLE’s by definition. If all the CLE’s are negative, the synchronization is achieved since limt→∞E(t) = 0. Otherwise, if positive CLE exists, the difference system grows apart as t goes to infinity. Figure 10 shows the first CLE’s of (4.7) against ˆδ for δ = 0, in which simplex transmission is achieved. All the second and third CLE’s are negative and are not shown in the figure. If the first CLE’s are positive, it implies the difference system grows apart as t → ∞. For ˆδ ≥ 4.2 × 10−3, all CLE’s are negative and synchronization are achieved.
4.2 Synchronization of duplex transmission
For duplex transmission where both ˆδ and δ are nonzero, synchronization can be established. If Npand ˆNpare directly solved from (4.1) to (4.6) for ˆδ = δ, ˆNp/Np = 1 implies that the drive and response systems are synchronized. Figure 11 presents ˆNp/Np against ˆδ for ˆδ = δ. Figure 12 shows the first CLE against ˆδ when ˆδ = δ. For ˆδ ≥ 2.2 × 10−3, no
0 0.005 0.01 3.5
0 3.5 7
The first CLE (103)
Fig. 10. The first conditional Lyapunov exponent of the drive-response system versus ˆδ for δ = 0 (simplex transmission).
0 0.005 0.01
3 0 3
Np/Np (log)
Fig. 11. The value of NpˆNp against ˆδ for ˆδ = δ in the duplex drive-response system.
0 0.005 0.01
0 1
The first CLE (103)
Fig. 12. The first conditional Lyapunov exponent of the drive-response system versus ˆδ for ˆδ = δ (duplex transmission).
positive CLE exists. The drive and response systems become synchronized in the sense that the difference between the two systems vanishes.
5. Conclusion
Bistable laser diodes can be used as chaotic masking de- vices to change sinusoidal electronic signals into optical chaotic light. Chaotic behaviors of the bistable laser diodes under an external drive were numerically examined. Period- doubling routes to chaos were observed in bifurcation dia- grams. Poincar´e maps and Lyapunov exponents were com- puted to investigate the chaotic phenomena solutions to the rate equations. The results are all consistent with each other nicely. The chaotic light also suggests a possible light source for secure optical communication. A chaotic state variable is transmitted between the drive and response systems in sim- plex or duplex transmission, by applying Pecora and Carroll’s theory. Synchronization can be established between the drive and response systems by applying strong coupling. The drive and response systems quickly become uncorrelated with weak coupling. The ability to design synchronizing systems using bistable laser diodes offers opportunities for applying chaos on optical information transmissions.
Acknowledgement
This work was partially supported by the National Science Council, R.O.C.
1) R. Reinisch and G. Vitrant: Prog. Quantum Electron. 18 (1994) 1.
2) H. Kawaguchi: IEEE J. Sel. Top. Quantum Electron. 3 (1997) 1254.
3) J. Chrostowski, R. Vallee and C. Delisle: Can. J. Phys. 61 (1983) 1143.
4) F. A. Narducci, G. L. Oppo, D. W. Bromley and J. R. Tredicce: Opt.
Commun. 75 (1989) 184.
5) J. Zemmouri, W. Sergent and B. S´egard: Opt. Commun. 84 (1991) 199.
6) H. Kawaguchi: Appl. Phys. Lett. 45 (1984) 1264.
7) H. Kawaguchi: Proc. IEE 129 (1982) 141.
8) H. Kawaguchi: Opt. Quantum Electron. 19 (1987) S1.
9) S. Bennett, C. M. Snowden and S. Iezekiel: IEEE J. Quantum Electron.
33 (1997) 2076.
10) E. Hemery, L. Chusseau and J. M. Lourtioz: IEEE J. Quantum Electron.
26 (1990) 633.
11) L. M. Pecora and T. L. Carroll: Phys. Rev. Lett. 64 (1990) 821.
12) L. M. Pecora and T. L. Carroll: IEEE Trans. Circuit Syst. 38 (1991) 453.
13) K. M. Cuomo, A. V. Oppenheim and S. H. Strogatz: IEEE Trans. Circuit Syst. II 40 (1993) 626.
14) C. Juang, T.-M. Hwang, J. Juang and W.-W. Lin: IEEE J. Quantum Electron. 36 (2000) 300.
15) J.-Y. Wang, M. Cada, R. Van Dommelen and T. Makino: IEEE J. Sel.
Top. Quantum Electron. 3 (1997) 1271.
16) K. Alligood, T. Sauer and J. A. Yorke: Chaos — an introduction to dynamical systems (Springer-Verlag, Heidelberg, 1996).
17) T. S. Parker and L. O. Chua: Practical Numerical Algorithms for Chaotic Systems (Springer-Verlag, Heidelberg, 1989).
