Harmonic Oscillator Simulation

Free oscillation and the three damping regimes

Newton's second law for a mass m on a spring of stiffness k, with a damper that resists motion in proportion to velocity (coefficient c) and an optional driving force:

Two numbers describe the free motion (F₀ = 0). The natural frequency ω₀ = √(k/m) is how fast the spring would swing the mass without friction; the period is 2π/ω₀, independent of the amplitude. The damping ratio ζ = c / (2√(km)) compares friction with the spring:

• ζ = 0, undamped: a pure sine wave forever; energy ½mv² + ½kx² is conserved and the phase portrait is a closed ellipse.
• 0 < ζ < 1, underdamped: the mass still oscillates, at the slightly lower frequency ω_d = ω₀√(1 − ζ²), inside an envelope that shrinks like e^−ζω₀t. The phase portrait spirals inward.
• ζ = 1, critically damped: no oscillation; the fastest return to rest without overshooting. Door closers and analog meter needles are tuned near it.
• ζ > 1, overdamped: no oscillation and a slower, creeping return — the damper resists so strongly that the spring cannot pull the mass back quickly.

With the default m = 1 kg and k = 4 N/m, critical damping is at c = 2√(km) = 4 N·s/m; try it, then 8.

Driving force and resonance

Switch on F₀ and the motion becomes the sum of two parts: a transient that behaves like the free oscillation and dies out like e^−ζω₀t, and a steady oscillation at the driving frequency ω with amplitude

The resonance curve plots A·k/F₀, the amplitude relative to the static stretch F₀/k. Far below ω₀ the mass follows the force; far above it, inertia wins and the motion shrinks. Near ω₀ the curve peaks, higher and narrower the smaller ζ is. Driven exactly at ω₀ the amplitude is F₀/(cω₀): only damping limits it. With c = 0 and ω = ω₀ there is no steady state at all — the amplitude grows linearly for as long as the drive lasts.

The lab measures the steady-state amplitude from the last three driving periods and compares it with the formula. Set F₀ = 1 N with the defaults: after 30 s the simulated amplitude is close to the predicted 1.25 m, with a small leftover transient; lengthen the run and the two agree to more digits.

How the simulation checks itself

The equation is integrated with the classical fourth-order Runge–Kutta method (RK4), with a fixed step h = 0.02/ω_max, where ω_max is the largest of the natural frequency, the drive frequency and the damping rate c/m, and never more than 0.01 s. Alongside position and velocity the integrator accumulates the energy dissipated by the damper, ∫cv² dt, and the work done by the drive, ∫F·v dt. Energy balance requires

and the lab reports the largest violation relative to the largest energy involved. For the undamped, undriven oscillator this is simply the energy drift. RK4 does not conserve energy exactly: for a pure oscillator each step multiplies the energy by about 1 − (ω₀h)⁶/72, a tiny but systematic loss. With ω₀h ≤ 0.02 that is below 10⁻¹² per step. The simpler Euler method instead multiplies the energy by 1 + (ω₀h)² every step, so an undamped oscillator spirals outward — see Euler vs Runge–Kutta.

Example result

With m = 1 kg, k = 4 N/m and c = 0.4 N·s/m the natural frequency is ω₀ = 2.000 rad/s and the damping ratio ζ = 0.100, so the motion is underdamped. Released from x₀ = 1 m it oscillates with a simulated period of 3.157 s, matching the damped period 2π/ω_d = 3.157 s rather than the undamped 3.142 s. The energy balance holds to 3.2e-9.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator. The excerpts show the core of each.

// m·x'' + c·x' + k·x = F₀·cos(ωt), as a first-order system in (x, v). Two extra
// components integrate the energy lost to damping, ∫c·v² dt, and the work done by the
// drive, ∫F·v dt, so the energy balance can be checked at every step.
export function derivatives(t, s, p) {
  const [x, v] = s;
  const f = p.f0 * Math.cos(p.omega * t);
  return [v, (f - p.c * v - p.k * x) / p.m, p.c * v * v, f * v];
}

export function rk4Step(t, s, h, p) {
  const k1 = derivatives(t, s, p);
  const k2 = derivatives(t + h / 2, s.map((y, i) => y + (h / 2) * k1[i]), p);
  const k3 = derivatives(t + h / 2, s.map((y, i) => y + (h / 2) * k2[i]), p);
  const k4 = derivatives(t + h, s.map((y, i) => y + h * k3[i]), p);
  return s.map((y, i) => y + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]));
}

// Mechanical energy: kinetic plus spring potential.
export const energy = (x, v, p) => 0.5 * p.m * v * v + 0.5 * p.k * x * x;

def derivatives(t, s, p):
    """m x'' + c x' + k x = F0 cos(omega t), plus the running integrals of c v^2 and F v."""
    x, v = s[0], s[1]
    f = p["f0"] * math.cos(p["omega"] * t)
    return [v, (f - p["c"] * v - p["k"] * x) / p["m"], p["c"] * v * v, f * v]


def rk4_step(t, s, h, p):
    k1 = derivatives(t, s, p)
    k2 = derivatives(t + h / 2, [y + h / 2 * k1[i] for i, y in enumerate(s)], p)
    k3 = derivatives(t + h / 2, [y + h / 2 * k2[i] for i, y in enumerate(s)], p)
    k4 = derivatives(t + h, [y + h * k3[i] for i, y in enumerate(s)], p)
    return [y + h / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i, y in enumerate(s)]


def energy(x, v, p):
    """Kinetic plus spring potential energy."""
    return 0.5 * p["m"] * v * v + 0.5 * p["k"] * x * x

Run the full model from a terminal:

node simulation.js --m 1 --k 4
python simulation.py --m 1 --k 4

Download the model

Assumptions & limitations

• Linear spring. The restoring force is exactly proportional to displacement (Hooke's law) at every amplitude; real springs stiffen or soften, and pendulums are only linear for small angles.
• Viscous damping. Friction is proportional to velocity. Dry (Coulomb) friction and air drag at speed behave differently.
• Point mass, one dimension. The spring is massless and the motion is along a line.
• Ideal drive. The driving force is a perfect sinusoid unaffected by the motion it causes.

Related simulations

• OS-036Double pendulumTwo nearly identical double pendulums that part ways: sensitive dependence on initial conditions.
• OS-034Projectile motionThe vacuum parabola against a trajectory with quadratic air drag, and why 45° stops being the best angle.
• OS-030Lotka–VolterraTwo coupled equations whose predator and prey populations circle forever on closed orbits.

Background: Euler vs Runge–Kutta: Time Stepping in Simulations

References

1. Taylor, J. R. (2005). Classical Mechanics, chapter 5 (oscillations). University Science Books.
2. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos, 2nd ed. Westview Press.
3. Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007). Numerical Recipes, 3rd ed., chapter 17. Cambridge University Press.