---
title: Numerical Dynamics
module: Dynamics
moduleNumber: 2
lessonNumber: 4
order: 204
summary: >
  Most force laws — quadratic drag, coupled oscillators, anything nonlinear — admit no
  closed-form trajectory, so we advance the motion one small time step at a time and
  let arithmetic do what algebra cannot. This lesson turns $\d\vec y/\d t=f(t,\vec y)$
  into a marching rule. We derive the Euler, Euler--Cromer, midpoint, and Verlet
  updates, weigh their accuracy and stability, watch a drifting energy expose a bad
  scheme, and use step-halving and conserved quantities to separate the error of the
  method from the error of the model.
topics: [Dynamics]
draft: false
sources:
  - book: Tipler & Mosca
    ref: "Ch. 5 — Additional Applications of Newton's Laws; §5-4"
---

## Initial-value models and Euler updates

Many force laws lead to differential equations whose exact solutions are unavailable
or unsuitable for the geometry and measured input data. A projectile with quadratic
drag, a spring with a measured nonlinear restoring force, and a vehicle whose thrust
changes with time all require the same numerical structure. At one time, record the
position $\vec r$, velocity $\vec v$, time $t$, mass $m$, and any
parameters entering the force law. The model yields acceleration through

$$
\vec a(t,\vec r,\vec v)=\frac{\vec F(t,\vec r,\vec v)}{m}.
$$

Initial values at $t_0$ determine one numerical trajectory once a time-step rule
is selected. Reducing the step size cannot compensate for missing initial conditions.
The launch position, initial velocity components, force parameters, coordinate
orientation, and stopping condition belong in the model statement before the first
update. A drag coefficient with units inconsistent with the chosen drag law can
generate a smooth-looking but physically meaningless trajectory.

Represent the state at sample $n$ by
$(t_n,\vec r_n,\vec v_n)$. A time interval $h$ advances the clock to
$t_{n+1}=t_n+h$. An update rule uses force information at the old state to estimate
the new position and velocity. A smaller interval usually improves the approximation,
but it also increases computation, roundoff accumulation, and the amount of measured
input that must be interpolated.

**Forward Euler updates.**

The forward Euler method uses the old acceleration throughout one short interval:

$$
\vec v_{n+1}=\vec v_n+\vec a_nh,
\qquad
\vec r_{n+1}=\vec r_n+\vec v_nh,
\qquad
\vec a_n=\frac{\vec F(t_n,\vec r_n,\vec v_n)}{m}.
$$

Velocity changes by acceleration times time, and position changes by the old velocity
times time. Every right-side quantity is known at step $n$, so the rule is explicit.
Calculate force, divide by mass, update velocity, and then update position. Updating
position with the new velocity defines a different method and must be recorded as
such.

For constant acceleration, forward Euler reproduces velocity exactly at the step
times but underestimates position during each positive-acceleration interval because
it uses the interval's starting velocity. The local position error scales as $h^2$;
after a fixed total time, the accumulated global position error scales approximately
as $h$. Halving $h$ should therefore roughly halve the error when Euler behavior
dominates and the solution remains smooth.

$$
% caption: Forward Euler follows the tangent velocity at the start of each time
% interval. Under positive acceleration the straight position advance lies below the
% curved exact position record, so repeated steps accumulate a systematic position
% error whose size depends on the chosen interval.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.30,.15)--(6.05,.15) node[right] {time};
  \draw[->,black] (.55,.00)--(.55,2.80) node[above] {position};
  \draw[acc,very thick] (.65,.30) .. controls (2.20,.45) and (4.20,1.45) .. (5.65,2.55);
  \draw[black,thick] (.65,.30)--(2.15,.44)--(3.65,1.02)--(5.15,2.02);
  \draw[black,dashed] (3.65,1.02)--(3.65,1.30);
  \draw[black,dashed] (5.15,2.02)--(5.15,2.42);
  \foreach \x/\y in {.65/.30,2.15/.44,3.65/1.02,5.15/2.02} {\draw[fill=white,draw=black,thick] (\x,\y) circle (1.6pt);}
  \node[acc] at (2.55,2.35) {smooth path};
  \node[black] at (4.10,.70) {Euler steps};
\end{tikzpicture}
$$

Velocity-dependent drag illustrates the need to evaluate the force at the current
state. In one dimension, a common model is

$$
m\frac{\d v}{\d t}=mg-cv|v|.
$$

With downward chosen positive, the drag term changes sign with velocity and
opposes motion. At every time step, calculate $v_n|v_n|$ from the signed old
velocity; replacing it by $v_n^2$ would apply downward drag even while the object
moves upward. A numerical state table should include time, position, velocity,
acceleration, force components, and any event flag such as ground contact. These
columns expose sign errors and prevent a plotted curve from becoming the only record
of the calculation.

> **Worked example.** A $0.50\ \mathrm{kg}$ ball falls from rest with quadratic drag
> $c=0.20\ \mathrm{kg\,m^{-1}}$ ($g=9.81\ \mathrm{m\,s^{-2}}$, downward positive), so
> $a=g-\tfrac{c}{m}v|v|$ with $c/m=0.40\ \mathrm{m^{-1}}$. The terminal speed is
> $$
> v_t=\sqrt{\frac{mg}{c}}
> =\sqrt{\frac{(0.50\ \mathrm{kg})(9.81\ \mathrm{m\,s^{-2}})}{0.20\ \mathrm{kg\,m^{-1}}}}
> =4.95\ \mathrm{m\,s^{-1}}.
> $$
> Take $h=0.10\ \mathrm s$ and $x_0=v_0=0$. Drag vanishes at $v_0=0$, so the first
> forward Euler step uses $a_0=9.81\ \mathrm{m\,s^{-2}}$:
> $$
> v_1=v_0+a_0h=0+(9.81)(0.10)=0.981\ \mathrm{m\,s^{-1}},\qquad
> x_1=x_0+v_0h=0.
> $$
> The position does not advance on the first step, because forward Euler moves it with
> the old velocity $v_0=0$. The second step evaluates drag at the new speed,
> $$
> a_1=g-\tfrac{c}{m}v_1|v_1|=9.81-(0.40)(0.981)^2=9.43\ \mathrm{m\,s^{-2}},
> $$
> $$
> v_2=0.981+(9.43)(0.10)=1.92\ \mathrm{m\,s^{-1}},\qquad
> x_2=0+(0.981)(0.10)=0.098\ \mathrm m.
> $$
> Drag has already cut the acceleration from $9.81$ to $9.43\ \mathrm{m\,s^{-2}}$; as
> $v\to v_t$ it falls to zero.

