From ODEs to State-Space

Chapter 1: Three Forces, One Equation

In Chapter 0 you felt the system. Now we make it precise — not by quoting a formula, but by writing down the forces one at a time and adding them up. By the end of this chapter you will have built the equation of motion from scratch, and every symbol in it will mean something physical to you.

Let x be the position of the block, measured from its rest position — the spot where the spring is neither stretched nor squished. So x = 0 means "at rest position," x > 0 means "pulled to the right," x < 0 means "pushed to the left." This choice (measuring from rest) is what keeps the math clean; it's worth a sentence of thought.

We need two more quantities, and they are both just rates of change of x. The velocity is how fast position changes, written x′ (read "x-prime") or dx/dt. The acceleration is how fast velocity changes, written x″ (read "x-double-prime") or d²x/dt². If those primes feel scary, hold this anchor: x′ is the speedometer, x″ is the gas-and-brake pedal. That's all a derivative is here — a rate.

Force 1: The spring (the restoring force)

Stretch a spring twice as far, it pulls back twice as hard. This proportionality is Hooke's Law, and the constant of proportionality is the stiffness k (units: newtons per meter — how many newtons of pull you get per meter of stretch). A stiff spring (large k) yanks hard; a floppy spring (small k) barely tugs.

The crucial detail is the direction. The spring always pulls back toward rest. Pull the block right (x > 0) and the spring pulls left (negative force). Push it left (x < 0) and the spring pulls right (positive force). The force always opposes the displacement — that opposing sign is what makes it a restoring force. So:

Read that minus sign as a sentence: "the further you go from home, the harder home pulls you back." That single restoring term is the entire reason the system can oscillate at all. Without it (k = 0) there's nothing to pull the block back, so it never returns — it just sits wherever the damper lets it stop.

Force 2: The damper (the friction-like force)

A damper is a piston pushed through thick oil — exactly like the closer at the top of a slow-shutting door. The faster you try to move it, the harder it resists. It doesn't care where the block is; it cares how fast the block is moving. So the damping force is proportional to velocity x′, with constant b (the damping coefficient, units newton-seconds per meter).

And again the direction matters: the damper always opposes the motion. Move right (x′ > 0) and it pushes left; move left (x′ < 0) and it pushes right. The force opposes velocity, so it carries a minus sign:

Here is the physical heartbeat of the whole lesson: the damper is the only term that removes energy. The spring stores energy and gives it all back. The mass carries energy but doesn't create or destroy it. Only the damper, by always pushing against motion, drains the system's energy into heat. That is why b decides whether the system rings forever or settles.

Force 3: Inertia (the mass term)

The last "force" isn't really an applied force at all — it's the left side of Newton's law. The mass m resists changes in velocity. Whatever the net force is, it produces acceleration through m·x″. A heavy block (large m) is sluggish; it accelerates slowly. A light block snaps around. We'll move this term to the same side as the others in a moment.

Assembling the equation of motion

Newton: mass × acceleration = sum of forces. The only two forces on our block (it's a frictionless table; the only contacts are the spring and the damper) are the spring force and the damping force:

Now we do one tiny piece of algebra: move every term to the left side so the equation equals zero. Add k·x and b·x′ to both sides:

That's it. That is the equation of motion, and you just built it from three physical sentences. Read it left to right as a story: m·x″ (inertia fighting the change) + b·x′ (damping fighting the speed) + k·x (spring fighting the displacement) = 0 (they're in balance at every instant). The "= 0" means no one is pushing the block from outside — it's just the three internal effects negotiating among themselves. This is called the free response or homogeneous equation ("homogeneous" is jargon for "right side is zero, no external driving force").

A worked sanity check with real numbers

Let's make sure every term is concrete. Take m = 2 kg, b = 3 N·s/m, k = 8 N/m. Suppose at some instant the block is at x = 0.5 m (pulled right half a meter) and moving at x′ = −1.2 m/s (moving left at 1.2 m/s). What is its acceleration right now?

Compute the spring force: F_spring = −k·x = −8 × 0.5 = −4 N. The spring pulls left (negative), because the block is to the right of rest. Good — that matches intuition.

