The Duffing oscillator is a well-known and widely studied system in nonlinear dynamics, often used to illustrate complex behaviors such as bifurcations and chaos. It serves as a quintessential example of how simple nonlinear systems can produce rich and sometimes unpredictable dynamics. Typically, it models a damped oscillator with a nonlinear restoring force and an external periodic driving force. Although the specific notation and parameter names can vary slightly across different texts and studies, the underlying structure of the equation remains fundamentally the same. The system's sensitivity to initial conditions and parameter values makes it a powerful demonstration tool for exploring the transition from regular motion to chaotic regimes in physics and engineering.

• Build Up Intuition
• How it works ?

Build Up Intuition

Before diving into theory, I would like you to have a chance to play with this simple simulation and build your own intuition about the dynamics of the Duffing oscillator. Building up intuition is important because it grounds abstract mathematical concepts in visual and experiential understanding. Especially in nonlinear systems like this, where small parameter changes can lead to dramatically different outcomes, intuition helps you recognize patterns, anticipate behaviors, and grasp the qualitative nature of stability, bifurcation, and chaos—long before equations reveal them formally. By interacting with the system directly, you train your mind to connect cause and effect in dynamic environments, which is essential not just for learning, but for creativity and problem-solving in complex systems.

Following is usage of this program.

Starting the Visualization

Both simulators run inside this page, so you do not need to install anything. The first bullet below comes from the standalone version of the tool, and here the simulator is already embedded above. It integrates the equation with a fourth-order Runge-Kutta step, and nothing moves until you press Start.

• Open DuffingOscillator.html in a modern web browser
• Click the "Start" button to begin the simulation
• Click "Pause" to temporarily stop the simulation
• Use "Reset" to return to initial conditions

Parameter Controls

The sliders map one to one onto the equation x'' + δx' + αx + βx³ = γcos(ωt), which is the form the simulator integrates. The sign of α matters most. A positive α gives a single potential well, while a negative α with a positive β gives two wells. Four of the seven presets use α = -1, so they start from the double-well shape.

• α (Alpha): Controls the linear stiffness term
• β (Beta): Controls the nonlinear stiffness term
• δ (Delta): Controls the damping coefficient
• γ (Gamma): Controls the amplitude of periodic forcing
• ω (Omega): Controls the frequency of periodic forcing
• Initial Conditions: Set x₀ and ẋ₀ for different starting points

Visualization Options

The two plots show the same trajectory in two different ways. The time series keeps time on the horizontal axis, so it shows when the motion changes. The state space plot drops time and draws velocity against displacement. A periodic motion then closes on itself, and a chaotic motion keeps filling a band.

• Time Series Plot: Shows x and ẋ vs time
• State Space Plot: Displays the phase portrait (x vs ẋ)
• Trail Length: Adjust how many points are shown in the plots (100 to 100,000)
• Color Scheme: Choose from:
  • Rainbow
  • Fire
  • Ocean
  • Grayscale

Preset Configurations

Each preset loads a full parameter set, and its name tells you the behaviour to look for. Let the transient die away first, and then watch the state space plot. The Default preset settles onto a closed periodic orbit. The Chaotic Attractor preset keeps drawing new loops that never repeat.

• Default: Standard parameter set
• Chaotic Attractor: Shows chaotic behavior
• Period Doubling: Demonstrates period-doubling bifurcation
• Soft Spring: Models a softening spring system
• Hard Spring: Models a hardening spring system
• Escape from Well: Shows escape from potential well
• Period-3 Orbit: Displays a period-3 solution

Zoom Controls

A chaotic trajectory spreads over a wide range, but the detail you want is often a small fold near one of the wells. Zoom lets you enlarge one region of a plot without restarting the simulation.

Both Time Series and State Space plots have independent zoom options:

• Both Axes: Zoom affects both x and y axes
• X-Axis Only: Horizontal zoom only
• Y-Axis Only: Vertical zoom only

Additional Features

These extras support the experiments above, and they do not change the physics. The time step is the one to watch. A step that is too large lets the Runge-Kutta error grow, and in a chaotic regime that error changes the trajectory you see.

• Mouse coordinates are displayed in the bottom-right corner
• Save the current view as an image using the "Save Image" button
• Adjust the time step for simulation accuracy (0.001 to 0.05)
• Real-time parameter adjustment
• Smooth animation with configurable trail length