**Step size, convergence, and physical checks.**

Run the same model with several step sizes. Compare position, velocity, event time,
and conserved quantities at common physical times rather than comparing arrays with
different sample indices. For a ball with drag, compare its height and speed at a
fixed clock time and its computed ground-contact time. For an undamped oscillator,
compare amplitude and total energy. Report a result only when step reduction changes
the stated quantities by less than the required tolerance.

Euler's method can add or remove energy from an oscillator even when the modeled
force is conservative. The drift is numerical rather than physical. Under step
refinement, steady error reduction supports convergence; erratic changes can indicate
an event-handling mistake, an unstable step, or a discontinuous force law. Close
encounters, impacts, and sharp force transitions need smaller steps or an explicit
event solver.

$$
% caption: Step-refinement records compare the same physical time for several
% intervals. The smaller-step trajectory approaches a stable value, while a large
% interval can cross a contact boundary late or miss a short force pulse. Convergence
% is established from reported quantities, not from the visual smoothness of a plot.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.35,.35)--(6.05,.35) node[right] {time};
  \draw[->,black] (.55,.18)--(.55,2.72) node[above] {height};
  \draw[black,thin] (.70,.35) .. controls (1.60,2.30) and (3.30,2.30) .. (4.80,.35);
  \draw[acc,very thick] (.70,.35)--(1.30,1.55)--(1.95,2.16)--(2.70,2.24)--(3.45,1.80)--(4.15,1.02)--(4.80,.35);
  \draw[black,thick] (.70,.35)--(2.05,2.02)--(3.60,1.42)--(5.35,.35);
  \foreach \x/\y in {.70/.35,2.05/2.02,3.60/1.42,5.35/.35} {\draw[fill=white,draw=black,thick] (\x,\y) circle (1.7pt);}
  \draw[black,dashed] (4.80,.35)--(4.80,.12);
  \draw[black,dashed] (5.35,.35)--(5.35,.12);
  \node[black,below left] at (4.86,.14) {on time};
  \node[black,below right] at (5.30,.14) {late};
  \node[acc] at (3.05,2.45) {small step};
  \node[black] at (4.42,1.66) {large step};
\end{tikzpicture}
$$

Store enough information to reproduce the calculation: force equation, parameter
values and units, update order, step size, initial state, event conditions, software
precision, and convergence comparison. A numerical answer without those details is
an uncheckable trajectory, even when its final number has many displayed digits.

Use a convergence record with one row for each step size and one column for every
reported observable. Include the comparison time, event time when relevant, and the
difference from the next refinement. A final coordinate can appear stable while an
event time or phase still moves appreciably; recording both prevents an endpoint
agreement from masking a time-resolution failure.

**Forward Euler error at one step and many steps.**

Forward Euler replaces the exact state change over an interval by the tangent at the
interval start. Taylor expansion of a smooth position record gives

$$
x(t_n+h)=x_n+v_nh+\frac12a_nh^2+O(h^3).
$$

Euler keeps only the first two terms. The omitted position contribution over one
step is proportional to $h^2$. Velocity has the analogous omitted contribution from
the rate at which acceleration changes. These are local truncation errors: they
compare one numerical step with an exact step that begins from the exact state.

Over a fixed physical duration, the number of steps is proportional to $1/h$.
Accumulation then gives a global state error proportional to $h$ for a smooth
first-order Euler calculation. Halving the time step should approximately halve a
reported position or velocity error once the calculation is in its asymptotic
convergence range. This ratio need not appear at coarse steps, where the trajectory
may still be a poor approximation to the continuous solution.

Local and global error should be reported separately. A tiny local defect does not
guarantee a small final error after many steps, and a final error can partly cancel
for one special time without indicating a reliable method. Compare whole trajectories
or several stated sample times. For problems with a known analytic solution, measure
the numerical difference directly. For problems without one, use a successively
refined numerical result as a reference only after its own convergence has been
demonstrated.

The proportionality constant in an $O(h)$ global error depends on the force law,
initial state, and requested duration. Two first-order calculations can therefore
have very different errors at the same step size. The order predicts how error changes
under refinement; it does not predict whether a chosen step meets a physical
tolerance. A dimensional step check also catches inconsistent step units. The product of a characteristic
frequency and $h$ must be small enough to resolve the fastest smooth change in the
model. Reporting $h$ alone is incomplete when the force has a natural time scale.

## Euler--Cromer and convergence

Euler--Cromer changes one line of the update. Evaluate acceleration from the old
state, update velocity first, then use that new velocity for position:

$$
\vec v_{n+1}=\vec v_n+\vec a_nh,
\qquad
\vec r_{n+1}=\vec r_n+\vec v_{n+1}h.
$$

The method is still first order in global accuracy, but its long-time behavior can
differ sharply from forward Euler. Update order is part of the method definition.
Writing the two equations in the opposite order returns to forward Euler; averaging
the two velocities defines yet another scheme with different error properties.

> **Worked example.** Repeat the falling ball ($m=0.50\ \mathrm{kg}$,
> $c/m=0.40\ \mathrm{m^{-1}}$, $g=9.81\ \mathrm{m\,s^{-2}}$, $h=0.10\ \mathrm s$,
> $x_0=v_0=0$) with Euler--Cromer. The velocity step is identical,
> $v_1=0+(9.81)(0.10)=0.981\ \mathrm{m\,s^{-1}}$, but the position step uses the new
> velocity:
> $$
> x_1=x_0+v_1h=0+(0.981)(0.10)=0.098\ \mathrm m.
> $$
> Forward Euler gave $x_1=0$ from the same state. Continuing with
> $a_1=9.43\ \mathrm{m\,s^{-2}}$ and $v_2=1.92\ \mathrm{m\,s^{-1}}$,
> $$
> x_2=x_1+v_2h=0.098+(1.92)(0.10)=0.29\ \mathrm m,
> $$
> against $x_2=0.098\ \mathrm m$ for forward Euler. Both schemes are first-order in
> global error; the one-line change in update order shifts the position record from
> the first step onward.