Compute the damping force: F_damp = −b·x′ = −3 × (−1.2) = +3.6 N. The damper pushes right (positive), because the block is moving left and the damper opposes that motion. Also matches intuition.

Net force = −4 + 3.6 = −0.4 N. Now Newton: acceleration x″ = (net force) / m = −0.4 / 2 = −0.2 m/s². So at this instant the block is accelerating gently to the left. We computed it from raw forces — no solution formula needed. The equation of motion is just bookkeeping for forces, evaluated instant by instant. (Note this matches plugging into m·x″ + b·x′ + k·x = 0: 2×(−0.2) + 3×(−1.2) + 8×0.5 = −0.4 − 3.6 + 4 = 0. The equation balances. ✓)

The widget below lets you do this for any state. Drag the block to set x, drag the velocity arrow to set x′, and watch the three force arrows update live, along with the resulting acceleration. This is the equation of motion made visible: three arrows that must sum to m·x″.

Chapter 2: Guessing the Answer, and the Four Regimes

We have an equation: m·x″ + b·x′ + k·x = 0. But an equation tells us the rule the motion obeys at every instant; it doesn't yet tell us the trajectory — the actual x(t) curve of position over time. Solving for that curve is the heart of this chapter, and the method is gloriously sneaky: we're going to guess.

The educated guess: what function survives differentiation?

Look at the equation again. We need a function x(t) such that adding up a copy of itself (the kx term), a copy of its first derivative (the bx′ term), and a copy of its second derivative (the mx″ term) gives exactly zero. For that cancellation to be possible at every instant, the function and its derivatives must all have the same shape — otherwise they could never cancel for all t.

What function is its own derivative (up to a constant)? The exponential. If x(t) = e^λt, then x′ = λ·e^λt and x″ = λ²·e^λt. Every derivative is just the original times another factor of λ. That self-similarity is precisely the property we need. So we make the guess (called an ansatz, German for "an attempt"):

Here λ (lambda) is an unknown number we're solving for. It controls everything: a negative λ makes e^λt decay toward zero (good — the block settles); a positive λ makes it blow up (a system that explodes); and as we'll see, a complex λ secretly hides a sine wave (oscillation). One symbol, all the behaviors.

Plug in and watch the magic cancellation

Substitute the guess into m·x″ + b·x′ + k·x = 0:

Every single term has a factor of e^λt. Factor it out:

Now the key observation: e^λt is never zero — an exponential never touches the axis. So the only way the product can be zero is if the other factor is zero:

This is the characteristic equation, and it is a thing of beauty. We started with a differential equation (the hard, calculus-flavored kind) and, by the single act of guessing an exponential, we've reduced it to an ordinary quadratic equation — the kind you solved in high school. The calculus has evaporated. All the dynamics are now hidden in the roots of a quadratic.

Solving the quadratic: the roots

A quadratic aλ² + bλ + c = 0 has the famous solution λ = (−b ± √(b² − 4ac)) / (2a). Here a = m, the linear coefficient is b, and c = k. Substituting:

Two roots, λ₁ and λ₂ (the ± gives two). Everything — whether the car bounces or settles, whether the robot overshoots — is decided by one piece of this formula: the quantity under the square root, called the discriminant:

The sign of Δ is a fork in the road. It compares b² (damping, squared) against 4mk (spring stiffness scaled by mass). When damping wins, Δ > 0. When the spring wins, Δ < 0. When they exactly balance, Δ = 0. Those three cases — plus the special case b = 0 — are exactly the four personalities you felt in Chapter 0. Let's nail each one with real numbers.

Regime 1 — Undamped (b = 0)

Set m = 1, b = 0, k = 4. The characteristic equation is λ² + 4 = 0, so λ² = −4, giving λ = ±√(−4) = ±2i. The roots are purely imaginary (i = √−1). A purely imaginary root is the mathematical signature of pure oscillation that never decays. Using Euler's identity (e^iωt = cosωt + i·sinωt), the solution is x(t) = A·cos(2t) + B·sin(2t) — a sine wave of constant amplitude, ringing forever at angular frequency ω = 2 rad/s. This is the car with no shock absorber. No real part means no decay.

Regime 2 — Underdamped (Δ < 0, b small)

