Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

Derive the equation of motion $m\ddot{x} + c\dot{x} + kx = 0$ from Newton's second law for a mass on a spring with viscous damping.
Re-cast any free mechanical vibration into the dimensionless standard form $\ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x = 0$ and read off the natural frequency $\omega_0$ and damping ratio $\zeta$ .
Classify the four regimes (undamped, underdamped, critically damped, overdamped) directly from $\zeta$ .
Write the closed-form solution in each regime and apply both initial conditions $x(0), \dot{x}(0)$ .
Connect the mathematical formulas to physical intuition: kinetic energy, potential energy, the damped period $T_d$ , and the logarithmic decrement $\delta$ .
Implement a vibration simulator in Python and a system-identification routine in PyTorch that recovers $(\omega_0, \zeta)$ from measured peaks.

The Big Picture: Why Vibrations?

"The whole of engineering vibration analysis is a one-page story: mass × acceleration plus damping × velocity plus stiffness × position equals zero — and the discriminant decides everything."

Look around. A guitar string after you pluck it. The bridge deck of a skyscraper in a gust. The needle of a galvanometer settling on its reading. A car wheel after you drive over a pothole. A drug molecule rattling around inside a binding pocket. Every one of these systems obeys the SAME second-order linear ODE with the SAME three coefficients. Get the math right once and you have understood half of mechanical engineering and a sizeable slice of physics.

The previous five sections built the toolkit: characteristic equation, three root cases, undetermined coefficients, variation of parameters. This section is the moment those abstract tools meet a physical reality you can feel under your hand.

The mechanical vibration ODE

A mass $m$ attached to a spring of stiffness $k$ , with a dashpot supplying viscous damping $c$ , satisfies

m\,\ddot{x} + c\,\dot{x} + k\,x = 0,\qquad x(0)=x_0,\;\; \dot{x}(0)=v_0.

Everything in this section follows from this single equation. No new physics is needed — only careful reading of the three coefficients.

🚗 Automotive

Suspension & shock absorbers
Engine mounts and exhaust hangers
Steering-wheel shimmy

🏗️ Civil & structural

Earthquake response of buildings
Bridge oscillations (Tacoma Narrows redux)
Tuned mass dampers in skyscrapers

🎵 Acoustics & instruments

String, drum, plate vibrations
Loudspeaker cone dynamics
Tuning-fork ring-down

🔬 Precision instruments

Atomic-force microscope cantilevers
MEMS gyroscopes and accelerometers
Galvanometer needles, seismometers

From Newton's Law to the ODE

Set up coordinates so that $x = 0$ is the spring's natural (unstretched) length. Pull the mass to the right by $x$ . Three forces act:

Force	Expression	Why this form
Spring (restoring)	F_spring = −k x	Hooke's law: spring force is proportional to displacement and points back toward equilibrium.
Viscous damping	F_damp = −c x'	Resistance from a dashpot or air drag at low speed: proportional to velocity, opposing motion.
External (zero for now)	F_ext = 0	Free vibration: no driving force. Section 22.7 will add F_ext = F₀ cos(ω t).

Plug into Newton's second law $m\ddot{x} = \sum F$ :

m\,\ddot{x} = -k\,x - c\,\dot{x}.

Move everything to the left-hand side to obtain the canonical form:

m\,\ddot{x} + c\,\dot{x} + k\,x = 0.

Sign of the damping term

The damping force is $-c\dot{x}$ because it always opposes the direction of motion. Moving to the right ( $\dot{x}>0$ ) ⇒ damping points left. Moving to the left ( $\dot{x}<0$ ) ⇒ damping points right. The minus sign is what guarantees that damping always drains energy from the system, never adds it.

Standard Form: Reading $\omega_0$ and $\zeta$

Divide the equation by $m$ :

\ddot{x} + \frac{c}{m}\,\dot{x} + \frac{k}{m}\,x = 0.

The two dimensional groups $k/m$ and $c/m$ have natural interpretations. Define

\omega_0 = \sqrt{\tfrac{k}{m}}\quad\text{(natural frequency, rad/s)}

\zeta = \dfrac{c}{2\sqrt{mk}}\quad\text{(damping ratio, dimensionless)}

With these substitutions the ODE becomes

\boxed{\;\ddot{x} + 2\zeta\omega_0\,\dot{x} + \omega_0^2\,x = 0.\;}

Two parameters now describe every linear free vibration in the universe: $\omega_0$ sets the time scale (how fast it wants to oscillate), and $\zeta$ sets the dissipation (how strongly the oscillation is suppressed).

Why this rescaling is worth doing

Two systems with very different masses, springs and dampers can have the SAME $\omega_0$ and $\zeta$ . They will then behave identically up to a change of time units. This is the same idea as the Reynolds number in fluid dynamics: collapse a high-dimensional design space onto a low-dimensional one. Engineers TUNE $\zeta$ directly.

The characteristic equation of the rescaled ODE is

r^2 + 2\zeta\omega_0\,r + \omega_0^2 = 0,

whose discriminant simplifies beautifully:

