W05: From differential equations to trustworthy simulation

Prerequisites and observable goals#

The central question: what is the simulation supposed to reproduce?#

A simulation consists of at least three objects: a mathematical model, a numerical method, and a particular implementation. The model specifies a continuous law under assumptions. The method approximates that law using finite steps. The implementation stores finite-precision state and executes updates. A defect in any layer can give wrong results, and testing one layer does not automatically validate the others.

Use a unit-mass harmonic oscillator as a bounded example. Position is $x(t)$, velocity $v(t)$, and angular frequency $\omega>0$. The equations are

The exact solution is $x(t)=\cos(\omega t)$ and $v(t)=-\omega\sin(\omega t)$ in the chosen position units. Position and velocity do not have interchangeable units; $\omega$ has units inverse time, so $\omega^2x$ is an acceleration. This model idealizes a restoring force proportional to displacement, with no damping, driving force, collision, or changing parameter.

Define energy per unit mass $E=\tfrac12(v^2+\omega^2x^2)$. Differentiating along the exact solution gives

Conservation follows from the specified equations. It does not prove that an arbitrary animation or numerical integrator conserves that energy.

Explicit Euler: a reasonable local approximation with a global weakness#

Using the old state for both right sides yields

For small $h$, each resembles the first terms of a Taylor expansion. But substitute both into the discrete energy:

The cross terms cancel; the remaining squared-step terms are positive. Thus $E_{n+1}=(1+h^2\omega^2)E_n$. For every nonzero positive step, this method gains energy on a nonzero state. Reducing $h$ reduces the drift over a fixed interval, but no positive fixed step makes this scheme energy-stable for arbitrarily many oscillations. “The force points inward” is not enough to prevent numerical growth.

This is not a claim that Euler is useless. It is a precise warning about one model and long-time behavior. For a decaying scalar equation $y'=-ky$, its stability condition is different. Analyze the differential equation and the update matrix actually being used.

Velocity Verlet: split the acceleration update#

An acceleration depending on position suggests a half velocity update, a position update, then the remaining half velocity update:

The second half uses the new position. Store the intermediate explicitly; casually overwriting values may define another method. Expanding gives

This matrix has determinant 1 and trace $2-q^2$. Its eigenvalues lie on the unit circle with distinct oscillatory modes for $0<q<2$. For $q>2$, a real eigenvalue has magnitude greater than 1 and some states grow. At $q=2$, a repeated eigenvalue can have a nontrivial Jordan block, so the endpoint does not supply a safe bounded-oscillation guarantee. Use the strict interval.

For $0<q<2$, this scheme exactly preserves, in real arithmetic, the positive quadratic expression

It generally does not preserve the original $E$ exactly. Because $\widetilde E$ is positive definite in this range, it bounds $x$ and $v$ over arbitrarily many exact-arithmetic steps. Floating-point roundoff and changing step sizes introduce further issues. The distinction between a modified invariant and exact physical energy is important, not a technicality.

A small complete step#

With $\omega=1,h=.1,x_0=1,v_0=0$, the half velocity is $-.05$. Position becomes $.995$. Final velocity is $-.05-.05(.995)=-.09975$. Compare with exact $\cos(.1)\approx.9950041653$ and $-\sin(.1)\approx-.0998334166$. These errors are small but nonzero. The computed physical energy is

not exactly $.5$. The modified energy is $.49875$ both before and after the step in exact rational arithmetic. Check the equations rather than describing any tiny numerical change as “energy conservation.”

Accuracy, stability, and model adequacy are separate#

A stable method can still have a phase error: the oscillation stays bounded but peaks at the wrong times. Under the Verlet stability relation, the numerical phase angle satisfies $\cos\theta=1-q^2/2$. Its effective frequency $\theta/h$ differs from $\omega$ for nonzero step. Energy bounds alone do not guarantee accurate event times or positions.

For a fixed final time $T$, compute with steps $h,h/2,h/4$ and compare all methods at that same time. A second-order method in its asymptotic regime often has error ratios approaching 4, because the global error scales like $Ch^2$. A ratio is evidence under a model, not proof of the order for every input. Roundoff, a near-zero leading error coefficient, instability, or a changed final time can spoil the naive interpretation. Use known analytic solutions where available; use independent methods and convergence studies for harder cases.

Even a highly accurate solution of this oscillator model would not validate a real suspension, vehicle, robot, or orbital system. The physical model may lack damping, nonlinearity, external forces, measurement uncertainty, and safety constraints. The lab is a reproducible numerical experiment, not safety-critical controller software. Later T04/T05 and GP2 reuse the same discipline of explicit assumptions and fault boundaries in other domains.

Independent practice#

F foundation C standard course work T transfer S optional stretch L runnable lab or project. Attempt the complete set before opening its solutions.

Complete solutions: W05#

Compare the first step where your work differs. A corrected answer is not yet an independent success.

Open every solution

"""W05: velocity Verlet for x''=-omega^2*x, unit mass, fixed time step."""
import math


def verlet(x=1.,v=0.,omega=1.,dt=.1,steps=100):
    if not all(math.isfinite(z) for z in (x,v,omega,dt)) or omega<=0 or dt<=0:
        raise ValueError("finite state, positive frequency and step required")
    if type(steps) is not int or not 0<=steps<=100000:
        raise ValueError("bounded nonnegative step count required")
    energy0=.5*(v*v+omega*omega*x*x)
    if not math.isfinite(energy0):
        raise ValueError("initial energy is not finite; rescale the problem")
    rows=[(0.,x,v,energy0)]
    for k in range(1,steps+1):
        half=v-.5*dt*omega*omega*x
        new_x=x+dt*half
        new_v=half-.5*dt*omega*omega*new_x
        x,v=new_x,new_v
        energy=.5*(v*v+omega*omega*x*x)
        if not all(math.isfinite(z) for z in (x,v,energy,k*dt)):
            raise ValueError("nonfinite trajectory; check step and scale")
        rows.append((k*dt,x,v,energy))
    return rows


def exact(t,x0=1.,v0=0.,omega=1.):
    if not all(math.isfinite(z) for z in (t,x0,v0,omega)) or omega<=0:
        raise ValueError("finite inputs and positive frequency required")
    phase=omega*t
    if not math.isfinite(phase):
        raise ValueError("nonfinite phase; rescale the problem")
    result=(x0*math.cos(phase)+(v0/omega)*math.sin(phase),
            -omega*x0*math.sin(phase)+v0*math.cos(phase))
    if not all(math.isfinite(z) for z in result):
        raise ValueError("nonfinite analytic state; rescale the problem")
    return result


if __name__=='__main__':
    for dt in (.1,.05,.025):
        t,x,v,e=verlet(dt=dt,steps=round(1/dt))[-1]
        target=exact(t)
        print(dt,math.hypot(x-target[0],v-target[1]),e)

Retention and readiness#

Explain why a locally plausible derivative step can create global energy growth. Name the conserved quantity actually checked. Before trusting a simulation, state model error, discretization error, and finite-precision error separately.