Now add a little damping: m = 1, b = 2, k = 5. Discriminant Δ = b² − 4mk = 4 − 20 = −16. Negative! So √Δ = √(−16) = 4i, and the roots are λ = (−2 ± 4i) / 2 = −1 ± 2i. Complex, with a negative real part (−1) and a nonzero imaginary part (±2).

Decode those two pieces. The real part −1 is a decay rate: it produces an envelope e^−1·t that shrinks the motion over time. The imaginary part 2 is an oscillation frequency: it produces sin and cos waves at ω = 2 rad/s. Together they give a decaying oscillation — bounces that shrink. The general underdamped solution is:

where here λ = 1 is the decay rate (the magnitude of the real part) and ω = 2 is the ringing frequency (the imaginary part). The constants A and B are fixed by where you start (initial position and velocity). This is the car with worn-but-not-dead shocks, the under-tuned robot arm that overshoots. Underdamped is the only regime that oscillates while still settling.

Regime 3 — Critically damped (Δ = 0)

Tune the damping until the discriminant hits exactly zero. With m = 1, k = 5, we need b² = 4mk = 20, so b = √20 ≈ 4.47. Then Δ = 0, and the ±√0 collapses both roots onto a single value: λ = −b/(2m) = −4.47/2 ≈ −2.24 (a repeated real root). No imaginary part, so no oscillation; the root is real and negative, so it decays. This is the critically damped case: the fastest possible return to rest without a single overshoot. It is the exact boundary — one drop less damping and it would overshoot; one drop more and it would crawl.

Critical damping is the holy grail of mechanical design: car suspensions, screen-door closers, scale needles, and well-tuned servo loops are all designed to sit right here. Fast, but clean.

Regime 4 — Overdamped (Δ > 0)

Pile on more damping: m = 1, b = 6, k = 5. Discriminant Δ = 36 − 20 = +16. Positive, so √16 = 4 is a real number, and the roots are λ = (−6 ± 4)/2 = {−1, −5} — two distinct real, negative roots. No imaginary part anywhere, so no oscillation. Two pure decaying exponentials: x(t) = A·e^−t + B·e^−5t. The motion is so thoroughly smothered by friction it can't even complete one swing. It just oozes home. Slower than critical — the price of over-insuring against overshoot.

The widget below is the regime explorer. As you drag b, watch the discriminant Δ flip sign, watch the two roots move in the complex plane (the right panel), and watch the resulting x(t) trajectory (the left panel). The whole story of this chapter lives in how those two dots move as you turn one knob.

Chapter 3: State-Space — Turning One Hard Equation Into Two Easy Ones

We can now solve the spring-damper by hand. But there's a problem lurking that you'll hit the moment you try to put this on a computer, or chain many of these together, or build a controller around it: the equation has a second derivative in it. m·x″ + b·x′ + k·x = 0 is a second-order ODE, and almost every tool in numerical computing, control theory, and machine learning is built for first-order systems — equations where the highest derivative is just a single prime.

So we need a way to rewrite one second-order equation as a system of first-order equations. The trick that does this is the single most important reframing in all of control theory and robotics, and it's called the state-space form. Once you see it you'll use it forever.

Step 1: Name the state variables

The state of a system is the minimal set of numbers you need to predict its entire future. For the spring-damper, knowing only the position is not enough — a block at x = 0 could be flying right, flying left, or sitting still, and those have wildly different futures. You also need the velocity. So the state is the pair (position, velocity). We give each a clean name:

This is the entire creative leap, and it's deceptively simple: we promoted velocity from a derivative to a first-class variable with its own name. Instead of one quantity x and its hidden derivatives, we now have two quantities y₁ and y₂ standing side by side. Everything that follows is just turning the crank.

Step 2: Differentiate each new variable

We want a first-order system, which means we need to express y₁′ and y₂′ — the rates of change of our two state variables — purely in terms of y₁ and y₂ (no second derivatives allowed). Let's compute them one at a time.

The first equation is almost free. What is y₁′? Since y₁ = x, its derivative is y₁′ = x′. But we just named x′ as y₂! So:

Read it as plain English: "the rate of change of position is the velocity." Obvious, almost insultingly so — and yet it is half of our system. The definition of velocity is itself one of our two equations. This is the gift the change of variables hands us for free.