\Delta = (2\zeta\omega_0)^2 - 4\omega_0^2 = 4\omega_0^2(\zeta^2 - 1).

Sign of $\Delta$ ⇔ sign of $\zeta^2-1$ ⇔ position of $\zeta$ relative to 1. The famous $\zeta=1$ boundary is exactly the $\Delta=0$ knife edge from Section 22.1.

Free Undamped Vibrations ( $\zeta=0$ )

Set $c = 0$ . The equation collapses to

\ddot{x} + \omega_0^2 x = 0.

The characteristic roots are $r = \pm i\omega_0$ : pure imaginary. The general real solution is

x(t) = C_1 \cos(\omega_0 t) + C_2 \sin(\omega_0 t),

which any trigonometry textbook lets us repackage in amplitude–phase form:

x(t) = A\cos(\omega_0 t - \phi),\quad A = \sqrt{C_1^2 + C_2^2},\quad \tan\phi = C_2/C_1.

This is simple harmonic motion (SHM): a pure sinusoid of constant amplitude $A$ , angular frequency $\omega_0$ , and phase $\phi$ . The mass would oscillate forever — a physically idealised limit, never quite achievable in any real system because some damping is always present.

Period and frequency

The period is $T = 2\pi/\omega_0 = 2\pi\sqrt{m/k}$ . Heavy mass on weak spring ⇒ long period. Light mass on stiff spring ⇒ short period (high-pitched). The frequency in Hz is $f = 1/T = \omega_0/(2\pi)$ .

Damped Free Vibrations

Turn damping back on ( $\zeta>0$ ). The characteristic roots are

r_{1,2} = -\zeta\omega_0 \;\pm\; \omega_0\sqrt{\zeta^2 - 1}.

Notice the structure: there is always a real, negative common piece $-\zeta\omega_0$ that drives the envelope. The square root then decides whether the second piece is real or imaginary — exactly the discriminant classification.

ζ	Roots r₁, r₂	Solution form	Name
ζ = 0	±iω₀	C₁ cos(ω₀t) + C₂ sin(ω₀t)	Undamped (SHM)
0 < ζ < 1	−ζω₀ ± iω_d, with ω_d = ω₀√(1−ζ²)	e^(−ζω₀t)·(C₁ cos ω_d t + C₂ sin ω_d t)	Underdamped
ζ = 1	−ω₀ (repeated)	(C₁ + C₂t) e^(−ω₀t)	Critically damped
ζ > 1	−ζω₀ ± ω₀√(ζ²−1) (both real, negative)	C₁ e^(r₁ t) + C₂ e^(r₂ t)	Overdamped

Three observations worth memorising:

The exponential envelope $e^{-\zeta\omega_0 t}$ appears in EVERY stable regime. The time constant is $\tau = 1/(\zeta\omega_0)$ .
Damping slows oscillation: $\omega_d = \omega_0\sqrt{1-\zeta^2} < \omega_0$ . For light damping this shift is small (≈ $\zeta^2/2$ relative).
Critical damping ( $\zeta=1$ ) returns the system to equilibrium as fast as possible without overshoot. Door-closers and oven-thermostat positioners aim for this regime.

The Four Regimes — Side by Side

Each regime corresponds to one of the three discriminant cases, with $\zeta=0$ split out for clarity:

ζ = 0 — Undamped

Pure SHM. Pendulum in a vacuum, ideal LC tank, frictionless mass on a spring. Energy is conserved. The phase-portrait trajectory is a closed ellipse.

0 < ζ < 1 — Underdamped

Oscillates while decaying. Car suspensions, guitar strings, RLC circuits with light R. Phase trajectory: inward spiral. Most physical systems live here.

ζ = 1 — Critically damped

The boundary. Fastest non-oscillating return to equilibrium. Door closers, instrument indicators, the design target for precision positioning.

ζ > 1 — Overdamped

Sluggish, no oscillation. Heavy hydraulic dampers, screen-door closers in winter, signal probes designed for slew-rate-limited response.

Why the slow root dominates overdamped systems

For $\zeta>1$ the two real roots $r_1, r_2$ satisfy $r_1 r_2 = \omega_0^2$ and $r_1 + r_2 = -2\zeta\omega_0$ . The root closer to zero (smaller in magnitude) has the LONGER time constant and dominates the late-time behaviour. The faster root dies off quickly and matters mostly near $t=0$ .

Energy in Vibrations

Multiply the equation of motion by $\dot{x}$ :

m\dot{x}\ddot{x} + c\dot{x}^2 + k x \dot{x} = 0.

Recognise the first and third terms as time derivatives of the kinetic and potential energies:

\frac{d}{dt}\!\left(\tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2\right) = -c\dot{x}^2.

On the left is the total mechanical energy $E(t) = \tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2$ . On the right is the dissipation rate $-c\dot{x}^2 \le 0$ — never positive, zero only when the mass is instantaneously at rest. So:

Undamped (c = 0): $E$ is conserved. The phase-plane trajectory is the ellipse $\tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2 = E_0$ .
Damped: $E(t)$ decreases monotonically toward zero. Energy bleeds into heat at the rate $c\dot{x}^2$ .