For the harmonic oscillator, forward Euler expands the phase-space area slightly at
every step. The numerical orbit spirals outward and its calculated energy grows even
though the model force is conservative. Euler--Cromer has unit determinant for the
linear oscillator map. Its phase-space curve remains bounded for sufficiently small
$h\omega$, and its energy usually oscillates around the exact value instead of
drifting steadily upward. It does not reproduce the exact ellipse or conserve the
exact energy at each step.

$$
% caption: Phase-space behavior for a harmonic oscillator. Forward Euler produces an expanding spiral and secular energy growth, whereas Euler--Cromer keeps a bounded orbit whose energy varies about the physical value for a sufficiently small step.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.32,1.60)--(5.95,1.60) node[right] {position};
  \draw[->,black] (3.00,.20)--(3.00,3.10) node[above] {velocity};
  \draw[black,thick] plot[smooth] coordinates {(3.55,1.60)(3.43,2.03)(3.00,2.27)(2.48,2.12)(2.21,1.60)(2.40,1.00)(3.00,.69)(3.69,.91)(4.03,1.60)(3.77,2.37)(3.00,2.75)(2.14,2.46)(1.73,1.60)};
  \draw[acc,very thick] plot[smooth cycle] coordinates {(3.82,1.60)(3.58,2.18)(3.00,2.42)(2.42,2.18)(2.18,1.60)(2.42,1.02)(3.00,.78)(3.58,1.02)};
  \node[black] at (1.58,2.42) {Euler};
  \node[acc] at (4.42,.72) {Cromer};
\end{tikzpicture}
$$

The linear oscillator shows how update order reshapes phase-space structure; both
schemes stay explicit, and their accuracy is set by the step size. A rapidly changing
force, a high oscillator frequency, or a long simulation requires a refinement study.
Energy behavior is one diagnostic; position, velocity, and phase require separate
checks.

The linear oscillator's Euler--Cromer phase curve remains bounded only within
its stability range. Increasing $h\omega$ eventually produces a distorted numerical
orbit and then an unstable recurrence. Bounded energy at one coarse step is therefore
not enough evidence for accuracy. Compare the numerical period with the physical
period, inspect the phase offset after many cycles, and repeat with a smaller step.
The bounded phase-space structure separates a method with qualitatively appropriate
long-time behavior from one whose energy grows systematically, even when both have
first-order global state error.

> **Worked example.** For the oscillator $\ddot x=-\omega^2x$ with
> $\omega=2.0\ \mathrm{s^{-1}}$ ($m=1.0\ \mathrm{kg}$, $k=m\omega^2=4.0\ \mathrm{N\,m^{-1}}$),
> start at $x_0=1.0\ \mathrm m$, $v_0=0$, $h=0.10\ \mathrm s$. The exact energy is
> $E=\tfrac12kx_0^2=2.00\ \mathrm J$. Forward Euler
> ($x_{n+1}=x_n+v_nh$, $v_{n+1}=v_n+a_nh$, $a_n=-\omega^2x_n$) gives
> $$
> \begin{aligned}
> (x_1,v_1)&=(1.00,\,-0.40): & E_1&=\tfrac12(0.40)^2+\tfrac12(4)(1.00)^2=2.08\ \mathrm J,\\
> (x_2,v_2)&=(0.96,\,-0.80): & E_2&=\tfrac12(0.80)^2+\tfrac12(4)(0.96)^2=2.16\ \mathrm J.
> \end{aligned}
> $$
> Euler--Cromer from the same start ($x_{n+1}=x_n+v_{n+1}h$) gives
> $E_1=\tfrac12(0.40)^2+\tfrac12(4)(0.96)^2=1.92\ \mathrm J$ and $E_2=1.86\ \mathrm J$.
> Forward Euler's energy grows every step (the outward spiral above); Euler--Cromer's
> stays near $2.00\ \mathrm J$ and oscillates, returning over a full period
> $T=2\pi/\omega=3.14\ \mathrm s$.

**Convergence evidence.**

Use a fixed final time and compute the same state with steps $h$, $h/2$, and $h/4$.
Let $D_h$ be a norm of the difference between the $h$ and $h/2$ results at common
times. For a first-order method in its convergence range,

$$
D_{h/2}\simeq\frac12D_h,
\qquad
p\simeq\log_2\!\left(\frac{D_h}{D_{h/2}}\right)\simeq1.
$$

The norm can be an absolute position difference, a velocity difference, or a
combined state norm with declared units and scaling. Sampling at matching physical
times matters: comparing step number one hundred across runs compares different
times when the step sizes differ.

> **Worked example.** Integrate linear drag $\dot v=-\gamma v$ with
> $\gamma=1.0\ \mathrm{s^{-1}}$, $v_0=10\ \mathrm{m\,s^{-1}}$, to $T=1.0\ \mathrm s$.
> Forward Euler gives $v(T)=v_0(1-\gamma h)^{T/h}$; the exact value is
> $v_0e^{-\gamma T}=3.679\ \mathrm{m\,s^{-1}}$. With $E$ the absolute error in
> $\mathrm{m\,s^{-1}}$,
> $$
> \begin{aligned}
> h=0.500\ \mathrm s:&\quad v=10(0.500)^{2}=2.500,\quad E=1.179;\\
> h=0.250\ \mathrm s:&\quad v=10(0.750)^{4}=3.164,\quad E=0.515;\\
> h=0.125\ \mathrm s:&\quad v=10(0.875)^{8}=3.436,\quad E=0.243.
> \end{aligned}
> $$
> The error ratios $1.179/0.515=2.29$ and $0.515/0.243=2.12$ approach $2$ as $h$
> halves, giving observed orders $p=\log_2(2.29)=1.20$ and $\log_2(2.12)=1.08$. The
> global error is first-order in $h$, and the ratio settles toward $2$ only once the
> steps are fine enough to be in the asymptotic range.

$$
% caption: First-order convergence evidence. On logarithmic axes, error decreases with slope near one as the time step is reduced; a flattened or irregular sequence indicates that the selected steps have not reached a reliable refinement regime.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.45,.35)--(5.70,.35) node[right] {step size};
  \draw[->,black] (.45,.35)--(.45,2.95) node[above] {error};
  \draw[acc,very thick] (.80,2.45)--(4.95,.60);
  \foreach \x/\y in {1.05/2.34,1.95/1.94,2.90/1.51,3.85/1.10,4.75/.69} {\draw[fill=white,draw=black,thick] (\x,\y) circle (2pt);}
  \draw[black,dashed] (2.90,.35)--(2.90,1.51)--(.45,1.51);
  \node[acc] at (3.60,2.14) {slope one};
