Mass-Spring-Damper System

Governing Equation

The simulation models a mass \(m\) attached to a spring with natural (rest) length \(L_0\) [m] and stiffness \(k\), subject to viscous damping \(c\). Let \(x\) denote the total extension of the spring from the fixed support. The forces acting on the mass are the spring restoring force \(-k(x-L_0)\), viscous damping \(-c\dot{x}\), and gravity \(mg\) (downward). Newton's second law gives the full equation of motion:

The static equilibrium extension satisfies \(k(x_{\rm eq}-L_0)=mg\), so \(x_{\rm eq}=L_0+mg/k\). Defining the displacement from static equilibrium \(y=x-x_{\rm eq}\) and noting \(\ddot{y}=\ddot{x}\), the equation reduces to the homogeneous form

The rest length \(L_0\) does not appear in the reduced equation — it sets the visual scale of the simulation and, together with \(m\), \(g\), and \(k\), determines the equilibrium position shown by the dashed line. Larger \(L_0\) provides more canvas space for large-amplitude oscillations.

Characteristic Equation and Eigenvalues

Substituting \(y=e^{\lambda t}\) yields the characteristic equation \(m\lambda^2+c\lambda+k=0\), with roots

The discriminant \(\Delta=c^2-4mk\) determines the qualitative behaviour:

Three Regimes

• Underdamped \((\Delta<0,\;\zeta<1)\): The roots are complex conjugates \(\lambda=-\alpha\pm i\omega_d\) where \(\alpha=c/(2m)\) and \(\omega_d=\sqrt{k/m-\alpha^2}\). The solution oscillates with exponentially decaying amplitude: \(y(t)=e^{-\alpha t}(A\cos\omega_d t+B\sin\omega_d t)\).
• Critically damped \((\Delta=0,\;\zeta=1)\): Repeated real root \(\lambda=-c/(2m)\). The mass returns to equilibrium as fast as possible without oscillating: \(y(t)=(A+Bt)e^{-ct/(2m)}\).
• Overdamped \((\Delta>0,\;\zeta>1)\): Two distinct negative real roots. The return to equilibrium is slower than critical and without oscillation.

The damping ratio is \(\zeta=c/(2\sqrt{mk})\) and the natural frequency is \(\omega_0=\sqrt{k/m}\).

Numerical Integration

• Velocity Verlet (Störmer): Symplectic integrator used for non-stiff systems (\(\Delta\le5\)). Preserves the symplectic structure and provides excellent long-time energy conservation in undamped or lightly damped cases. Fixed step \(\Delta t=0.005\,\text{s}\), advanced by real elapsed time so playback is frame-rate independent.
• Backward Euler: A-stable implicit method used for the first step when the system is stiff. Avoids the exponential blow-up that explicit methods suffer for large damping.
• BDF-2 (Backward Differentiation Formula, order 2): Second-order A-stable implicit method used for subsequent steps in the stiff regime. Provides better accuracy than backward Euler while retaining unconditional stability.
• Gauss–Legendre (2-stage): A 4th-order, symplectic, A-stable implicit Runge–Kutta method, selectable from the Numerical method menu. Because the free system is linear, each step reduces to one fixed \(2\times2\) matrix product (the \((2,2)\)-Padé propagator \(M=(I+\tfrac{X}{2}+\tfrac{X^2}{12})(I-\tfrac{X}{2}+\tfrac{X^2}{12})^{-1}\), \(X=\Delta t\,A\)), so it is as fast as Verlet yet essentially exact.

How to Use

• Adjust \(m,c,k,L_0\) and initial conditions using the value boxes. The live equation updates automatically.
• Increase \(L_0\) to lengthen the spring and allow larger initial displacements \(y_0\).
• Click Reset to restart from the current parameters. Use Pause/Resume to freeze the animation.
• Watch the displacement graph to distinguish the three regimes — decaying oscillations (underdamped), smooth exponential decay (critically/overdamped).
• Try \(c=0\) for pure simple harmonic motion, or increase \(c\) until \(c^2=4mk\) for critical damping.