For an underdamped vibration the envelope of $x(t)$ decays like $e^{-\zeta\omega_0 t}$ , so the envelope of $E(t)$ — quadratic in amplitude — decays like $e^{-2\zeta\omega_0 t}$ . Energy decays at twice the amplitude rate.

Logarithmic Decrement: Measuring Damping from Data

In the lab you cannot read off $\zeta$ directly. You can, however, hit the structure with a hammer and watch successive peaks of the response $x_0, x_1, x_2, \ldots$ at times $t_0, t_1=t_0+T_d, t_2=t_0+2T_d, \ldots$ .

Because the underdamped solution is an exponentially decaying sinusoid,

\frac{x_n}{x_{n+1}} = \frac{A e^{-\zeta\omega_0 t_n}}{A e^{-\zeta\omega_0 (t_n+T_d)}} = e^{\zeta\omega_0 T_d}.

Take the natural log:

\delta = \ln\!\frac{x_n}{x_{n+1}} = \zeta\,\omega_0\,T_d = \frac{2\pi\zeta}{\sqrt{1-\zeta^2}}.

This is the logarithmic decrement. It is independent of which two peaks you pick — every consecutive pair gives the same $\delta$ . From a measurement of two peaks you recover

\zeta = \frac{\delta}{\sqrt{4\pi^2 + \delta^2}}.

Practical estimator

For more accuracy use $N$ peaks and average: $\delta \approx \tfrac{1}{N}\ln(x_0/x_N)$ . This averages out individual measurement errors and is the basis of all standard ring-down damping tests (ASTM E756, ISO 2017-1).

Worked Example (Step-by-Step)

A 200 kg quarter-car (one wheel's share of an 800-kg sedan's sprung mass) is supported by a spring of stiffness 18 000 N/m and a shock absorber providing viscous damping $c=1200$ N·s/m. The car is lifted by 5 cm and released from rest. Describe the motion.

Click to expand the full hand calculation

Step 1 — Read off the coefficients.

From $m\ddot{x}+c\dot{x}+kx=0$ : $m=200$ , $c=1200$ , $k=18000$ . Initial conditions: $x(0)=0.05$ m, $\dot{x}(0)=0$ .

Step 2 — Compute $\omega_0$ .

$\omega_0 = \sqrt{k/m} = \sqrt{18000/200} = \sqrt{90} \approx 9.4868$ rad/s. In hertz this is $f_0 = \omega_0/(2\pi) \approx 1.510$ Hz.

Step 3 — Compute $\zeta$ .

$\zeta = c/(2\sqrt{mk}) = 1200/(2\sqrt{200 \cdot 18000}) = 1200/(2\sqrt{3{,}600{,}000}) = 1200/3794.7 \approx 0.3162$ . Since $0<\zeta<1$ the system is underdamped — exactly the regime a passenger car suspension is tuned for.

Step 4 — Damped frequency.

$\omega_d = \omega_0\sqrt{1-\zeta^2} = 9.487\sqrt{1-0.1} \approx 9.487 \cdot 0.9487 \approx 9.000$ rad/s. Damped period $T_d = 2\pi/\omega_d \approx 0.698$ s.

Step 5 — Real part of the roots.

$\alpha = -\zeta\omega_0 = -0.3162\cdot 9.487 \approx -3.000$ 1/s. The amplitude envelope is $e^{-3t}$ : every second the oscillation amplitude shrinks by a factor $e^{-3}\approx 0.0498$ — about 20×.

Step 6 — Plug into the underdamped form.

x(t) = e^{-3t}\bigl(C_1\cos(9t) + C_2\sin(9t)\bigr).

Step 7 — Apply $x(0)=0.05$ .

At $t=0$ : $x(0)=C_1=0.05$ .

Step 8 — Apply $\dot{x}(0)=0$ .

Differentiate: $\dot{x}(t) = e^{-3t}\bigl[(-3C_1 + 9C_2)\cos(9t) + (-3C_2 - 9C_1)\sin(9t)\bigr]$ . At $t=0$ : $\dot{x}(0) = -3C_1 + 9C_2 = 0$ ⇒ $C_2 = C_1/3 = 0.05/3 \approx 0.01667$ .

Step 9 — Particular solution.

\boxed{\,x(t) = e^{-3t}\bigl(0.05\cos(9t) + 0.01667\sin(9t)\bigr).\,}

Step 10 — Amplitude–phase form.

$A = \sqrt{0.05^2 + 0.01667^2} \approx 0.0527$ m. Phase $\phi = \arctan(0.01667/0.05) \approx 0.3217$ rad ≈ 18.4°. So $x(t) \approx 0.0527\,e^{-3t}\cos(9t - 0.32)$ m.

Step 11 — Logarithmic decrement check.

$\delta = 2\pi\zeta/\sqrt{1-\zeta^2} = 2\pi\cdot 0.3162/0.9487 \approx 2.094$ . Ratio of successive peaks: $e^{2.094} \approx 8.12$ . Each peak is about 8× smaller than the previous one — by the third bounce the oscillation is essentially gone.