\end{tikzpicture}
$$

Convergence evidence should include the observable, final time, step sequence, and
the criterion used to stop refinement. A statement that two plots look similar does
not quantify numerical reliability. Report the change between the last two
refinements and compare it with the physical tolerance required by the calculation.
For conserved systems, also report the range of numerical energy over the interval;
a bounded oscillation and a steady drift can have similar endpoint values while
representing different integration behavior.

When no exact solution is available, self-convergence should use the same initial
state and force parameters in every run. Interpolate the finer record to the stated
comparison times only when those times do not already coincide. The interpolation
error must be smaller than the step-refinement difference being reported. A table of
step size, final position, final velocity, energy range, and adjacent-run difference
often communicates the evidence more clearly than several overlaid trajectories.
Convergence is a property of stated observables over a stated interval, not of a
single visually smooth curve.

## Events, impacts, and discontinuities

An ordinary update advances a smooth differential equation over one selected time
interval. A physical model can change before that interval ends. A particle can reach
a wall, cross a switching location, enter a constrained region, or pass a measurement
plane. These are events. Locate each event time as part of the numerical solution.

Represent an event by a scalar function $g(t,\vec r,\vec v)$. A contact with
a plane at $x=x_w$, for example, can use $g=x-x_w$. The event occurs at
$g=0$. After a trial step, compare the event-function values at the step endpoints.
Opposite signs bracket a crossing when the function is continuous over that segment.
A crossing can also be detected when a known monotonic coordinate reaches a stated
threshold. The bracket gives a time interval that contains the event, rather than
an assumption that the event occurred at the trial step endpoint.

Refine the bracket by bisection, by a safeguarded root method, or by reintegrating
substeps until the event-time tolerance is met. Then advance the state to that event
time, apply the contact or regime rule there, and integrate the unused remainder of
the original time interval with the new rule. Applying a contact response at the end
of a long step permits unphysical penetration and shifts every later state in time.
The error can remain visible even if the ordinary force update is otherwise accurate.

$$
% caption: Event bracketing during a position update. A trial step begins before a wall and ends beyond it, so the signed gap changes sign. Refining the interval locates contact before the state update is continued under the post-contact rule.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.35,.35)--(6.02,.35) node[right] {time};
  \draw[->,black] (.58,.18)--(.58,2.95) node[above] {gap};
  \draw[black,dashed] (.58,1.40)--(5.55,1.40);
  \node[black,right] at (5.55,1.40) {wall};
  \draw[black,very thick] (.85,2.45)--(4.65,.55);
  \draw[fill=white,draw=black,thick] (1.20,2.28) circle (1.8pt);
  \draw[fill=white,draw=black,thick] (4.30,.80) circle (1.8pt);
  \draw[black,dashed] (2.95,.35)--(2.95,2.20);
  \draw[fill=acc,draw=acc] (2.95,1.40) circle (1.8pt);
  \node[above] at (2.95,2.12) {contact};
  \draw[<->,black] (1.20,.74)--(4.30,.74);
  \node[above] at (2.75,.74) {bracket};
  \node[below] at (1.20,.35) {step start};
  \node[below] at (4.30,.35) {step end};
\end{tikzpicture}
$$

An event detector must state its direction convention. A contact surface can be
crossed from either side, while a detector may be intended to record only approach
from the allowed side. A near-tangent trajectory can touch $g=0$ without a sign
change. In that case, monitor a local minimum of the gap or use a detector based on
the normal velocity as well as the gap. The rule must match the physical boundary;
adding a sign-change test alone does not identify every possible contact.

**Impact and constrained contact states.**

At an ideal hard impact, position is continuous while normal velocity changes
abruptly. Let $\hat n$ point from the contact surface into the allowed
region and let $v_n^-=\vec v^-\mathbin{\cdot}\hat n$ be the incoming
normal velocity. A coefficient of restitution $e$ defines the post-impact normal
velocity by

$$
v_n^+=-e v_n^-,
\qquad 0\le e\le1.
$$

State the tangential velocity rule separately. A frictionless contact leaves the
tangential component unchanged; a contact model with friction requires an additional
impulse or force law. An impact model needs a normal direction, a restitution
convention, and a rule for every velocity component.

After locating contact, store the pre-impact state, apply the velocity update at the
event time, and continue from the same position. A small position offset into the
allowed region may be used to prevent repeated detection of the identical contact
because of roundoff. That offset must be recorded and must be much smaller than the
spatial resolution required by the problem. A large arbitrary offset changes flight
time and energy, especially in a small apparatus.

Retain the event bracket, the accepted contact time, and both one-sided velocities in
the output record. These values distinguish an integration error from a restitution
model choice when a rebound height or impulse is checked later. Repeated near-zero
events may indicate contact chatter, a tolerance mismatch, or a force law that needs
a persistent-contact regime rather than another isolated impact update.

For persistent contact, an impact rule is insufficient. A block resting on a support
requires a constraint force that prevents further motion into the support. The normal
force is determined together with the acceleration constraint; it is not a fixed
force inserted after penetration. Numerical contact methods often solve a complement
condition: either the gap is positive with zero normal force, or the gap is zero with
a nonnegative normal force. The active regime can change during the calculation.

$$
% caption: Contact-state update at a rigid boundary. The pre-contact velocity has a negative normal component toward the surface; the post-contact state reverses that component according to the stated restitution while position remains at the located contact point.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[black,very thick] (4.70,.35)--(4.70,3.05);
  \node[below] at (4.70,.35) {boundary};
  \draw[thick] (1.20,.95) circle (.12);
  \draw[thick] (4.52,2.05) circle (.12);
  \draw[->,acc,very thick] (1.40,1.02)--(4.28,1.96);
  \node[acc,below] at (2.55,1.22) {incoming};
  \node[below] at (1.20,.72) {pre-contact};
  \draw[->,black,very thick] (4.40,2.16)--(2.10,2.76);
  \node[black,above] at (2.80,2.72) {outgoing};
  \draw[->,black,thick] (4.52,2.05)--(3.55,2.05);
  \node[black,left] at (3.55,2.05) {normal};
\end{tikzpicture}
$$

**Discontinuous force regimes.**

A force law can be discontinuous even when no hard contact occurs. A control force
may switch at a position threshold, a material law may change at an interface, or a
constraint may become active when a gap closes. Treat each smooth regime with its own
force expression. Locate the boundary crossing, stop at the crossing, and restart
the integration with the new expression. Do not evaluate an average of incompatible
forces across a step unless the model explicitly defines such an average.

The usual smooth-solution error order does not automatically apply across a
discontinuity. Refinement can still improve event time and state accuracy, but the
observed ratios can be irregular until the event location is resolved. Report event
times and post-event states separately from ordinary smooth-interval errors. A
trajectory that looks continuous on a coarse plot can conceal a force switch displaced
by several time steps.

Velocity-dependent forces require the velocity at the state where the force is
evaluated. For quadratic drag in one dimension,

$$
F_{\rm drag}=-c\,v|v|,
$$

the sign reverses with velocity and the magnitude grows quadratically. Evaluating
drag from an outdated velocity gives an explicit approximation; evaluating it from a
new unknown velocity produces an implicit equation. Both choices are valid numerical
methods when stated, but they have different stability and error behavior. A force
table should record the time, position, velocity, regime identifier, and force
components used for each accepted step.

## Numerical reporting and reproducibility

Reproducible reporting begins with the mathematical state and ends with the accepted
output. State the coordinate definitions, initial position and velocity, mass,
parameters with units, force law, and every regime condition. Identify the integrator
and its update order, the nominal step size, adaptive-step limits if used, arithmetic
precision, and the interpolation method used for event roots. These details determine
the calculated path as directly as the numerical parameter values.

For each event, report the detector function, direction condition, time tolerance,
position tolerance, pre-event state, post-event rule, and any constraint-force
convention. A contact time rounded to a displayed digit can be insufficient if a
later measurement depends on phase or travel duration. Keep internal event times at
the solver precision and round only the final reported quantities.

Verification should include more than a plotted path. Compare event times after
tightening the root tolerance, compare post-event states after reducing the time step,
and check momentum or energy against the stated contact law. For velocity-dependent
forces, record the evaluated force alongside the state used to compute it. A result
can then be regenerated, inspected, and modified without guessing which discontinuous
rule was applied between two saved samples.

Event-time uncertainty can dominate a reported observable even when the ordinary
trajectory between events is accurate. Report the final bracket width or root residual
alongside the event time, and propagate that timing uncertainty into any speed, phase,
or duration derived from it. A detector tolerance smaller than the arithmetic noise
of the state evaluation gives no additional physical resolution. Conversely, a loose
event tolerance can mask improvements from a smaller integration step. The step size,
root tolerance, and state-interpolation rule must be varied separately when assessing
an event-driven calculation.

## Second-order and coupled updates

Higher-order methods estimate the force at a state closer to the middle of the time
interval. For the first-order state system

$$
\dot{\vec r}=\vec v,
\qquad
\dot{\vec v}=\vec a(t,\vec r,\vec v),
$$

the explicit midpoint method first forms a trial half-step from the old state:

$$
\begin{aligned}
\vec r_{n+1/2}&=\vec r_n+\frac h2\vec v_n,\\
\vec v_{n+1/2}&=\vec v_n+\frac h2\vec a_n,\\
\vec a_{n+1/2}&=\vec a(t_n+h/2,\vec r_{n+1/2},\vec v_{n+1/2}).
\end{aligned}
$$

The full update then uses the midpoint velocity and acceleration,

$$
\vec r_{n+1}=\vec r_n+h\vec v_{n+1/2},
\qquad
\vec v_{n+1}=\vec v_n+h\vec a_{n+1/2}.
$$

All components of the midpoint state must belong to the same intermediate time. A
force evaluated from a midpoint position and an old velocity is a different method.
For smooth forces, midpoint sampling captures the leading variation of acceleration
over the interval and has second-order global state accuracy. Its local truncation
error is proportional to $h^3$, while the accumulated error over a fixed duration is
proportional to $h^2$.

Midpoint is explicit because the half-step force uses quantities already computed
from the old state. It works directly with forces that depend on velocity, position,
and time. A drag force, a driven force, and a coupled spring force can all be
evaluated at the trial midpoint without solving an algebraic equation. The method
still assumes that the force remains smooth enough over the selected interval for one
midpoint sample to represent its variation.

$$
% caption: Midpoint construction for one coordinate. The old state advances to a trial half-step, where acceleration is evaluated before the full state update. The force sample is centered in time rather than held at the interval start.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.35,.35)--(5.92,.35) node[right] {time};
  \draw[->,black] (.60,.18)--(.60,3.00) node[above] {state};
  \draw[black,dashed] (1.05,.35)--(1.05,.78);
  \draw[black,dashed] (2.90,.35)--(2.90,1.60);
  \draw[black,dashed] (4.75,.35)--(4.75,2.44);
  \draw[->,black,thick] (1.15,.84)--(2.74,1.52);
  \draw[->,very thick] (3.02,1.66)--(4.60,2.36);
  \draw[fill=white,draw=black,thick] (1.05,.78) circle (2pt);
  \draw[fill=acc!14,draw=acc,thick] (2.90,1.60) circle (2.2pt);
  \draw[fill=white,draw=black,thick] (4.75,2.44) circle (2pt);
  \node[below] at (1.05,.35) {old state};
  \node[below] at (2.90,.35) {half step};
  \node[below] at (4.75,.35) {new state};
  \node[acc,above] at (2.98,1.76) {force sample};
\end{tikzpicture}
$$

**Velocity--Verlet for position forces.**

Velocity--Verlet is appropriate when force depends on position but not
velocity. Begin with $\vec a_n=\vec F(\vec r_n,t_n)/m$, advance position,
evaluate the new acceleration, and complete the velocity update:

$$
\begin{aligned}
\vec r_{n+1}&=\vec r_n+\vec v_nh+\frac12\vec a_nh^2,\\
\vec a_{n+1}&=\frac{\vec F(\vec r_{n+1},t_{n+1})}{m},\\
\vec v_{n+1}&=\vec v_n+\frac h2(\vec a_n+\vec a_{n+1}).
\end{aligned}
$$

The average acceleration in the last line is the update that makes velocity and
position consistent to second order. In a continuing calculation, the final
acceleration of one step becomes the initial acceleration of the next. The method
therefore needs one new force evaluation per subsequent step.

Velocity-dependent force requires a modified algorithm
because $\vec a_{n+1}$ also depends on the unknown $\vec v_{n+1}$. One option
is an implicit solve. Another is a method designed for the full first-order state
system, such as midpoint. Substituting the old velocity into a new-force evaluation
changes the stated method and its accuracy.

Velocity--Verlet has favorable long-time behavior for conservative position forces.
It is time reversible and keeps the phase-space structure close to the continuous
Hamiltonian flow. Its calculated energy is not exactly constant, but for a stable
step it commonly remains in a bounded oscillatory band rather than developing a
steady one-way drift.

> **Worked example.** Advance the oscillator ($\omega=2.0\ \mathrm{s^{-1}}$,
> $x_0=1.0\ \mathrm m$, $v_0=0$, $h=0.10\ \mathrm s$, $a=-\omega^2x$) by one
> velocity--Verlet step. With $a_0=-\omega^2x_0=-4.0\ \mathrm{m\,s^{-2}}$,
> $$
> x_1=x_0+v_0h+\tfrac12a_0h^2=1.0+0-\tfrac12(4.0)(0.10)^2=0.980\ \mathrm m,
> $$
> $$
> a_1=-\omega^2x_1=-3.92\ \mathrm{m\,s^{-2}},\qquad
> v_1=v_0+\tfrac12(a_0+a_1)h=\tfrac12(-7.92)(0.10)=-0.396\ \mathrm{m\,s^{-1}}.
> $$
> The exact solution $x=\cos\omega t$ gives $x(0.1)=0.98007\ \mathrm m$ and
> $v(0.1)=-0.39734\ \mathrm{m\,s^{-1}}$, so the step errors are
> $7\times10^{-5}\ \mathrm m$ and $1\times10^{-3}\ \mathrm{m\,s^{-1}}$, second-order and
> far below forward Euler's. The energy is
> $$
> E_1=\tfrac12(0.396)^2+\tfrac12(4)(0.980)^2=1.9992\ \mathrm J,
> $$
> a drift of $8\times10^{-4}\ \mathrm J$ ($0.04\%$) from $E_0=2.00\ \mathrm J$ that
> stays bounded, against forward Euler's $+4\%$ over the same step.

**Coupled coordinates at one time level.**

For several coordinates, write the state as vectors. With a mass matrix
$\vec M$,

$$
\vec M\ddot{\vec q}=\vec F(t,\vec q,\dot{\vec q}).
$$

Every force component must be evaluated from one consistent vector state. Consider
two masses coupled by springs:

$$
F_1=-k_1x_1-k_c(x_1-x_2),
\qquad
F_2=-k_2x_2-k_c(x_2-x_1).
$$

The coupling terms require both coordinates from the same time level. Updating
$x_1$ first and using that new value with old $x_2$ changes the force law during the
step and introduces an unintended ordering bias. Form the complete acceleration
vector, then apply the chosen update to all coordinates together.

$$
% caption: Coupled-coordinate force evaluation. Each spring force depends on the simultaneous positions of both masses, so one complete state vector is used before either coordinate is advanced to the next time level.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[black,very thick] (.28,.72)--(.90,.72);
  \draw[black,thick] (.90,.72)--(1.14,.94)--(1.38,.50)--(1.62,.94)--(1.86,.50)--(2.10,.72);
  \draw[black,fill=black!8,thick] (2.10,.30) rectangle (3.06,1.14);
  \draw[black,thick] (3.06,.72)--(3.32,.94)--(3.58,.50)--(3.84,.94)--(4.10,.50)--(4.36,.72);
  \draw[black,fill=black!8,thick] (4.36,.30) rectangle (5.32,1.14);
  \draw[->,black,thick] (2.58,1.44)--(2.58,2.14) node[above] {coordinate one};
  \draw[->,black,thick] (4.84,1.44)--(4.84,2.14) node[above] {coordinate two};
  \node[above] at (3.71,.94) {coupling};
  \node[below] at (2.58,.30) {mass one};
  \node[below] at (4.84,.30) {mass two};
\end{tikzpicture}
$$

Coupled systems have several natural frequencies. The highest resolved frequency
usually sets the restrictive step scale, even when the initial motion mainly excites
a slow collective coordinate. A weakly visible fast mode can accumulate phase error
over a long calculation. Inspect the normal-mode frequencies of a linearized system,
or estimate the fastest time scale from the force derivatives, before selecting a
nominal time step.

## Stability, adaptive steps, and error

Stability limits are method and model dependent. For velocity--Verlet applied to a
harmonic oscillator of angular frequency $\omega$, bounded discrete motion requires

$$
h\omega<2.
$$

The inequality marks the stability boundary. Accuracy requires a smaller step.
Near-boundary integration can remain bounded while producing substantial phase error
and a distorted energy oscillation. Choose $h\omega$ comfortably below the limit,
then assess period and state error at the physical duration of interest. In coupled
systems use the largest relevant frequency in this check.

$$
% caption: Stability and accuracy for an oscillator integration. Below the method-specific stability boundary the numerical orbit remains bounded, but phase error grows as the step approaches the boundary; beyond it, the discrete update is unstable.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.45,.40)--(5.95,.40) node[right] {step frequency product};
  \draw[->,black] (.68,.22)--(.68,2.96) node[above] {response};
  \draw[acc,very thick] (.95,1.15)--(3.95,1.15);
  \draw[acc,very thick] (3.95,1.15) .. controls (4.32,1.42) and (4.58,1.95) .. (4.82,2.72);
  \draw[black,dashed] (3.95,.40)--(3.95,2.62);
  \node[above] at (1.80,1.15) {bounded};
  \node[above] at (3.28,1.18) {phase error};
  \node[black,above] at (4.62,2.60) {unstable};
  \node[below] at (3.95,.40) {limit};
\end{tikzpicture}
$$

Dimension checks catch update mistakes before a trajectory is plotted. In
velocity--Verlet, $\vec v h$ and $\vec a h^2$ both have position units, while
$h(\vec a_n+\vec a_{n+1})$ has velocity units. In a coupled system,
$\vec M^{-1}\vec F$ must have coordinate acceleration units. Unit checks also
apply to nondimensional variables: state the length, time, and mass scales used to
recover physical results.

Exact limits provide compact validation. Under constant acceleration,
velocity--Verlet reproduces the exact position and velocity at all step times because
the acceleration at both ends is identical. With zero force, every method should
retain constant velocity exactly within arithmetic roundoff. For uncoupled
coordinates, a coupled-code implementation should reduce to independent
one-coordinate updates when the coupling constant is set to zero. These tests locate
indexing, force-evaluation, and unit errors without relying on an unknown reference
trajectory.