The second equation needs the physics. What is y₂′? Since y₂ = x′, its derivative is y₂′ = x″ — the acceleration. And acceleration is exactly what the equation of motion tells us. Go back to m·x″ + b·x′ + k·x = 0 and solve it for x″:

We just isolated acceleration: it's a weighted combination of position and velocity. Now substitute our state names — x is y₁, and x′ is y₂:

And there is no second derivative left in sight — we expressed x″ entirely in terms of y₁ and y₂. We've successfully traded one second-order equation for two first-order ones. Let's set them side by side:

Step 3: Take m = 1 and write it as a matrix

To see the structure cleanly, set m = 1 (you can always rescale time later to recover any mass). The system becomes:

Notice I padded the first equation with "0·y₁ + 1·y₂" to make the pattern explicit: each y′ is a linear combination of y₁ and y₂. Whenever you see "outputs are linear combinations of inputs," you are looking at a matrix-times-vector in disguise. Stack y₁ and y₂ into a state vector y, collect the four coefficients into a 2×2 matrix, and the whole system collapses to one line:

Read off the matrix entries against the two scalar equations to confirm them. Row 1 is [0, 1], because y₁′ = 0·y₁ + 1·y₂. Row 2 is [−k, −b], because y₂′ = −k·y₁ − b·y₂. That matrix A is called the state matrix (or system matrix), and it is the single most important object in linear control theory. Every property of the system — stability, oscillation, speed of response — is encoded in those four numbers.

A fully worked change of variables with numbers

Let's grind one all the way through so nothing is abstract. Take the underdamped system from Chapter 2: m = 1, b = 2, k = 5. The state matrix is A = [[0, 1], [−5, −2]]. Now suppose the block is at state y = [y₁, y₂] = [0.4, −0.3] (40 cm right of rest, drifting left at 0.3 m/s). What is the instantaneous rate of change of the whole state, y′ = A·y?

Matrix-vector multiply, row by row. Row 1: y₁′ = 0×0.4 + 1×(−0.3) = −0.3. (Of course — the rate of change of position equals the velocity, −0.3, exactly as our top row promised.) Row 2: y₂′ = (−5)×0.4 + (−2)×(−0.3) = −2.0 + 0.6 = −1.4. So the velocity is changing at −1.4 m/s² — the block is accelerating leftward.

So y′ = [−0.3, −1.4]. In one matrix multiply we got both the velocity and the acceleration of the system at this instant. That single line, y′ = A·y, repeated over and over, is a numerical simulator for the spring-damper. The simulation you played with in Chapter 0 was secretly doing exactly this multiply, thousands of times per second.