Step 12 — Settling time.

Define settling when the envelope falls below 2% of its initial value: $e^{-3 t_s} = 0.02$ ⇒ $t_s = \ln(50)/3 \approx 1.30$ s. So the car stops bouncing in just over a second — exactly the comfort target of an automotive damper engineer.

Interactive Mass–Spring–Damper

Drag the sliders for $m, c, k$ , the initial displacement $x_0$ , and initial velocity $v_0$ . Watch the mass move and the trajectory plot itself in real time. The damping classification updates live as you cross the $\Delta = 0$ boundary.

Loading mass–spring–damper simulation…

Things to try

Set $c = 0$ . You will see a perfectly constant amplitude oscillation — undamped SHM.
Slowly raise $c$ . The amplitude envelope shrinks faster and faster, but the system still oscillates — that is the underdamped regime.
Tune $c$ to make $c^2 \approx 4mk$ . The oscillation just disappears. That is critical damping.
Push $c$ well past critical. The return to equilibrium gets slower — the overdamped regime is not faster, contrary to a common misconception.

Interactive: All Four Regimes Side-by-Side

The same natural frequency $\omega_0$ is kept across four curves with different damping ratios $\zeta \in \{0,\,0.2,\,1,\,2\}$ . Watch how the same initial perturbation evolves so differently. The right panel shows the phase portrait — position on the horizontal axis, velocity on the vertical axis.

Loading damping-regime comparator…

What to notice in the phase portrait

The undamped curve closes into an ellipse — energy conserved. The underdamped curve spirals inward, slowly. The critically damped and overdamped curves head straight for the origin without looping — no rotation, just attraction. The four geometries are the visual signatures of the four regimes.

Interactive: Logarithmic Decrement

Watch successive peaks of an underdamped trace shrink by a constant ratio $e^\delta$ . The table below the plot confirms that $\ln(x_n/x_{n+1})$ is the same for every consecutive pair — this is the property that makes $\delta$ a reliable experimental fingerprint of the damping.

Loading logarithmic-decrement visualizer…

Computation in Python

Math first, code second. The function below implements the three analytic branches (over-, critically-, and underdamped) using the closed-form solutions we just derived.

Closed-form vibration simulator

Direct closed-form evaluation of x(t) given (m, c, k, x₀, v₀)

🐍free_vibration.py

Explanation(15)

Code(68)

1Import NumPy

NumPy vectorises the math: np.exp, np.cos, and np.sin all act elementwise on the time grid t, so the full trajectory is built without an explicit Python loop. We also rely on np.sqrt for the discriminant.

4Function header

Four physical inputs (m, c, k, x0, v0) plus the time grid t. We need BOTH x(0) and x′(0) because a second-order ODE has a two-parameter family of solutions — one initial condition would leave the system under-determined.

22Natural angular frequency

ω₀ = √(k/m) is the angular frequency the system would oscillate at if there were no damping (c = 0). Its physical meaning: stiffer spring or lighter mass → faster oscillations. Units: rad / s. Multiply by 1/(2π) to get the frequency in Hz.

23Damping ratio ζ

ζ = c / (2√(mk)) is the dimensionless damping ratio. It is the single most informative number about the system. ζ = 0 → pure oscillation; 0 < ζ < 1 → underdamped (oscillates while decaying); ζ = 1 → critically damped (fastest return without overshoot); ζ > 1 → overdamped (sluggish). Engineers reach for ζ first, before m, c, or k individually.

26Discriminant

Δ = c² − 4mk is the discriminant of the characteristic polynomial m r² + c r + k = 0. Its SIGN — not its value — picks the qualitative regime. Notice that ζ² − 1 = Δ / (4mk), so Δ and ζ − 1 share a sign: ζ < 1 ⇔ Δ < 0, ζ = 1 ⇔ Δ = 0, ζ > 1 ⇔ Δ > 0. Two languages for the same boundary.

28Critical-damping branch

When Δ ≈ 0 the characteristic equation has a repeated root r = −ω₀. One exponential e^(rt) is not enough; the second linearly independent solution is t·e^(rt). The general solution is (C₁ + C₂t)e^(rt) — the SAME amplitude that grows linearly in t before the exponential beats it down.

32Apply both initial conditions (critical case)

Plug t = 0: C₁ = x(0). Differentiate (C₁ + C₂t)e^(rt) and set t = 0: x′(0) = C₂ + r·C₁, so C₂ = v0 − r·x0. Two unknowns, two equations, solved by inspection.

35Overdamped branch

When Δ > 0 the roots are two distinct real negatives r₁, r₂ (both negative since c, m, k are all positive and the polynomial is stable). Solution: x(t) = C₁ e^(r₁ t) + C₂ e^(r₂ t). No oscillation — the slower mode (root closer to zero) dominates after a few time constants.

38Quadratic formula

Standard quadratic roots r = (−c ± √Δ) / (2m). Because m > 0, dividing by 2m doesn't flip signs. Because Δ > 0 here, sqrtD is a real number.

