Notebooks · posts/chuas-circuit-double-scroll

The double scroll of Chua's circuit

Adam Zettel · September 25, 2026

Two capacitors, an inductor, a resistor and one diode whose current is a three-piece straight line. The voltages never settle and never repeat: they wind around two centers, hopping between them at irregular times. Nearby starts separate like $e^{λ t}$ with $λ = 0.33$ , which is what a separate scipy calculation gives within the scatter expected of a run this long.

x ˙ = α (y - x - g (x)), y ˙ = x - y + z, z ˙ = - β y

The exponent is measured here by following two copies of the circuit started $1 0^{- 8}$ apart and pulling them back to that distance after every step (Benettin's method). The reference, in posts/refs/chuas-circuit-double-scroll.py, uses a different method and a different integrator: scipy's DOP853 at a tolerance of $1 0^{- 11}$ , carrying the linearized equations along the orbit, over 24000 time units. It gives $λ = 0.3277 \pm 0.0026$ . Split into 40 windows of 600 time units (the length of this post's run), the estimates scatter with a standard deviation of 0.0165, so the test allows three of those: 0.05.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

The circuit

Chua's circuit is an electronic circuit designed by Leon Chua to be chaotic: its voltages follow a definite rule, with no randomness in it, yet never settle into a repeating pattern. It has two capacitors, an inductor, a resistor, and one nonlinear part, a Chua diode: a small circuit whose current $g (x)$ is a function of the voltage $x$ across it made of three straight pieces. In the middle the slope is $m_{0}$ and outside $∣ x ∣ > 1$ it is $m_{1}$ , and both slopes are negative, so the part feeds energy into the circuit instead of absorbing it.

In scaled units, with $x$ and $y$ the two capacitor voltages and $z$ the inductor current, the equations are the three above. The program uses the standard double-scroll values $α = 9$ , $β = 100/7$ , $m_{0} = - 8/7$ and $m_{1} = - 5/7$ .

[1]lines 39–68

ALPHA = 9.0
BETA = 14.285714285714286
M0 = -1.1428571428571428
M1 = -0.7142857142857143
# points on the diode's curve
N_DIODE = 81

#| The Chua diode's current at voltage x: slope M0 for |x| < 1, slope
#| M1 outside, joined continuously at x = ±1.
def diode(x) {
    return M1 * x + 0.5 * (M0 - M1) * (abs(x + 1.0) - abs(x - 1.0))
}

v = linspace(-3.0, 3.0, N_DIODE)
i_v = zeros(N_DIODE)
for k to N_DIODE {
    i_v[k] = diode(v[k])
}
show diode(0.5)
show diode(2.0)

figure("The Chua diode's current: three straight pieces", "voltage x", "current g(x)")
plot(v, i_v, "g(x)")
vline("x = -1", -1.0)
vline("x = 1", 1.0)

test "the diode is m0 x inside |x| < 1 and m1 x + (m0 - m1) outside" {
    assert_near(diode(0.5), -0.5714285714285714, 1.0e-15)
    assert_near(diode(2.0), -1.8571428571428572, 1.0e-15)
}

Inside $∣ x ∣ < 1$ the current is $m_{0} x$ , so $g (0.5) = - 4/7 = - 0.5714$ , and outside it is $m_{1} x + (m_{0} - m_{1})$ for $x > 1$ , so $g (2) = - 13/7 = - 1.8571$ . The negative slopes make the origin unstable: a small voltage grows instead of dying away. Far out, the circuit's other parts win and push it back. With nowhere stable to rest, the state keeps moving.

Stepping the equations

The program advances the state $(x, y, z)$ with the classical fourth-order Runge-Kutta method at a step of $0.005$ , then runs 20000 steps (100 time units) so the state reaches the set it keeps returning to, the attractor, before anything is measured.

[2]lines 85–124

DT = 0.005
SETTLE = 20000
# how long the settling run lasts, in time units
settle_time = SETTLE * DT
show settle_time

test "the settling run is 20000 steps of 0.005, 100 time units" {
    assert_near(settle_time, 100.0, 1.0e-10)
}

#| The right-hand side of Chua's equations at (x, y, z), written into out.
def rhs(x, y, z, out: Arr) {
    out[0] = ALPHA * (y - x - diode(x))
    out[1] = x - y + z
    out[2] = -BETA * y
    return 0.0
}

#| One classical Runge-Kutta step of length dt on the state st = [x, y, z];
#| k1 to k4 are scratch arrays for the four slopes.
def rk4(st: Arr, dt, k1: Arr, k2: Arr, k3: Arr, k4: Arr) {
    rhs(st[0], st[1], st[2], k1)
    rhs(st[0] + 0.5 * dt * k1[0], st[1] + 0.5 * dt * k1[1], st[2] + 0.5 * dt * k1[2], k2)
    rhs(st[0] + 0.5 * dt * k2[0], st[1] + 0.5 * dt * k2[1], st[2] + 0.5 * dt * k2[2], k3)
    rhs(st[0] + dt * k3[0], st[1] + dt * k3[1], st[2] + dt * k3[2], k4)
    for c to len(st) {
        st[c] = st[c] + dt / 6.0 * (k1[c] + 2.0 * k2[c] + 2.0 * k3[c] + k4[c])
    }
    return 0.0
}

sa = [0.7, 0.0, 0.0]
k1 = zeros(len(sa))
k2 = zeros(len(sa))
k3 = zeros(len(sa))
k4 = zeros(len(sa))
for s to SETTLE {
    rk4(sa, DT, k1, k2, k3, k4)
}
show sa

DT = 0.005
SETTLE = 20000
# how long the settling run lasts, in time units
settle_time = SETTLE * DT
show settle_time

test "the settling run is 20000 steps of 0.005, 100 time units" {
    assert_near(settle_time, 100.0, 1.0e-10)
}

#| The right-hand side of Chua's equations at (x, y, z), written into out.
def rhs(x, y, z, out: Arr) {
    out[0] = ALPHA * (y - x - diode(x))
    out[1] = x - y + z
    out[2] = -BETA * y
    return 0.0
}

#| One classical Runge-Kutta step of length dt on the state st = [x, y, z];
#| k1 to k4 are scratch arrays for the four slopes.
def rk4(st: Arr, dt, k1: Arr, k2: Arr, k3: Arr, k4: Arr) {
    rhs(st[0], st[1], st[2], k1)
    rhs(st[0] + 0.5 * dt * k1[0], st[1] + 0.5 * dt * k1[1], st[2] + 0.5 * dt * k1[2], k2)
    rhs(st[0] + 0.5 * dt * k2[0], st[1] + 0.5 * dt * k2[1], st[2] + 0.5 * dt * k2[2], k3)
    rhs(st[0] + dt * k3[0], st[1] + dt * k3[1], st[2] + dt * k3[2], k4)
    for c to len(st) {
        st[c] = st[c] + dt / 6.0 * (k1[c] + 2.0 * k2[c] + 2.0 * k3[c] + k4[c])
    }
    return 0.0
}

sa = [0.7, 0.0, 0.0]
k1 = zeros(len(sa))
k2 = zeros(len(sa))
k3 = zeros(len(sa))
k4 = zeros(len(sa))
for s to SETTLE {
    rk4(sa, DT, k1, k2, k3, k4)
}
show sa

The double scroll

Now follow the state for 210 time units (42000 steps). Record $x$ , and count how often it changes sign: the attractor has two lobes, one around $x > 0$ and one around $x < 0$ , and a sign change is a hop from one to the other. Before looking: if the motion were periodic, the hops would come at a fixed rhythm. Do they?

[3]lines 136–200

FILM = 42000
# keep every twentieth step (0.1 time units) for the pictures
EVERY = 20
# the x(t) figure shows the first 100 time units
SHOW_STEPS = 20000
# picture coordinates are rounded to 1e-3, well below a pixel, to keep the figure small
ROUND = 1000.0

#| x rounded to three decimals, for a figure.
def r3(x) {
    return round(x * ROUND) / ROUND
}

n_keep = FILM / EVERY
px = zeros(n_keep)
py = zeros(n_keep)
pz = zeros(n_keep)
n_time = SHOW_STEPS / EVERY
t_x = zeros(n_time)
x_t = zeros(n_time)
jumps = 0.0
prev_sign = 1.0
if sa[0] < 0.0 {
    prev_sign = -1.0
}
for s to FILM {
    rk4(sa, DT, k1, k2, k3, k4)
    sign = 1.0
    if sa[0] < 0.0 {
        sign = -1.0
    }
    if sign * prev_sign < 0.0 {
        jumps = jumps + 1.0
    }
    prev_sign = sign
    if mod(s, EVERY) == 0.0 {
        n = s / EVERY
        px[n] = r3(sa[0])
        py[n] = r3(sa[1])
        pz[n] = r3(sa[2])
        if n < n_time {
            t_x[n] = r3((s + 1.0) * DT)
            x_t[n] = px[n]
        }
    }
}
show jumps

figure("The voltage x over 100 time units", "time", "x")
plot(t_x, x_t, "x(t)")
hline("x = 0", 0.0)

figure("The double scroll, seen along z", "x", "y")
plot(px, py, "(x, y) over 210 time units")

plot3d(px, py, pz, "the double scroll in (x, y, z): drag to rotate")

# scipy's long run (posts/refs/chuas-circuit-double-scroll.out): 40.87 hops per 210 time units on average,
# with a standard deviation of 3.07 over 114 windows of that length
HOPS_SCIPY = 40.87
HOPS_SPREAD = 3.07

test "the orbit hops between the scrolls about 41 times in 210 time units, as scipy's long run does" {
    assert_near(jumps, HOPS_SCIPY, 3.0 * HOPS_SPREAD)
}

FILM = 42000
# keep every twentieth step (0.1 time units) for the pictures
EVERY = 20
# the x(t) figure shows the first 100 time units
SHOW_STEPS = 20000
# picture coordinates are rounded to 1e-3, well below a pixel, to keep the figure small
ROUND = 1000.0

#| x rounded to three decimals, for a figure.
def r3(x) {
    return round(x * ROUND) / ROUND
}

n_keep = FILM / EVERY
px = zeros(n_keep)
py = zeros(n_keep)
pz = zeros(n_keep)
n_time = SHOW_STEPS / EVERY
t_x = zeros(n_time)
x_t = zeros(n_time)
jumps = 0.0
prev_sign = 1.0
if sa[0] < 0.0 {
    prev_sign = -1.0
}
for s to FILM {
    rk4(sa, DT, k1, k2, k3, k4)
    sign = 1.0
    if sa[0] < 0.0 {
        sign = -1.0
    }
    if sign * prev_sign < 0.0 {
        jumps = jumps + 1.0
    }
    prev_sign = sign
    if mod(s, EVERY) == 0.0 {
        n = s / EVERY
        px[n] = r3(sa[0])
        py[n] = r3(sa[1])
        pz[n] = r3(sa[2])
        if n < n_time {
            t_x[n] = r3((s + 1.0) * DT)
            x_t[n] = px[n]
        }
    }
}
show jumps

figure("The voltage x over 100 time units", "time", "x")
plot(t_x, x_t, "x(t)")
hline("x = 0", 0.0)

figure("The double scroll, seen along z", "x", "y")
plot(px, py, "(x, y) over 210 time units")

plot3d(px, py, pz, "the double scroll in (x, y, z): drag to rotate")

# scipy's long run (posts/refs/chuas-circuit-double-scroll.out): 40.87 hops per 210 time units on average,
# with a standard deviation of 3.07 over 114 windows of that length
HOPS_SCIPY = 40.87
HOPS_SPREAD = 3.07

test "the orbit hops between the scrolls about 41 times in 210 time units, as scipy's long run does" {
    assert_near(jumps, HOPS_SCIPY, 3.0 * HOPS_SPREAD)
}

In 210 time units the state changed lobes 41 times (jumps): a few turns around one center, then a few around the other, with the number of turns changing from visit to visit. That irregular count is the visible sign of chaos. A scipy run a hundred times longer (posts/refs/chuas-circuit-double-scroll.py) averages 40.9 hops per 210 time units, and single stretches of that length range from 34 to 50, with a standard deviation of 3.1; the test allows three of those.

How fast nearby starts separate

Chaos has a measurable rate. Start two copies of the circuit a distance $d_{0} = 1 0^{- 8}$ apart. On average the distance grows like $d_{0} e^{λ t}$ ; $λ$ is the largest Lyapunov exponent, and a positive $λ$ is what chaos means in numbers. The two copies cannot be left alone for long, since once their distance is as large as the attractor it stops growing. So after every step the program measures the new distance $d$ , adds $ln (d / d_{0})$ to a running sum, and moves the second copy back along the same direction to distance $d_{0}$ . The exponent is the sum divided by the elapsed time.

The run is 120000 steps, 600 time units. The reference value from scipy is $0.3277$ , and a run of this length scatters about the true value with a standard deviation of $0.0165$ .

[4]lines 227–281

LYAP_STEPS = 120000
D0 = 1.0e-8
# points on the running-estimate figure
N_TRACE = 200
LAMBDA_SCIPY = 0.3277

sb = sa + 0.0
sb[0] = sa[0] + D0
trace_every = LYAP_STEPS / N_TRACE
trace_t = zeros(N_TRACE)
trace_lambda = zeros(N_TRACE)
log_sum = 0.0
for s to LYAP_STEPS {
    rk4(sa, DT, k1, k2, k3, k4)
    rk4(sb, DT, k1, k2, k3, k4)
    dx = sb[0] - sa[0]
    dy = sb[1] - sa[1]
    dz = sb[2] - sa[2]
    d = sqrt(dx * dx + dy * dy + dz * dz)
    log_sum = log_sum + log(d / D0)
    f = D0 / d
    sb[0] = sa[0] + dx * f
    sb[1] = sa[1] + dy * f
    sb[2] = sa[2] + dz * f
    if mod(s + 1.0, trace_every) == 0.0 {
        n = (s + 1.0) / trace_every - 1.0
        trace_t[n] = (s + 1.0) * DT
        trace_lambda[n] = log_sum / ((s + 1.0) * DT)
    }
}
lyap = log_sum / (LYAP_STEPS * DT)
show lyap

figure("The running estimate of the largest Lyapunov exponent", "time", "$\\lambda$ estimate", 0.0, 0.0, 0.0, 0.6)
plot(trace_t, trace_lambda, "ln(d / d0) summed, over elapsed time")
hline("scipy: 0.3277", LAMBDA_SCIPY)
# how far a 600-time-unit estimate scatters about scipy's value: one standard deviation
LAMBDA_SPREAD = 0.0165
fill_between("scipy's value, give or take one 600-unit run's spread", trace_t, zeros(N_TRACE) + (LAMBDA_SCIPY - LAMBDA_SPREAD), zeros(N_TRACE) + (LAMBDA_SCIPY + LAMBDA_SPREAD))

test "the largest Lyapunov exponent is scipy's 0.3277 within three standard deviations of a 600-time-unit estimate" {
    # posts/refs/chuas-circuit-double-scroll.out: 0.3277 +- 0.0026
    # over 24000 time units; one 600-unit window scatters by 0.0165
    assert_near(lyap, LAMBDA_SCIPY, 0.05)
}

horizon = log(1.0 / D0) / lyap
show horizon

test "a 1e-8 difference reaches order one in ln(1e8) / lambda = 56 time units, scipy's 56.2 within 15 percent" {
    # ln(1e8) / 0.3277 = 56.21; a 0.05 error in lambda moves it by about 15
    # percent, and the test allows the 3-sigma band of the measured lambda
    assert_near(1.0 / D0, 1.0e8, 0.0)
    assert_near(horizon, 18.420680743952367 / LAMBDA_SCIPY, 0.15 * 56.21)
}

LYAP_STEPS = 120000
D0 = 1.0e-8
# points on the running-estimate figure
N_TRACE = 200
LAMBDA_SCIPY = 0.3277

sb = sa + 0.0
sb[0] = sa[0] + D0
trace_every = LYAP_STEPS / N_TRACE
trace_t = zeros(N_TRACE)
trace_lambda = zeros(N_TRACE)
log_sum = 0.0
for s to LYAP_STEPS {
    rk4(sa, DT, k1, k2, k3, k4)
    rk4(sb, DT, k1, k2, k3, k4)
    dx = sb[0] - sa[0]
    dy = sb[1] - sa[1]
    dz = sb[2] - sa[2]
    d = sqrt(dx * dx + dy * dy + dz * dz)
    log_sum = log_sum + log(d / D0)
    f = D0 / d
    sb[0] = sa[0] + dx * f
    sb[1] = sa[1] + dy * f
    sb[2] = sa[2] + dz * f
    if mod(s + 1.0, trace_every) == 0.0 {
        n = (s + 1.0) / trace_every - 1.0
        trace_t[n] = (s + 1.0) * DT
        trace_lambda[n] = log_sum / ((s + 1.0) * DT)
    }
}
lyap = log_sum / (LYAP_STEPS * DT)
show lyap

figure("The running estimate of the largest Lyapunov exponent", "time", "$\\lambda$ estimate", 0.0, 0.0, 0.0, 0.6)
plot(trace_t, trace_lambda, "ln(d / d0) summed, over elapsed time")
hline("scipy: 0.3277", LAMBDA_SCIPY)
# how far a 600-time-unit estimate scatters about scipy's value: one standard deviation
LAMBDA_SPREAD = 0.0165
fill_between("scipy's value, give or take one 600-unit run's spread", trace_t, zeros(N_TRACE) + (LAMBDA_SCIPY - LAMBDA_SPREAD), zeros(N_TRACE) + (LAMBDA_SCIPY + LAMBDA_SPREAD))

test "the largest Lyapunov exponent is scipy's 0.3277 within three standard deviations of a 600-time-unit estimate" {
    # posts/refs/chuas-circuit-double-scroll.out: 0.3277 +- 0.0026
    # over 24000 time units; one 600-unit window scatters by 0.0165
    assert_near(lyap, LAMBDA_SCIPY, 0.05)
}

horizon = log(1.0 / D0) / lyap
show horizon

test "a 1e-8 difference reaches order one in ln(1e8) / lambda = 56 time units, scipy's 56.2 within 15 percent" {
    # ln(1e8) / 0.3277 = 56.21; a 0.05 error in lambda moves it by about 15
    # percent, and the test allows the 3-sigma band of the measured lambda
    assert_near(1.0 / D0, 1.0e8, 0.0)
    assert_near(horizon, 18.420680743952367 / LAMBDA_SCIPY, 0.15 * 56.21)
}

The estimate wanders early, while one or two passes near the center of a lobe can move it, and settles as the time grows. This run gives $λ = 0.327$ (lyap), against scipy's $0.3277$ . At that rate a difference of $1 0^{- 8}$ in the starting voltage grows to order one, the size of the attractor, in $ln (1 0^{8}) / λ \approx 56$ time units (horizon). Past that time the equations are still exact, but a measurement of the starting voltage to eight digits no longer predicts which lobe the circuit is on.

[5]lines 294–295

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors