Integrate harmonic, damped, Van der Pol, Duffing, and Lorenz systems using RK4. Explore attractors, limit cycles, and chaos — click the canvas to launch trajectories from any initial condition.
The state of a dynamical system is completely determined by its position $x$ and velocity $\dot{x}$. The plane with these as coordinate axes is called phase space, and the time evolution of the system traces out a curve (trajectory) in this space.
Energy-conserving system. Phase portrait shows closed ellipses — nested concentric orbits that never cross. The ellipse shape encodes the amplitude and frequency of oscillation.
Energy-dissipating system. Trajectories spiral inward toward the origin (stable fixed-point attractor). With $\gamma < \omega$: underdamped spiral; with $\gamma > \omega$: overdamped node.
For $|x|<1$: negative damping (energy pumped in); for $|x|>1$: positive damping (energy removed). All trajectories converge to a stable limit cycle. Models vacuum tube oscillators and cardiac rhythms.
Nonlinear spring with periodic forcing. At suitable parameters produces a strange attractor and deterministic chaos. Launch two nearby trajectories to observe exponential divergence (butterfly effect).
Simplified atmospheric convection model. The 3D trajectory is projected onto the x-z plane, revealing the iconic "butterfly" strange attractor. The classic demonstration of deterministic chaos.
Student: How is a phase portrait different from a regular time-series plot?
Professor: In a time-series plot you see position versus time — a sine wave for the harmonic oscillator. In phase space, time is eliminated and you plot velocity against position. What you gain is a clear view of the system's qualitative structure: closed orbits mean periodicity, spirals mean decay, and complicated tangles mean chaos. A single glance at the phase portrait tells you more about the system's long-term behavior than pages of time-series data.
Student: So a closed ellipse in phase space always means the motion is periodic?
Professor: Exactly. The Poincare-Bendixson theorem tells us that in 2D phase space, if a trajectory stays bounded and doesn't approach a fixed point, it must approach a closed curve. That's why chaos is impossible in 2D autonomous systems — you need at least 3D (like Lorenz) or an explicit time-dependence (like Duffing with forcing) to get chaotic trajectories.
Student: Why does the Van der Pol limit cycle attract trajectories from both inside and outside?
Professor: The key is the $\mu(1-x^2)y$ term. When $|x|<1$, the factor $(1-x^2)$ is positive, so this term acts like \lt strong style="color:#7dd3fc"\gt negative damping — it injects energy and trajectories spiral outward. When $|x|>1$, the factor becomes negative, giving positive damping — trajectories spiral inward. The limit cycle is the unique closed orbit where these two effects exactly balance. It's a self-sustaining oscillation that doesn't need external forcing.
Student: Is that really how heartbeats work?
Professor: It's a useful model. The cardiac pacemaker cells fire autonomously at a steady rhythm, and if you perturb them slightly, they return to that rhythm — just like a limit cycle attractor. Try cranking up $\mu$ in the simulator. The waveform becomes sharper and more sawtooth-like, resembling a relaxation oscillation. Real ECG signals have a similar character during arrhythmia analysis.
Student: Is chaos the same as randomness? They both seem unpredictable.
Professor: This is a crucial distinction. Chaos is deterministic — the equations are completely fixed, and given identical initial conditions you get identical trajectories every time. The unpredictability comes solely from sensitivity to initial conditions: a tiny difference in starting point grows exponentially. Randomness, on the other hand, is intrinsic — even with perfect knowledge of initial conditions you can't predict the outcome. Weather forecasting fails beyond ~2 weeks not because the atmosphere is random, but because it's chaotic.
The core of the simulator is solving a set of first-order ordinary differential equations (ODEs) that define the system's dynamics. The general form for each system is:
$$ \dot{\mathbf{x}}= f(\mathbf{x}, t) $$where $\mathbf{x}$ is the state vector (e.g., $[x, y, z]$) and $\dot{\mathbf{x}}$ is its time derivative. The function $f$ contains the specific physics of the chosen model.
For example, the famous Lorenz system (modeling atmospheric convection) is defined by three coupled, nonlinear ODEs:
$$ \begin{align*}\dot{x}&= \sigma (y - x) \\ \dot{y}&= x (\rho - z) - y \\ \dot{z}&= xy - \beta z \end{align*} $$Here, $x$ is convection intensity, $y$ horizontal temperature variation, $z$ vertical temperature variation. The parameters $\sigma$ (Prandtl number), $\rho$ (Rayleigh number), and $\beta$ are constants. The nonlinear terms $xz$ and $xy$ are the source of the chaotic behavior.
Climate & Weather Modeling: The Lorenz system was originally derived from a simplified model of atmospheric convection. It demonstrates why long-term weather prediction is fundamentally limited—tiny errors in initial conditions grow exponentially, a hallmark of chaos studied with tools like this simulator.
Engineering Vibration Control: The Duffing oscillator (a nonlinear spring) models structures like aircraft wings or building frames that experience large deflections. Engineers use phase portraits to identify dangerous resonance conditions and chaotic vibrations that could lead to failure.
Electrical Circuit Design: The Van der Pol oscillator originally modeled circuits with vacuum tubes. Today, it's used to design and analyze oscillators in radios and pacemakers, where understanding the stability of limit cycles is critical for reliable operation.
Biological Rhythms: Nonlinear dynamics models the firing of neurons, the beating of a heart, and population cycles in ecology. Phase space analysis helps determine if a system will settle into a steady rhythm (limit cycle) or exhibit erratic, potentially pathological behavior.
When you start using this tool, there are a few common pitfalls. First is the case where "the trajectories look the same even when you change the initial conditions". For example, try turning off the forced oscillation and setting damping to zero for the "Duffing oscillator". Even if you change the initial conditions, all trajectories will lie on the same closed curve. This is because it's a conservative system. Since energy is conserved, the system keeps orbiting along a curve of constant energy (like a contour line), determined by the initial "push". The trajectories aren't overlapping; independent trajectories with different energies are drawing separate closed curves. If you add a little "damping" in the parameter settings, you can observe how these trajectories converge in a spiral pattern, so give it a try.
Next is over-reliance on numerical calculation accuracy. While this tool uses a high-precision method called RK4, especially in chaotic systems (like the Lorenz system), even tiny computational errors can be amplified explosively over time (that's the butterfly effect). Therefore, don't consider the detailed trajectories from long-term simulations as "absolute truth". In practice, it's fundamental to use it for observing qualitative behaviors like the "shape" of an attractor or its "properties" (whether it's stable or unstable).
Finally, there's "binge-changing" parameters. If you change a parameter drastically all at once, the system's behavior can change completely, often leaving you wondering "what just happened?". For instance, if you change the nonlinearity parameter $\mu$ for the van der Pol oscillator from 0.1 to 5.0 in one go, the shape changes dramatically. To understand the change, the trick is to change it little by little—0.1, 0.5, 1.0, 2.0…—and track how the limit cycle deforms continuously. This lets you get a feel for the "physical meaning of the parameter".
For the Lorenz system with σ=10, ρ=28, β=8/3: starting from (x₀=1, y₀=1, z₀=1), the RK4 integrator (dt=0.01s, 10,000 steps) converges to the famous butterfly attractor with fractal dimension ~2.06. At ρ=24.74, the system exhibits a period-doubling cascade; at ρ=28, fully chaotic motion dominates with Lyapunov exponent λ≈0.906 bits/s. The Duffing oscillator (α=1, β=1, γ=0.3, ω=1.2) shows period-1 limit cycles at low forcing amplitudes, transitioning to chaos around γ≈0.4 with basin-of-attraction boundary fractality.