41Solve 2×2 linear system for C₁, C₂

Two conditions: C₁ + C₂ = x0 and r₁C₁ + r₂C₂ = v0. Eliminate C₁ via C₁ = x0 − C₂ and solve for C₂. The denominator r₂ − r₁ is non-zero because the roots are distinct in this branch.

45Underdamped branch

When Δ < 0 the roots are α ± iω_d. The real part α = −c/(2m) = −ζω₀ controls the exponential ENVELOPE; the imaginary part ω_d controls the oscillation frequency. Euler's formula converts the complex exponential into a real combination of cos(ω_d t) and sin(ω_d t) wrapped by e^(αt).

47Damped natural frequency

ω_d = √(−Δ) / (2m) = ω₀·√(1 − ζ²). Always strictly less than ω₀ for ζ > 0: damping slows oscillations a little. For light damping (ζ small) the slowdown is tiny — to first order, ω_d ≈ ω₀(1 − ζ²/2).

50Initial conditions for underdamped case

At t = 0 the cosine is 1 and the sine is 0, so x(0) = C₁ directly gives C₁ = x0. Differentiating x(t) and setting t = 0 yields v0 = αC₁ + ω_d·C₂, so C₂ = (v0 − α·x0) / ω_d. ω_d is non-zero here because we are inside the Δ < 0 branch.

57Worked example — quarter-car suspension

We model one corner of a car: 200 kg of sprung mass, a spring of stiffness 18000 N/m, and a shock absorber of c = 1200 N·s/m. Release from rest 5 cm above equilibrium. Compute ω₀ = √(18000/200) = √90 ≈ 9.487 rad/s. ζ = 1200 / (2·√(200·18000)) = 1200 / 3794.7 ≈ 0.3162 — clearly underdamped, as a passenger-car suspension should be. The system will oscillate a few cycles and then settle.

62Sanity-check outputs

Print the regime classification so we can confirm the branch the code took. Also print max|x| (no overshoot above the initial pull-up, as expected since v0 = 0) and a value half-a-second in — by then about half a damped period has passed and x should already be on the negative side of equilibrium.

53 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def free_vibration(m, c, k, x0, v0, t):
5    """
6    Closed-form solution of  m*x'' + c*x' + k*x = 0
7    for free mechanical vibration.
8
9    Parameters
10    ----------
11    m  : float   mass (kg)
12    c  : float   viscous damping coefficient (N·s/m)
13    k  : float   spring stiffness (N/m)
14    x0 : float   initial displacement x(0)   (m)
15    v0 : float   initial velocity     x'(0)  (m/s)
16    t  : np.ndarray   time grid where x(t) is evaluated.
17
18    Returns
19    -------
20    x   : np.ndarray   displacement at each t
21    info: dict with derived parameters and the damping regime.
22    """
23    # Natural angular frequency and damping ratio
24    omega0 = np.sqrt(k / m)                # rad / s
25    zeta   = c / (2.0 * np.sqrt(m * k))    # dimensionless
26
27    # Discriminant of m r^2 + c r + k = 0
28    disc = c * c - 4 * m * k
29
30    if abs(disc) < 1e-12:
31        # Critical damping: repeated root r = -omega0
32        r  = -omega0
33        C1 = x0
34        C2 = v0 - r * x0
35        x  = (C1 + C2 * t) * np.exp(r * t)
36        regime = "critical"
37    elif disc > 0:
38        # Overdamped: two distinct real roots
39        sqrtD = np.sqrt(disc)
40        r1 = (-c + sqrtD) / (2 * m)
41        r2 = (-c - sqrtD) / (2 * m)
42        C2 = (v0 - r1 * x0) / (r2 - r1)
43        C1 = x0 - C2
44        x  = C1 * np.exp(r1 * t) + C2 * np.exp(r2 * t)
45        regime = "overdamped"
46    else:
47        # Underdamped: complex conjugate roots alpha +/- i*omega_d
48        alpha   = -c / (2 * m)
49        omega_d = np.sqrt(-disc) / (2 * m)
50        C1 = x0
51        C2 = (v0 - alpha * x0) / omega_d
52        x  = np.exp(alpha * t) * (C1 * np.cos(omega_d * t) +
53                                  C2 * np.sin(omega_d * t))
54        regime = "underdamped" if zeta > 1e-9 else "undamped"
55
56    info = dict(omega0=omega0, zeta=zeta, disc=disc, regime=regime)
57    return x, info
58
59
60# ---- worked example: a car suspension impulse ----
61# 200 kg quarter-car, k = 18 000 N/m, c = 1200 N s/m
62# Released from x(0) = 0.05 m with v(0) = 0.
63t = np.linspace(0, 2.0, 600)
64x, info = free_vibration(m=200, c=1200, k=18000, x0=0.05, v0=0.0, t=t)
65
66print(info)
67print("max |x|   :", np.max(np.abs(x)))
68print("x at 0.5s :", x[150])

Phase portrait via the matrix exponential

For plotting trajectories in the $(x, v)$ plane and for sanity-checking energy conservation, it is sometimes cleaner to reduce the second-order ODE to a first-order vector ODE $Y' = AY$ and apply the matrix exponential. Same physics, different bookkeeping.

Phase-space trajectory and energy via expm

🐍phase_portrait.py

Explanation(11)

Code(35)

1Import NumPy

We use NumPy arrays as 2-vectors holding the state Y = [x, v]. NumPy makes the trajectory a (T, 2) array we can slice column-wise.

3Function header

Same physical inputs (m, c, k, x0, v0) and the time grid. We will return BOTH x(t) and v(t) — the full state — so the caller can draw a phase portrait or compute energies.

14Import scipy.linalg.expm

expm computes the dense matrix exponential. For a linear constant-coefficient system Y′ = AY the EXACT solution is Y(t) = e^(At) Y(0). No integration needed — this is closed form. The same routine powers most linear ODE solvers under the hood when they detect linearity.

15Build the companion matrix A

From m x″ + c x′ + k x = 0 introduce v = x′. Then x′ = v and v′ = −(k/m)x − (c/m)v. Stacking these gives the 2×2 matrix A. The top row [0, 1] is generic (because x′ = v); the bottom row encodes the original ODE. Key fact: eigenvalues of A are EXACTLY the roots r₁, r₂ of m r² + c r + k = 0. So all the physics we derived analytically lives in the spectrum of A.

19Initial state vector

Pack the two initial conditions into one 2-vector Y(0) = [x0, v0]. From here, the problem is one linear-algebra operation per time step.

22Vectorised matrix exponential

For each time t_i compute expm(A·t_i) — a 2×2 matrix — and multiply it by Y(0) to obtain the state Y(t_i). The list-comprehension builds a (T, 2) array. (For very long grids you would batch via scipy.linalg.expm_multiply or torch.linalg.matrix_exp, but for 600 points this is fast enough.)

23Return x and v separately

Slice column 0 for displacement and column 1 for velocity. The pair (x(t), v(t)) traces a curve in the phase plane — an inward spiral for underdamped systems, a direct attack on the origin for critically/overdamped systems, a closed ellipse for the undamped case.

31Compute kinetic energy along the trajectory

Classical KE = (1/2) m v². NumPy squares the v array elementwise, multiplies by m/2, and gives KE on the full grid in one expression. No loop.

32Compute potential energy along the trajectory

For a linear spring the elastic potential energy is (1/2) k x². The spring stores energy quadratically in displacement — pull twice as far, store four times the energy.

33Total mechanical energy

E(t) = KE(t) + PE(t). For an UNDAMPED system this would be exactly constant — a beautiful conservation law and the cleanest way to verify a numerical solver. For damped systems E(t) is monotonically decreasing because the dashpot dissipates energy at the rate c·v² (always non-negative).

35Read off energy loss

E(0.5) / E(0) tells us what fraction of the initial energy remains after 0.5 s. For the quarter-car with ζ ≈ 0.316 we expect roughly e^(−2ζω₀·0.5) = e^(−0.316·9.487·1) ≈ e^(−3) ≈ 0.05 — about 95% of the energy is gone, exactly the design goal of a passenger-car shock absorber: kill the energy in a fraction of a second so the cabin doesn't keep wobbling.

24 lines without explanation

1import numpy as np
2
3def vibration_state(m, c, k, x0, v0, t):
4    """
5    Return (x(t), v(t)) for the same ODE  m*x'' + c*x' + k*x = 0,
6    so we can plot the state in the phase plane (x, v).
7
8    Trick: rewrite as a first-order vector ODE
9        Y' = A Y  with  Y = [x, v],
10        A = [[ 0,    1   ],
11             [-k/m, -c/m ]]
12    Closed-form solution:  Y(t) = expm(A t) @ Y(0).
13    """
14    from scipy.linalg import expm
15    A  = np.array([[0.0,        1.0],
16                   [-k / m,    -c / m]])
17    Y0 = np.array([x0, v0], dtype=float)
18
19    # Stack expm(A*t_i) @ Y0 for each t_i  (vectorised over the grid)
20    states = np.array([expm(A * ti) @ Y0 for ti in t])
21    return states[:, 0], states[:, 1]   # x(t), v(t)
22
23
24# Same quarter-car parameters as before
25t   = np.linspace(0, 2.0, 600)
26x, v = vibration_state(m=200, c=1200, k=18000, x0=0.05, v0=0.0, t=t)
27
28# Energies along the trajectory
29KE = 0.5 * 200 * v**2          # 1/2 m v^2
30PE = 0.5 * 18000 * x**2        # 1/2 k x^2
31E  = KE + PE
32
33print("E(0)   :", E[0])
34print("E(0.5) :", E[150])
35print("Energy lost fraction in 0.5s :", 1 - E[150]/E[0])

PyTorch View: System Identification from Measurements

So far the recipe was parameters → trajectory. In practice you often have the opposite problem: given a measured trajectory, find the parameters. For an underdamped free vibration, the simplest such procedure is the logarithmic decrement fit: a linear regression of $\log|x_n|$ against $t_n$ . PyTorch tensors make this both fast and compatible with autograd, so the same routine can sit inside a larger learning pipeline (a neural physics model, a Kalman filter, a differentiable simulator).

Recover (ω₀, ζ) from successive peaks of a noisy measurement

🐍fit_damping.py

Explanation(12)

Code(58)

1Import PyTorch

We use torch tensors so the same code can later run on GPU and inside a learning loop. For a small least-squares fit this is overkill; the value of doing it in torch is that the function becomes drop-in for a larger differentiable pipeline (system-identification inside a neural model).

3Function header — system identification

Real engineers do not get to read off (m, c, k) — they observe a measurement of x(t) and have to INFER the parameters. Two well-known peaks (t_n, x_n) are enough to estimate ω₀ and ζ. This is the simplest non-trivial inverse problem of vibration analysis.

16Convert inputs to float64 tensors

torch.as_tensor avoids copying when the input is already a tensor. float64 gives the precision we want for log-arithmetic; with float32, log of a 10⁻⁴ value already loses several digits.

17Take absolute value of the peaks

Successive peaks alternate sign for an underdamped vibration (positive maxima and negative minima both count as 'peaks'). The envelope A·e^(−ζω₀t) is positive, so we fit |x_n|, not x_n.

20Linearise: take the log

If x_n ≈ A · e^(−ζω₀ t_n) then log|x_n| ≈ log A − (ζω₀) t_n — a STRAIGHT line in (t_n, log|x_n|) space. Slope = −(ζω₀), intercept = log A. The model is exact in the noiseless case and well-approximated under small additive noise.

22Closed-form least-squares slope

Classical formulas: t_mean = mean(t), y_mean = mean(y), b = S_xy / S_xx, a = y_mean − b·t_mean, where S_xy = Σ(t−t_mean)(y−y_mean) and S_xx = Σ(t−t_mean)². No need for torch.linalg.lstsq or autograd here — the equations have a clean closed form.

29Decode slope as damping rate

b is the slope of log|x| vs t. We modelled this as −ζω₀, so ζ·ω₀ = −b. This combined product is sometimes called the DAMPING RATE σ. Larger σ ⇒ faster amplitude decay, regardless of how that splits between ζ and ω₀.

30Recover damped period T_d from peak spacing

For an underdamped vibration the peaks are equally spaced by T_d = 2π/ω_d. We use the average inter-peak gap (more robust to noise than a single difference). From T_d we get ω_d immediately.

32Disentangle ω₀ and ζ

We know ω_d = ω₀·√(1−ζ²) and σ = ζω₀. Squaring and adding: ω_d² + σ² = ω₀² (1−ζ²) + ω₀² ζ² = ω₀². So ω₀ = √(ω_d² + σ²) — a Pythagorean relationship between the oscillation frequency, the decay rate, and the natural frequency. Then ζ = σ/ω₀ closes the loop.

41Simulate ground-truth data

We pretend we observed a real system with ω₀ ≈ 9.487 and ζ ≈ 0.316 (the quarter-car from earlier). We generate the noiseless peak times and amplitudes, sprinkle in some Gaussian measurement noise, then check whether the algorithm recovers the truth.

50Add measurement noise

0.0005 ≈ 1% of the starting amplitude. torch.randn_like draws standard normal noise of the same shape and dtype. This makes the test realistic: in the lab you never see clean exponentials. A good estimator must still pin down ω₀ and ζ within a few percent under noise.

54Run the estimator and print

We expect ω₀ to come back as ≈ 9.487 and ζ as ≈ 0.316. Any drift larger than the noise floor would indicate a bug. This is the same diagnostic an experimentalist runs after a 'ring-down' test: hit the structure, record the trace, fit the envelope, read off ζ and ω₀.

46 lines without explanation