From state-space to code (it's a five-line simulator)

Here is the entire spring-damper, simulated, in plain Python. Notice there is no special "ODE solver" magic — we just apply y′ = A·y and step forward a tiny bit in time (Euler's method). The matrix is the program.

python
import numpy as np

# Pick the system (m=1 baked in)
k, b = 5.0, 2.0
A = np.array([[0, 1],
              [-k, -b]])          # the state matrix

y  = np.array([1.0, 0.0])   # start: x=1, v=0 (pulled right, released)
dt = 0.01                       # tiny time step

for step in range(1000):
    ydot = A @ y                 # y' = A y   (the whole physics, one line)
    y    = y + ydot * dt         # Euler step: nudge state forward
    # y[0] is position now, y[1] is velocity now

Trace one iteration by hand to convince yourself it matches our arithmetic above. Start near y = [0.4, −0.3]. Then A @ y gives ydot = [−0.3, −1.4] — the exact numbers we computed. The new state is y + ydot·dt = [0.4, −0.3] + [−0.3, −1.4]×0.01 = [0.397, −0.314]. Position drifted slightly left, velocity got slightly more negative. Do that a thousand times and you've drawn the whole trajectory. That's it — that's how every physics engine and every robot simulator works under the hood.

Chapter 4: Eigenvalues — The Regimes Were Hiding in the Matrix

Here is the payoff that ties the whole lesson into a knot. In Chapter 2 we found that the system's personality is decided by the roots λ of the characteristic equation mλ² + bλ + k = 0. In Chapter 3 we packed the same system into a state matrix A = [[0, 1], [−k, −b]]. The stunning fact is: those roots λ are exactly the eigenvalues of A. The decay rates and oscillation frequencies you computed by hand are the eigenvalues of a 2×2 matrix. The two stories are one story.

Computing the eigenvalues by hand

Eigenvalues come from the condition det(A − λI) = 0, where I is the identity matrix. Let's just do it for our A and watch something beautiful fall out. First form A − λI by subtracting λ from the diagonal:

The determinant of a 2×2 matrix [[p, q], [r, s]] is ps − qr. Here p = −λ, q = 1, r = −k, s = −b−λ. So:

Expand the first product: (−λ)(−b−λ) = λb + λ² = λ² + bλ. And the second term: −(1)(−k) = +k. So the determinant is λ² + bλ + k. Setting it to zero:

Look at that. With m = 1, this is identically the characteristic equation from Chapter 2 (mλ² + bλ + k = 0 with m = 1). The "find the eigenvalues of the state matrix" calculation and the "guess e^λt and solve the quadratic" calculation produce the exact same equation. They were never two different things. The characteristic polynomial of the ODE is the characteristic polynomial of its state matrix.

Reading the four regimes off the eigenvalues

Now re-read Chapter 2's table through the eigenvalue lens. The eigenvalues live in the complex plane, and where they sit tells you everything:

The horizontal axis (real part) is the decay/growth axis: further left = faster settling, right of zero = explosion. The vertical axis (imaginary part) is the oscillation axis: further from the real axis = faster ringing. Damping pulls the eigenvalues down toward the real axis (killing oscillation) and the discriminant Δ = b² − 4mk is just the test for whether they've reached it yet.

A worked eigenvalue computation

Take the overdamped system m = 1, b = 6, k = 5, so A = [[0, 1], [−5, −6]]. The characteristic polynomial is λ² + 6λ + 5 = 0. Factor it: (λ + 1)(λ + 5) = 0, giving eigenvalues λ = −1 and λ = −5. Both real, both negative — so the system decays without oscillating (overdamped), exactly as Chapter 2 predicted. The −5 eigenvalue is the "fast mode" (decays as e^−5t, gone in a flash); the −1 is the "slow mode" (e^−t, the lingering tail). The overall settling speed is governed by the slowest eigenvalue, −1 — which is why overdamped systems are sluggish: one lazy eigenvalue drags out the whole response.

Now compare against critically damped (b = 2√5 ≈ 4.47): the polynomial is λ² + 4.47λ + 5 = 0 with the repeated root λ = −2.24. Both eigenvalues sit on top of each other at −2.24 — further left than the overdamped slow mode (−1)! That is the precise reason critical damping is faster than overdamping: it pushes the slowest eigenvalue as far left as possible before oscillation kicks in. Maximizing the worst-case decay rate — that's what "critically damped" means in eigenvalue language.

Where this goes next

You've now seen the full arc: physical forces → a second-order ODE → a quadratic and its four regimes → a first-order state-space system → eigenvalues of a matrix. That arc is the backbone of estimation and control. A few doors this opens:

In a Kalman filter, the system you're tracking is written in exactly this state-space form, y′ = A·y (often with added noise and control). The "state" is the same idea — position and velocity bundled together — and the matrix A propagates your belief forward in time. Everything you built here is the skeleton the Kalman filter hangs flesh on.

In control and reinforcement learning, designing a controller often means choosing where to place the eigenvalues of the closed-loop state matrix — literally dragging the dots in that complex-plane widget to a location that gives fast, non-oscillating behavior. That technique is called pole placement, and "poles" are just another word for these eigenvalues.

And the moment your system isn't linear (a real robot arm with gravity and joint limits), you'll linearize it around an operating point to get a local state matrix A, then read its eigenvalues to ask "is this configuration stable?" The spring-damper is the rehearsal for that move, performed exactly once and reused forever.

“What I cannot create, I do not understand.” You just created the spring-damper from three forces, solved it four ways, rebuilt it as a matrix, and found its soul in two eigenvalues. That is the whole method of physics, estimation, and control, rehearsed on the simplest thing that contains it all. — after R. Feynman