There is another way of visuallizing the chaotic behavior of a system. Poincare map is the one. A Poincaré map offers a powerful way to visualize the complex behavior of dynamical systems, especially in the context of chaos. Unlike traditional time series plots that show how a variable evolves over time, or state space trajectories that illustrate continuous paths through a system's phase space, a Poincaré map samples the system at discrete intervals—typically once per forcing period in a driven oscillator. This results in a scatterplot that condenses high-dimensional behavior into a two-dimensional slice, making it easier to identify periodicity, bifurcations, or chaotic regimes. While time series and state space plots help us see the continuous evolution of motion, the Poincaré map captures the underlying structure of that motion, offering a more distilled view of long-term dynamics and stability.

Followings are the description on how you play with this program

Parameters

The Poincare map tool integrates the same equation as the simulator above, again with a fourth-order Runge-Kutta step. Its default set is α = -1, β = 1, δ = 0.2, γ = 0.37 and ω = 1.2, which is a driven double-well oscillator. A numerical check gives a positive largest Lyapunov exponent, about 0.15 per unit time, for this set and for each of the other three presets. So all four presets produce a strange attractor on the map.

The Duffing oscillator is defined by the equation:

ẍ + δẋ + αx + βx³ = γcos(ωt)

Adjustable parameters:

• α (alpha): Stiffness coefficient [-2 to 2]
• β (beta): Nonlinear stiffness [-2 to 10]
• δ (delta): Damping coefficient [0 to 1]
• γ (gamma): Forcing amplitude [0 to 10]
• ω (omega): Forcing frequency [0.1 to 2]
• x₀: Initial position
• ẋ₀: Initial velocity

Simulation Controls

The integration settings decide how much of the attractor the map shows. The tool records one point at the end of every forcing period T = 2π/ω. With the default 10000 total periods and 500 transient periods, the map therefore carries 9500 points. Steps/Period sets the Runge-Kutta step as T divided by that number.

• Preset Settings: Choose from predefined configurations:
  • Default Chaotic
  • Period Doubling
  • Double Periodic
  • Single Scroll
• Integration Parameters:
  • Steps/Period: Number of integration steps per forcing period
  • Total Periods: Total number of periods to simulate
  • Transient: Number of initial periods to discard
• Display Options:
  • WebGL: Toggle WebGL acceleration
  • Grid: Toggle grid display

Interactive Features

Zoom and pan matter more here than in the time plots, because a chaotic Poincare map has fine structure. When you enlarge one fold of the attractor, you find thinner folds inside it. This layered structure is the fractal signature of a strange attractor.

• Mouse Control:
  • Zoom: Mouse wheel
  • Pan: Click and drag
  • Coordinates: Real-time position display
• Simulation Control:
  • Start/Pause: Control simulation
  • Reset: Return to initial conditions
  • Save: Export current view as image

How it works ?

The Duffing oscillator is a well-known example of a system used to demonstrate chaotic behavior and dynamics. It is represented mathematically by a second-order differential equation that incorporates terms for displacement, velocity, and nonlinear restoring forces. The general form of the equation is:

The meaning of each terms and coefficients are as follows :

x = displacement

dx/dt = velocity

d²x/dt² = acceleration

δ = damping coefficient

ω₀ = natural frequency of oscillation

β = nonlinearity coefficient

γ = forcing amplitude

ω = forcing frequency

φ = phase shift of forcing function

Here, each term and coefficient has a specific meaning. The variable x represents displacement, dx/dt represents velocity, and d²x/dt² represents acceleration. The parameter δ is the damping coefficient, which quantifies the system's resistance to motion. The terms β and ω₀ describe the nonlinear and natural frequency of oscillation, respectively. The forcing term γ cos(ωt + φ) accounts for an external periodic force, where γ is the forcing amplitude, ω is the forcing frequency, and φ is the phase shift of the forcing function.

Let's connect this form with the one the two simulators use. The coefficient of x in the equation above is ω₀², and the simulators call the same coefficient α. So the two forms are identical when α = ω₀² and φ = 0, and β keeps its name in both. The simulators also accept a negative α. In that case ω₀ is no longer a real frequency, and the potential changes shape, as the next paragraph shows.

A potential energy view explains why the sign of α matters so much. Without damping and forcing, the restoring force -αx - βx³ comes from the potential V(x) = αx²/2 + βx⁴/4. With α > 0 and β > 0, V has a single minimum at x = 0. With α < 0 and β > 0, the origin becomes a hill, and two wells appear at x = +/-√(-α/β). For α = -1 and β = 1, the wells sit at x = +/-1 and the hill between them is 0.25 high. Chaotic motion in this double-well case usually means that the trajectory jumps between the two wells in an irregular order.

