CHAOS IN CHUA’S CIRCUIT
December 6, 2007
Math 53: Chaos!
INTRODUCTION AND BRIEF HISTORY
This project analyzes Chua’s circuit, a simple electronic circuit that demonstrates chaotic behavior. The study consists of two main components: numerically exploring (in Matlab) the equations derived from the circuit, and constructing the actual circuit to replicate numerical results with experimental data. The physical circuit demonstrates that one can create chaotic behavior in a relatively simple system.
Leon Chua invented the circuit in 1983 to demonstrate chaos in an actual physical model
1and to prove that the Lorenz double-scroll attractor is chaotic. The electronic circuit suits the
study of chaos well because one can precisely control its parameters and can readily observe the results on an oscilloscope. The circuit became popular to study because it is easy to construct, and many people have built the circuit using ordinary electronic components. In fact, one can model the circuit using only resistors, capacitors, inductors, op-amps, and diodes. In 1986, Chua,
2Komuro, and Matsumoto proved rigorously that the Chua attractor is indeed chaotic.
CROSS’S VERSION OF THE CIRCUIT AND ITS EQUATIONS
As mentioned above, people have created versions of the circuit that do not rely on any specially produced parts (the original included a laboratory-made ―Chua diode‖); rather, one can
build versions of the circuit using off-the-shelf components. Michael Cross proposed the version
3of the Chua circuit used in this project and shown in Figure 1. The right-hand side of this circuit
(all the components to the right of C1 in the diagram) simulates the Chua diode, providing nonlinear negative resistance. The operational amplifier and its resistors have a negative
1 http://www.scholarpedia.org/article/Chua_Circuit. 2 Ibid. 3 Cross. http://www.cmp.caltech.edu/~mcc/chaos_new/Chua.html.
4resistance of size –R1, and the diodes provide nonlinearity. The left side of the circuit acts as an
RLC circuit, which would simply produce damped oscillations without the right-hand side.
5Figure 1: Cross’s Version of Chua’s Circuit
To derive the three differential equations for the system, we choose three variables that change over time: V1, the voltage across capacitor C1; V2, the voltage across capacitor C2; and I, the current through the inductor. We apply Kirchhoff’s first law—which states that the current
entering a node equals the current leaving a node—to the nodes above C1 and C2 in the diagram.
Next we apply Kirchhoff’s second law—which states that the sum of the voltage around a loop equals zero—to the loop containing the inductor and C2. Thus, Kirchhoff’s laws give the following equations for the circuit:
C1(dV1/dt) = (V2-V1)/R – g(V1)
C2(dV2/dt) = -(V2-V1)/R + I
L(dI/dt) = -rI – V2
4 Cross. http://www.cmp.caltech.edu/~mcc/chaos_new/Chua_docs/works.html. 5 This diagram comes from Cross’s website. http://www.cmp.caltech.edu/~mcc/chaos_new/Chua.html.
Note that g(V1) is the current through the nonlinear part of the circuit. Michael Cross transforms
6the above equations into the following nondimensional system of equations:
dX/dt = a(Y-X) – G(X)
dY/dt = s[-a(Y-X) + Z]
dZ/dt = -c(Y + pZ)
a = R1/R b = 1-R1/R2 c = (C1*R1^2)/L s = C1/C2 p = r/R1
7and G(X) is the normalized current-voltage characteristic for the right part of the circuit:
G(X) = -X, |X|?1
The current-voltage characteristic for the negative nonlinear resistance in normalized terms is:
To find the equilibrium points, we set the three differential equations equal to zero and let G(X)=+(1-b)-bX to derive (0,0,0) and:
X = + (1-b)/[a/(1+pa) – b] Y = [ap/(1+ap)]*X Z = -Y/p
These fixed points will be used to plot the bifurcation diagram and a 2D attractor.
6 Cross. http://www.cmp.caltech.edu/~mcc/chaos_new/Chua_docs/chua_eq.pdf. 7 Cloyd 9. Cloyd noticed sign errors in the G(X) function published by Cross. G(X) here includes Cloyd’s
NUMERICAL ANALYSES IN MATLAB
Matlab’s ode45 can numerically model the three normalized differential equations to produce attractors for the circuit. For the numerical analyses, we will use the values for R, R1, R2, C2, L, and r that are listed in Figure 1. We will vary C1 to change the behavior of the system. The default initial condition will be (0.1, 0.15, 0.05). Figure 2 shows a periodic orbit at
8C1=4.28e-9 and a chaotic orbit at C1=4.32e-9. In order to accurately model the attractors, I had
to adjust the ode45 settings of relative and absolute tolerance to 1e-6 and 1e-9, respectively.
Figure 2: Orbits when C1=4.28e-9 and C1=4.32e-9
To better understand the behavior of the circuit as C1 varies, we can examine a bifurcation diagram. To create the bifurcation diagram, I considered a projection of the orbit onto the X-Y plane (which corresponds to the normalized V1-V2 plane) and where the orbits
8 Cloyd identified these values for C1 as nonchaotic and chaotic; a bifurcation diagram also justifies these values.
intersect the projection onto the X-Y plane of the line which passes through the equilibrium
9points. Figure 3 shows a visualization of these intersection points. To find the points in Matlab, I used a flag variable that triggered if a given Y-value of the orbit was below the line between equilibrium points and the next Y-value was above it, or vice versa. The program then calculated the intersection of two lines projected onto the X-Y plane: the line between the equilibrium points and the line between the two consecutive points in the orbit that crossed the equilibrium-point line.
Figure 3: Procedure to Create Bifurcation Diagram
The bifurcation diagram used the X (normalized V1) values of the intersections and plotted them along the vertical axis. The varying values of C1 were plotted along the horizontal axis. To examine a wide range of behavior—including a period doubling route to chaos—I let
C1 vary between 4.2e-9 and 4.5e-9. In order to balance resolution with computing time, I used 1000 intervals. Figure 4 represents the end-result of an overnight computation: a clear bifurcation diagram for the values of C1. The diagram illustrates the periodic orbit at C1=4.28e-
9 The idea to use this method came from Cloyd’s paper; the Matlab code was my own.
9 and the chaotic orbit at C1=4.32e-9. Also, as C1 decreases from 4.5 to 4.4 nano-farads, we witness the period-doubling route to chaos.
Figure 4: Bifurcation Diagram for C1=4.2 to C1=4.5
We can also analyze a 2-dimensional surface of section for a chaotic orbit. For the SOS plane, I chose to extend the equilibrium line (used in creating the bifurcation diagram) in the Z dimension. This defines a plane whose normal is the cross product of (x1, y1, 0) and (0, 0, 1), where (x1, y1, z1) is a nontrivial equilibrium. Because the Z value should not affect the intersections, my code detects when the orbit crosses the plane in the same manner it calculated when the orbit crossed the equilibrium line in the X-Y plane. The code then calculates the
10intersection of the plane with the line between the two orbit points (on either side of the plane).
To plot a graph of the intersection points on the surface-of-section plane, I used the Z values of
10 Bourke. http://local.wasp.uwa.edu.au/~pbourke/geometry/planeline/. Bourke provides a nice formula for finding plane-line intersections using vectors. I used this method in my Matlab code.
the points as the horizontal axis. The vertical axis represents the distance of the intersection points from the origin in the X-Y plane, which is defined to be positive when X?0 and negative
when X<0. This works because the points, when projected onto the X-Y plane, lie on a straight line that passes through the origin. I analyzed the two-dimension attractor formed by the intersection of this surface-of-section with the chaotic orbit produced when C1=4.32e-9. For t=[0..1 million], the code found around 43 thousand points of intersection. Figure 5 shows the intersection points plotted on the surface of section plane. The Boxcount program calculated the
11boxcounting dimension of the 2D attractor to be about 0.9. The code also calculated the
12correlation dimension to be around 0.84.
Figure 5: 2D Surface of Section for C1=4.32e-9
A hallmark of chaos is sensitive dependence on initial conditions. In order to test this, I launched two trajectories with ICs very close together and plotted the distance between the
11 F. Moisy wrote a Matlab program to calculate the boxcounting dimension. http://www.mathworks.com/ matlabcentral/files/13063/content/boxcount/html/demo.html. 12 Professor Barnett wrote the code to find the correlation dimension.
points over time. Since ode45 does not allow the user to specify a time interval, I coded a loop that solved each initial condition problem for t=[0..1] then recorded the points at t=1 and used them as new initial conditions. My code then outputs a graph of the natural log of the Euclidean
distance between the orbits versus time. Figure 6 shows the graphs for C1=4.28e-9, which is a periodic orbit, and for C1=4.32e-9, which is a chaotic orbit. The initial conditions for the orbits were (0.1, 0.15, 0.05) and (0.1, 0.15, .05+1e-10), and the code ran 3000 iterations of the ode45 loops to show the distance between the orbits from t=0 to t=3000. For the periodic attractor, the two orbits remain close together; for the chaotic attractor, the two orbits grow exponentially far apart. To estimate a Lyapunov exponent, I exported the data for the chaotic attractor into Excel and found the slope of the line to be 0.011. The Chaos textbook reports the largest exponent of a
13chaotic Chua orbit to be 0.48. However, the book uses different equations to model the circuit, and we must keep in mind that we have normalized the time variable. In any case, we can clearly see that when the attractor is chaotic, the solutions are extremely sensitive to initial conditions.
Figure 6: Distance Between Orbits with Close ICs for C1=4.28e-9 (blue) and C1=4.32e-9
13 Alligood, et al. 390.
EXPERIMENTAL MODEL OF THE CIRCUIT
In order to collect experimental data on the Chua attractor, I built a working model of the circuit according to the Cross diagram in Figure 1. Due to the availability of parts, I had to change a few parameters: L=10mH, r=0.110K, and C2=68nF. For R, I used three potentiometers in series, whose maximum values were 1000, 500, and 25 ohms. For C1, I used a capacitance decade box (which had 0.1nF steps) in parallel with a 100pF tuning capacitor. To compensate for imperfect voltage sources, a small capacitor was placed across the V- and V+
14pins of the op-amp. Since the circuit operates at audible frequencies, it can output to a computer’s line-in jack. In order to do this, I used resistors to step down the voltage to an acceptable level for the input. The two outputs and the ground were connected to a 1/8‖ male
stereo connector. To prevent the output from altering the dynamics of the circuit, I used a pair of
15op-amps to create a buffer for each line out. Figure 7 shows the circuit diagram for each output
line, and Appendix A contains digital pictures of my circuit.
Figure 7: Output Buffer and Voltage Step-Down
The circuit was very sensitive to changes in R and C1, whose values I measured using a multimeter and an RCL-meter. Through trial and error, I found interesting behavior by setting
14 Professor Charles Sullivan proposed this to improve the circuit. I should also note that Professor Sullivan also helped me to debug my initial circuit. 15 Professor Barnett suggested this modification.