The simplest exact cases isolate distinct stages of the update. Their expected
records should be compared at identical physical times, not at array indices after
adaptive stepping.

| Limit | Exact record | Implementation defect exposed |
| --- | --- | --- |
| $\vec F=\vec 0$ | fixed velocity and linear position | time or state update |
| constant $\vec a$ | quadratic position and linear velocity | force centering or units |
| zero coupling | independent coordinate updates | coupled-force indexing |
| normal mode | analytic period and phase | mass matrix or force assembly |

Compare the numerical period, phase, and energy band of a harmonic oscillator with
the analytic solution over a stated duration. Compare a coupled linear system
against the normal-mode solution when available. A method can give a visually smooth
position trace while shifting phase by a significant fraction of a cycle. Exact-limit
tests and dimension checks therefore belong beside numerical output, together with
the integrator name, update equations, force model, step size, and parameter units.

For explicitly time-dependent forces, midpoint evaluation requires the midpoint time
as well as midpoint coordinates and velocities. A prescribed drive sampled at
$t_n$ instead of $t_n+h/2$ introduces a time-centering error even if the mechanical
state has been advanced correctly. Store the time associated with every force sample.
For tabulated force data, interpolate the drive to the time at which the numerical
method evaluates the state and record the interpolation rule.

Velocity--Verlet can be viewed as a position prediction followed by a force refresh
and a velocity correction. The new force must be computed after all coordinates have
been advanced to $\vec r_{n+1}$. Reusing the old force in the final velocity line
reduces the update to a different approximation. For a conservative force, monitor
the energy band, phase offset, and coordinate amplitude together. A narrow energy
band can coexist with a noticeable phase lag, so energy alone does not establish an
accurate time record.

Use complementary diagnostics for complementary numerical failures. A trajectory
can conserve a scalar quantity while locating the same physical event at the wrong
time.

| Diagnostic | What it constrains | What it can miss |
| --- | --- | --- |
| Energy band | conservative-update drift | phase offset |
| Phase difference | timing of oscillatory state | amplitude or energy error |
| Coordinate residual | sampled trajectory accuracy | force-model error between samples |
| Event-time residual | threshold or contact location | smooth-state accuracy away from the event |

Coupled coordinates also require consistent units in the mass matrix and force
vector. If one coordinate is an angle and another is a displacement, their entries
have different dimensions and the generalized forces are torques and forces,
respectively. Scaling coordinates before numerical integration can improve arithmetic
conditioning, but the scaling must be inverted before reporting physical amplitudes
or energy. The highest mode frequency should be recomputed after parameter changes,
because changing one coupling constant can alter the restrictive step scale for every
coordinate.

Validate in stages. Test zero force, constant acceleration, and one uncoupled
oscillator before the full coupled model. Excite one known normal mode at a time and
compare its period and phase with the analytic linear result. Then test a mode
superposition, where force evaluation must preserve the interaction terms. The staged
record locates errors in the integrator, coordinate transformation, coupled-force
expression, or parameter units.

Record the comparison interval and the norm used for each validation result. A maximum
coordinate error, an energy-band width, and a phase difference answer different
questions. Their numerical values cannot be compared unless the units, scales, and
sampling times are stated.

**Adaptive step control.**

A fixed time step spends the same computational effort in slowly varying and rapidly
varying parts of a trajectory. Adaptive integration changes the step size while
keeping a stated local error target. The solver proposes a step of size $h$, estimates
the numerical defect of that proposed update, and either accepts the state or repeats
the interval with a smaller step. Accepted states have unequal time spacing; a plot
must use their recorded times rather than their array indices.

A common estimator compares two approximations of different order or compares one
full trial step with two half steps. Let the resulting position and velocity
differences be $\delta\vec r$ and $\delta\vec v$. Normalize them with absolute
and relative tolerances:

$$
e=
\sqrt{
\left(\frac{|\delta\vec r|}
{r_{\rm abs}+r_{\rm rel}|\vec r|}\right)^2+
\left(\frac{|\delta\vec v|}
{v_{\rm abs}+v_{\rm rel}|\vec v|}\right)^2
}.
$$

A step with $e\leq1$ meets the selected state tolerance. A step with $e>1$ is
discarded and recomputed from the unchanged initial state with a shorter interval.
The normalization prevents a coordinate with large numerical units from hiding an
error in a smaller coordinate. Position and velocity tolerances should reflect the
physical quantity being reported. Available floating-point digits do not define the
required physical tolerance.

An estimator whose local defect scales as $h^{p+1}$ can use a trial next step:

$$
h_{\rm next}=s\,h\,e^{-1/(p+1)},
$$

where $s<1$ is a safety factor. Limit the growth and reduction factors so a single
small estimate does not produce an excessively large next step. A maximum step
protects output resolution and force sampling. A minimum step detects a demand for
more resolution than the specified precision or model can support. Repeated rejected
steps at the minimum should be reported as a tolerance or model failure, not silently
converted into an accepted inaccurate state.