The cubic term also makes the oscillation frequency depend on the amplitude. For a small β, a free undamped oscillation of amplitude A has a frequency close to ω₀ + 3βA²/(8ω₀). A positive β raises the frequency as A grows, and this is called a hardening spring. A negative β lowers it, and this is called a softening spring. For ω₀ = 1, β = 0.1 and A = 1, the formula gives 1.0375, and a direct numerical integration gives 1.0367.

The characteristics of the Duffing oscillator arise from the nonlinearity introduced by the x³ term. This nonlinearity enables the system to exhibit chaotic behavior under certain conditions. Depending on the parameter ranges, the oscillator may display bounded, periodic oscillations or chaotic dynamics. Varying the forcing frequency ω relative to the natural frequency ω₀ has a significant impact on the behavior of the system. As parameters are adjusted, the oscillations may switch between regular periodic motion and chaotic trajectories.

An important feature of the Duffing oscillator is its sensitivity to initial conditions. Even slight changes in the starting values can result in vastly different trajectories over time. The phase shift φ of the external forcing function plays a crucial role in determining when chaotic regimes emerge.

One point about φ needs care. If you replace t by t - φ/ω, cos(ωt + φ) becomes cos(ωt), so φ only moves the time origin of the forcing. It does not create or remove attractors. What φ changes is which attractor a given starting point x(0), dx/dt(0) reaches, and at which phase of the cycle the Poincare samples are taken.

To explore the behavior of the Duffing oscillator, tools like phase portraits and Poincaré sections can be used. These visualizations reveal the underlying attractors of the system and help illustrate its transition between order and chaos. By numerically solving the equation and experimenting with different parameters and initial conditions, one can observe how the system evolves and how its dynamics change in response to external influences.

The characteristics of this equation can be summarized as follows :

• The x³ term introduces nonlinearity which enables chaotic dynamics.
• The system exhibits chaotic, bounded oscillations for certain parameter ranges.
• Varying the forcing frequency ω relative to the natural frequency ω₀ impacts the behavior.
• The oscillations switch between periodic and chaotic as parameters change.
• Slight changes in initial conditions yield dramatically different trajectories over time.
• The phase shift φ of the external forcing alters when chaotic regimes emerge.
• Tools like phase portraits and Poincare sections reveal the underlying attractor.

You can plot out the solution of this equation using a simple script below. Play with parameters and initial conditions and see how the plot changes.

Let's read the plot above against the Matlab script that produced it. The script uses δ = 0.06, ω₀ = 1, β = 1, γ = 6, ω = 1 and φ = 0, and it starts from y(1) = 3.0 and y(2) = 4.1. Here y(1) is x and y(2) is dx/dt. The time series shows no repeating pattern up to t = 100. The phase portrait on the right keeps adding loops on both sides of y(1) = 0, even though this equation has a single potential well.

Be careful about how long that irregular motion lasts. If you run the same Matlab case to t = 400, the motion settles after roughly 40 forcing periods, near t = 240, onto a periodic orbit that repeats every three forcing periods. So the plot above shows transient chaos, not a chaotic attractor. Only the long-run behaviour tells you which of the two you have.

Also watch the sign of the linear term when you run the Python script. The Matlab line computes -(beta*y(1)^3 + w0^2*y(1)), which matches the equation above. The Python line computes - beta*y1**3 + w0**2*y1, so its linear term has the opposite sign. The Python script therefore solves x'' + δx' - ω₀²x + βx³ = γcos(ωt + φ), which is the double-well case with α = -1. That run stays chaotic, and its plot differs from the one above. To reproduce the Matlab case, change + w0**2*y1 to - w0**2*y1. Even then, expect the same general shape and not the same curve, because odeint and ode45 use different steps and tolerances.

• The sign of α sets the shape of the potential : a positive α gives one well, and a negative α with a positive β gives two wells at x = +/-√(-α/β).
• The cubic term makes the frequency depend on the amplitude : a positive β hardens the spring, and a negative β softens it.
• The phase φ only shifts the time origin : it changes which attractor a starting point reaches, but not which attractors exist.
• A Poincare map samples once per forcing period : a periodic orbit shows as a few points, and a chaotic one shows as a folded, layered set.
• A short run can look chaotic and still settle : the Matlab case above becomes periodic after about 40 forcing periods.
• Check the equation before you trust a plot : the Python listing flips the sign of the linear term, so it simulates a different oscillator from the Matlab listing.