1import torch
2
3def fit_damping_from_peaks(t_peaks, x_peaks):
4    """
5    Given the times t_n and amplitudes x_n of successive peaks
6    of an underdamped free vibration, recover (omega0, zeta) by
7    least squares.
8
9    Idea: the envelope of x(t) is A * exp(-zeta * omega0 * t).
10    Taking logs:
11        log |x_n| = log A  -  (zeta * omega0) * t_n.
12    A linear regression of log|x_n| against t_n recovers the
13    slope -(zeta * omega0). The damped period T_d = t_{n+1} - t_n
14    gives omega_d = 2 pi / T_d, then omega0 = sqrt(omega_d^2 +
15    (zeta*omega0)^2).
16    """
17    t  = torch.as_tensor(t_peaks, dtype=torch.float64)
18    x  = torch.as_tensor(x_peaks, dtype=torch.float64).abs()
19
20    # Linear fit: log|x| = a + b*t,  b = -zeta*omega0
21    y = torch.log(x)
22    n = t.numel()
23    t_mean = t.mean()
24    y_mean = y.mean()
25    Sxx = ((t - t_mean) ** 2).sum()
26    Sxy = ((t - t_mean) * (y - y_mean)).sum()
27    b = Sxy / Sxx
28    a = y_mean - b * t_mean
29
30    zeta_omega0 = -b                             # damping rate
31    T_d         = (t[-1] - t[0]) / (n - 1)       # mean inter-peak gap
32    omega_d     = 2 * torch.pi / T_d
33    omega0      = torch.sqrt(omega_d ** 2 + zeta_omega0 ** 2)
34    zeta        = zeta_omega0 / omega0
35    return dict(omega0=omega0, zeta=zeta, T_d=T_d,
36                A=torch.exp(a))
37
38
39# Synthesise data from a system we know, then recover its parameters
40torch.manual_seed(0)
41omega0_true = 9.487
42zeta_true   = 0.316
43omega_d     = omega0_true * (1 - zeta_true ** 2) ** 0.5
44T_d         = 2 * torch.pi / omega_d
45
46# 6 successive peaks
47n_peaks = 6
48t_peaks = torch.tensor([k * T_d for k in range(n_peaks)],
49                       dtype=torch.float64)
50x_peaks = 0.05 * torch.exp(-zeta_true * omega0_true * t_peaks)
51
52# Add a little noise to mimic measurement error (1% of amplitude)
53x_peaks = x_peaks + 0.0005 * torch.randn_like(x_peaks)
54
55fit = fit_damping_from_peaks(t_peaks, x_peaks)
56print("recovered omega0 :", fit["omega0"].item())
57print("recovered zeta   :", fit["zeta"].item())
58print("recovered T_d    :", fit["T_d"].item())

From log-decrement fitting to differentiable simulators

The same idea generalises: write the full underdamped expression $x(t;\omega_0,\zeta,A,\phi)$ , build a loss $\sum_i (x(t_i)-x_i^{\text{meas}})^2$ , and minimise via $\texttt{torch.optim.Adam}$ . Modern ML libraries treat physical parameters as just another set of weights — system identification becomes a special case of training. The starting point is always the analytic understanding you developed in this section.

Common Pitfalls

Confusing ω₀ with ω_d

They differ by a factor of $\sqrt{1-\zeta^2}$ . For lightly damped systems they are nearly equal, but for $\zeta = 0.5$ the damped frequency is already about 13% lower than $\omega_0$ . Always be clear which one your formula or measurement refers to.

Forgetting the second initial condition

A second-order ODE needs both $x(0)$ and $\dot{x}(0)$ . Many problems specify only $x(0)$ and tacitly assume $\dot{x}(0)=0$ ("released from rest"). Read the wording carefully.

Critical damping is not the fastest decay

The amplitude envelope $e^{-\zeta\omega_0 t}$ appears to grow faster (faster envelope decay) as $\zeta$ grows. But for $\zeta>1$ the SLOWER of the two real roots dominates and overall settling becomes slower. Critical damping gives the fastest non-oscillating return; lightly underdamped systems can reach equilibrium even faster, at the cost of a small overshoot.

Don't reach for natural log if x_n is negative

In the logarithmic-decrement formula $\delta = \ln(x_n/x_{n+1})$ , both $x_n$ and $x_{n+1}$ must be SAME-SIGN peaks of the envelope. If you mix a positive peak with a negative trough the ratio is negative and the log is undefined. Use $|x_n|$ and only pair maxima with maxima (or minima with minima).

Summary

Every free linear mechanical vibration is governed by the same two-parameter equation

\ddot{x} + 2\zeta\omega_0\,\dot{x} + \omega_0^2\,x = 0,

and the damping ratio $\zeta$ alone classifies the response:

ζ	Regime	x(t)	Key quantity
0	Undamped (SHM)	A cos(ω₀ t − φ)	Period T = 2π/ω₀
0 < ζ < 1	Underdamped	A e^(−ζω₀t) cos(ω_d t − φ)	Log decrement δ = 2πζ/√(1−ζ²)
ζ = 1	Critically damped	(C₁ + C₂t) e^(−ω₀t)	Fastest non-oscillating return
ζ > 1	Overdamped	C₁ e^(r₁t) + C₂ e^(r₂t)	Dominant time constant 1/\|r_slow\|

Three working tools to remember:

Standard form. Divide by $m$ and read off $\omega_0$ and $\zeta$ . These two numbers carry all the qualitative information.
Energy. $E(t) = \tfrac{1}{2}m\dot{x}^2 + \tfrac{1}{2}kx^2$ decays at rate $c\dot{x}^2$ ; the envelope of $E$ falls twice as fast as the envelope of $x$ .
Logarithmic decrement. Measure two peaks, take a log of their ratio, recover $\zeta$ . This is the bread-and-butter experimental test of damping.

The slogan

\omega_0

sets the rhythm,

\zeta

sets the fate."

Coming next: Section 22.07 reintroduces a driving force $F(t) = F_0\cos(\omega t)$ . The result is the magnificent phenomenon of resonance: when the drive frequency approaches $\omega_0$ the steady-state amplitude blows up — gently if $\zeta$ is large, catastrophically if $\zeta$ is small. The Tacoma Narrows bridge, a child on a swing, and an NMR machine are all variations on this single theme.