Duffing Oscillator

Definition

The Duffing oscillator is a driven, damped oscillator whose restoring force is not a straight line but a cubic. Its equation of motion is:

x'' + δ·x' + α·x + β·x³ = γ·cos(ω·t)

Read it term by term:

• x'' — acceleration (inertia).
• δ·x' — linear viscous damping; δ > 0 drains energy.
• α·x — the ordinary linear (Hooke's-law) spring stiffness.
• β·x³ — the cubic restoring term. This single nonlinear addition is what makes the oscillator a Duffing oscillator rather than a plain harmonic one.
• γ·cos(ω·t) — an external periodic drive of amplitude γ and frequency ω.

Set β = 0 and you have the familiar linear driven oscillator. Keep β and you unlock three behaviours the linear system can never produce: a bent resonance curve, hysteretic amplitude jumps, and — for a strong enough drive — deterministic chaos. Named after the German engineer Georg Duffing, who studied it in 1918 while analysing forced vibrations in machinery, it has become the canonical textbook example of a chaotic mechanical system.

How the cubic term changes everything

The potential energy whose gradient gives the restoring force is:

V(x) = ½·α·x² + ¼·β·x⁴      →      F = -dV/dx = -(α·x + β·x³)

The shape of V(x) depends on the signs of α and β, and that shape dictates the entire personality of the system:

• α > 0, β > 0 (single well, hardening): one parabola-like minimum, but the walls steepen faster than a parabola. The effective stiffness rises with amplitude, so the resonance frequency climbs as you push harder.
• α > 0, β < 0 (single well, softening): the walls flatten at large amplitude; resonance frequency drops as you push harder.
• α < 0, β > 0 (double well): the V(x) = -½x² + ¼x⁴ shape has two minima at x = ±√(−α/β) separated by a hump at the origin. This is the buckled-beam potential and the easiest road to chaos.

Because the restoring force is amplitude-dependent, the resonance peak bends. In a linear oscillator the response amplitude versus drive frequency is a single symmetric peak at ω = √α. In the Duffing oscillator that peak leans over until, in a band of frequencies, it folds back on itself — and a folded curve means more than one amplitude is possible at the same drive frequency. That fold is the seed of bistability, jumps, and ultimately chaos.

A worked example: the hysteretic jump

Take a hardening single-well oscillator with δ = 0.1, α = 1, β = 0.04, γ = 1.5, and slowly sweep the drive frequency ω. Using the harmonic-balance approximation, the steady-state amplitude A satisfies the cubic-in-A² frequency-response relation:

[ (α - ω² + ¾·β·A²)² + (δ·ω)² ] · A² = γ²

Now read off the numbers as ω rises from low to high:

• At ω = 0.9 the relation has one solution: A ≈ 1.6. The system sits on the low branch.
• By ω ≈ 1.25 the bent peak has formed and three solutions exist — say A ≈ 4.1 (high), A ≈ 1.0 (unstable middle), A ≈ 0.7 (low). The system, having ridden up the low branch, is now near the unstable region.
• At the upper fold ω ≈ 1.32 the low branch you were riding simply disappears (a saddle-node bifurcation). The amplitude has nowhere to go but the high branch — it jumps up by roughly a factor of three in a single drive cycle.
• Continue up and the high branch falls gently.

Sweep ω back down and the story is not the reverse: you ride the high branch past 1.32, all the way to a lower fold ω ≈ 1.05, where the high branch vanishes and you jump down. The up-jump (1.32) and down-jump (1.05) happen at different frequencies — the response depends on which way you came. That memory of the past is hysteresis, and the band 1.05 < ω < 1.32 is the bistable region where two stable amplitudes genuinely coexist.

From order to chaos

Switch to the double well: δ = 0.3, α = −1, β = 1, ω = 1.2, and dial the drive amplitude γ upward. This is the most-studied chaotic parameter family in the literature.

Drive γ	Behaviour	Poincaré section
0.20	Settles into one well, small periodic wobble	1 point
0.28	Period-1 limit cycle (one swing per drive period)	1 point
0.29	Period-2 — first period-doubling bifurcation	2 points
0.37	Chaotic — irregular hopping between both wells	Fractal cloud (strange attractor)
0.50	Robust chaos; sensitive dependence on initial conditions	Dense fractal scatter
0.65	Periodic window re-emerges (large-amplitude limit cycle)	Few points

The route to chaos is a period-doubling cascade: as γ grows, the period-1 cycle becomes period-2, then period-4, period-8, and so on, with the bifurcations crowding closer together (governed by the universal Feigenbaum constant δ_F ≈ 4.669) until the period becomes effectively infinite — that is chaos. In the chaotic regime the ball hops irregularly between the two wells, and two trajectories that start a millionth apart diverge exponentially: the largest Lyapunov exponent becomes positive. The motion is fully deterministic — the ODE has no random term — yet long-term prediction is impossible. That is the essence of the strange attractor the Poincaré section reveals.

Variants and regimes

• Undriven, undamped (γ = 0, δ = 0): a conservative anharmonic oscillator. Energy is conserved; orbits are closed; the period depends on amplitude (unlike a linear spring).
• Hardening spring (β > 0): resonance bends to higher frequency; the jump-up happens above the linear resonance. Common in stretched membranes and large-deflection beams.
• Softening spring (β < 0): resonance bends to lower frequency; the curve overhangs to the left. Seen in pendulums (sin θ ≈ θ − θ³/6 gives an effective β < 0) and in some MEMS devices.
• Double well (α < 0): two stable equilibria; supports inter-well hopping, fractal basin boundaries, and chaos at modest drive.
• Holmes / van der Pol–Duffing hybrids: add a velocity-dependent nonlinearity for self-sustained oscillations on top of the cubic stiffness.

Common pitfalls and misconceptions

• "The βx³ term is a small correction." It is small only at small amplitude. Near resonance the response grows until the cubic term dominates — that is exactly when the curve bends and the jump appears. Treating it perturbatively can mislead you about the global picture.
• "Chaos means the math has a random term." No. The Duffing equation is perfectly deterministic. Chaos is sensitivity to initial conditions, not randomness; the same initial state always gives the same trajectory.
• "Bigger drive always means bigger steady amplitude." In the bistable band, sweeping direction matters more than drive level — you can be on the high or low branch at the very same γ and ω.
• "You can solve it like a harmonic oscillator." There is no closed-form solution. Use harmonic balance / multiple scales for the response curve and numerical integration (RK4) for the chaotic regime.
• "Numerical chaos is just round-off error." A genuine positive Lyapunov exponent is reproducible across step sizes and integrators; round-off changes the specific trajectory but not the statistical structure of the attractor.
• "The Poincaré section is the trajectory." It is a stroboscopic snapshot — one sample per drive period. A periodic orbit shows finitely many points; chaos shows a fractal set.

Applications

• Buckled beams and structures. A column compressed past its buckling load has two bent equilibria; vibrating its base realises the double-well Duffing oscillator. Engineers use it to study snap-through and fatigue.
• MEMS and NEMS resonators. Micro-cantilevers driven hard enter the nonlinear regime; the bent resonance and jump set the usable dynamic range of sensors and oscillators.
• Energy harvesting. Bistable Duffing harvesters use inter-well hops to capture broadband ambient vibration far more efficiently than a tuned linear device.
• Electronics. An LC circuit with a nonlinear capacitor or a Josephson junction obeys a Duffing-type equation; the jump phenomenon underlies bifurcation amplifiers used to read out superconducting qubits.
• Ship and offshore dynamics. Roll motion with a softening righting moment is modelled as a Duffing oscillator to predict capsize.
• Pedagogy. It is the simplest mechanical system that exhibits the full period-doubling route to chaos, so it anchors the teaching of nonlinear dynamics alongside the logistic map.

Derivation and numerical analysis

To find the bent response curve, assume a steady single-harmonic response x(t) ≈ A·cos(ωt − φ) and substitute into the equation. The cubic term gives cos³ = ¾cos + ¼cos(3·) ; harmonic balance keeps only the fundamental, yielding the amplitude relation quoted above:

[ (α - ω² + ¾·β·A²)² + (δ·ω)² ] · A² = γ²

This is a cubic in A². Where it has three real positive roots, the system is bistable; the fold (jump) frequencies are where two roots merge — i.e. where the discriminant vanishes (a saddle-node bifurcation). To go beyond steady state, integrate numerically. Recast as a first-order system and step with RK4:

// Duffing oscillator: x'' + δ x' + α x + β x³ = γ cos(ω t)
// State y = [x, v];  v = x'
function duffing(t, y, p) {
  const [x, v] = y;
  const dxdt = v;
  const dvdt = -p.delta * v - p.alpha * x - p.beta * x * x * x
               + p.gamma * Math.cos(p.omega * t);
  return [dxdt, dvdt];
}

function rk4Step(f, t, y, h, p) {
  const k1 = f(t, y, p);
  const k2 = f(t + h / 2, [y[0] + h / 2 * k1[0], y[1] + h / 2 * k1[1]], p);
  const k3 = f(t + h / 2, [y[0] + h / 2 * k2[0], y[1] + h / 2 * k2[1]], p);
  const k4 = f(t + h,     [y[0] + h * k3[0],     y[1] + h * k3[1]],     p);
  return [
    y[0] + h / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]),
    y[1] + h / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]),
  ];
}