$$
% caption: Adaptive time grid along a varying trajectory. Short accepted intervals cluster where the state bends rapidly, while longer intervals are used where the local motion changes slowly; each accepted state retains its own recorded time.
\begin{tikzpicture}[>=stealth,font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[->,black] (.35,.35)--(6.02,.35) node[right] {time};
  \draw[->,black] (.60,.18)--(.60,2.98) node[above] {position};
  \draw[acc,very thick] (.76,.55) .. controls (1.85,.60) and (2.45,.90) .. (3.10,1.80) .. controls (3.46,2.40) and (4.30,2.62) .. (5.55,2.74);
  \foreach \x/\y in {.98/.56,1.86/.68,2.52/1.06,2.86/1.44,3.10/1.80,3.36/2.18,3.72/2.44,4.40/2.60,5.30/2.71} {\draw[fill=white,draw=black,thick] (\x,\y) circle (1.8pt);}
  \draw[<->,black] (.98,.82)--(1.86,.82);
  \node[above] at (1.42,.84) {large step};
  \node[acc] at (4.42,1.78) {dense samples};
\end{tikzpicture}
$$

Adaptive control estimates numerical truncation error. A force coefficient with
uncertain units, an omitted force, or a measurement bias can still produce a smooth
trajectory with a small step estimator. Keep step selection separate from physical
parameter uncertainty. Use adaptive stepping to resolve the stated differential
equation, then assess whether that equation represents the apparatus.

Accepted adaptive states are solver states, not automatically the desired measurement
times. Interpolate or request dense output with the method's stated order when a
camera frame, sensor timestamp, or comparison time falls between accepted steps.
Store that interpolation procedure with the data because its error contributes to
the reported observable independently of the accepted-step tolerance.

**Numerical error budgets.**

A reported numerical uncertainty has several sources. Truncation error arises from
approximating continuous motion over finite intervals. Floating-point roundoff enters
through arithmetic and subtraction of nearly equal values. Interpolation error enters
when states are requested at times not retained by the solver. Parameter uncertainty
comes from masses, force constants, initial conditions, and measured driving data.
Model discrepancy comes from a physical approximation rather than the numerical
algorithm.

Numerical and model errors require separate estimates. Reducing
the step size can reduce truncation error while leaving parameter uncertainty
unchanged. Raising arithmetic precision can reduce roundoff while leaving an
inaccurate force coefficient unchanged. A numerical report should state the estimated
integration contribution separately from the uncertainty assigned to physical inputs.

A conservative budget may add independent upper bounds. When the individual terms are
well characterized as independent random uncertainties, a root-sum-square estimate is
often appropriate:

$$
\sigma_{\rm total}\simeq
\sqrt{\sigma_{\rm num}^2+\sigma_{\rm par}^2+\sigma_{\rm meas}^2}.
$$

The notation is meaningful only when each term refers to the same observable at the
same physical time. An energy uncertainty and a position uncertainty cannot be added
directly. Propagate each source to a declared output, such as final height, period, or
arrival time, before forming a budget.

Tolerance selection should begin with the required output precision. If an experiment
measures a displacement to one millimetre, requesting nanometre numerical accuracy
rarely changes the physical conclusion. Conversely, a timing calculation that
distinguishes two nearby phases may require a much tighter state tolerance than a
position plot suggests. Use separate tolerances for coordinates with different scales
and state the criterion used for each accepted step.

## State logs and validation

A state log must allow another calculation to reconstruct both the trajectory and the
solver decisions. For every accepted state, record time, all coordinates, all
velocities, the step used to reach that state, the force or acceleration evaluated
there, and the estimated local error. Record rejected trials separately with their
proposed step and error estimate. Omitting rejected trials hides whether the solver
encountered regions that demanded much smaller intervals.

The log also needs a machine-readable description of the model: integrator name and
version, absolute and relative tolerances, maximum and minimum step, arithmetic
precision, parameter values with units, initial state, and force-function revision.
A plotted curve is a derived product of this information. It does not record a
changed unit, a different tolerance, or a force evaluation made with stale state
variables.

Use a fixed column order and include units in metadata rather than appending unit text
to every numerical cell. A row identifier and monotonically increasing time simplify
audits. If output is decimated for plotting, retain the complete accepted-state log
separately. Decimation can erase the short time intervals that demonstrate why an
adaptive solver used additional work in one portion of the trajectory.

**Analytic and limiting-case validation.**

Validate the complete numerical pipeline against cases with known answers. With zero
force, position must advance linearly and velocity must remain fixed. With constant
acceleration, position must follow a quadratic time record. For linear drag,
$\dot v=-\gamma v$, compare with the exponential solution over several decay times.
The tests check the state update, force interface, units, and time labels without
requiring an elaborate apparatus model.

A limiting case can isolate one parameter. Let a coupling constant tend to zero and
verify that coordinates evolve independently. Let a forcing amplitude tend to zero
and verify recovery of the unforced solution. Reduce a drag coefficient and verify
approach to the corresponding conservative trajectory over a declared interval. The
comparison should use the same initial state and output times; changing multiple
parameters at once obscures the source of a discrepancy.

Validation also needs a reference for scale. Check that each term in an acceleration
has units of length divided by time squared, that each logged energy has consistent
units, and that the reported time span covers the physical behavior of interest.
Compare analytic residuals, limiting-case residuals, and error-budget terms at the
same stated output points. Agreement in one case does not validate every force regime,
but a failure in a simple exact limit identifies a defect before a more complicated
calculation is interpreted.

Use a validation ledger that preserves the exact limit, the observable compared, and
the time points at which the residual was evaluated. A small final-position error
does not substitute for a phase or event-time check.

| Test case | Expected relation | Residual to retain |
| --- | --- | --- |
| Zero force | constant $\vec v$ and linear $\vec r(t)$ | position and velocity drift |
| Constant acceleration | quadratic position record | error at common output times |
| Linear oscillator | period and phase from the analytic solution | phase offset and energy band |
| Event or impact | stated bracket and transition rule | event-time and post-event state error |

Repeatability is a separate requirement from agreement with an analytic case. Given
the same input file, arithmetic mode, and force data, a deterministic run should
produce the same accepted times and state log. If parallel evaluation or randomized
sampling is used, record the execution settings and any random seed. A change in the
accepted-step sequence can alter interpolation times and downstream summaries even
when a final plotted curve appears unchanged. Retain the raw log used to create a
reported table or figure.

Reference comparisons should include both absolute and scaled residuals. An absolute
position difference is appropriate near a known fixed origin; a relative difference
can be more informative when the expected state spans several orders of magnitude.
Do not form a relative residual by dividing by an exact value that crosses zero.
State the alternative scale used near that point. These conventions make validation
results interpretable across trajectories with different coordinate ranges.

Archive the input parameters, accepted-state log, validation residuals, and plotting
script together. A numerical claim remains verifiable only when those records identify
the exact calculation that produced it.

Report a numerical trajectory at physically meaningful times and events. A landing
time, turning point, peak speed, or maximum load requires an interpolation rule and
an uncertainty that includes step control and model parameters. A final array entry
is rarely the desired physical quantity by itself. State whether an extremum was
sampled at an accepted state, interpolated between accepted states, or located by an
event search. Those methods can agree to the displayed precision in a well-resolved
calculation, but they have different error behavior near sharp force changes.

Keep numerical precision separate from physical accuracy. A solver may store fifteen
decimal digits while uncertainty in a drag coefficient, launch angle, or contact law
limits a predicted range to three significant figures. Round the reported result only
after convergence and parameter sensitivity have been assessed. Give the method,
step-control settings, validation case, and uncertainty source with the result. That
record connects a calculated number to the force model and measurement assumptions
that determine its scientific meaning.