// Classic chaotic double-well set
const p = { delta: 0.3, alpha: -1, beta: 1, gamma: 0.5, omega: 1.2 };
const T = 2 * Math.PI / p.omega;   // drive period
const h = T / 400;                 // 400 steps per drive period

let t = 0, y = [1, 0];
const poincare = [];
for (let n = 0; n < 400 * 4000; n++) {
  y = rk4Step(duffing, t, y, h, p);
  t += h;
  // Sample once per drive period -> Poincaré section
  if (n % 400 === 0 && t > 200) poincare.push([y[0], y[1]]);
}
// `poincare` traces a fractal strange attractor for gamma = 0.5
console.log(poincare.length, 'Poincaré points');

Cost and accuracy. RK4 is O(h⁴) accurate per step and O(h⁵) in local truncation error; for chaos you typically want 200–500 steps per drive period so the fastest motion is well resolved. Building a clean Poincaré section means discarding the initial transient (here t > 200, a few hundred drive periods) and then sampling once per period — so a 4000-period run at 400 steps each is 1.6 million RK4 evaluations, milliseconds in JavaScript. A bifurcation diagram repeats this across, say, 600 values of γ; estimating the largest Lyapunov exponent adds a second nearby trajectory whose separation you renormalise each period. None of it is expensive, which is precisely why the Duffing oscillator is the favourite numerical playground for first encounters with chaos.