Numerical Dynamics

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

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:

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.

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

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.

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.

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

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:

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.

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.

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.

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,

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.

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.

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}}^{-}\cdot \hat{n}$ be the incoming normal velocity. A coefficient of restitution $e$ defines the post-impact normal velocity by

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.

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,

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

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

The full update then uses the midpoint velocity and acceleration,

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.

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:

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.

Coupled coordinates at one time level.

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

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

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.

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

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.

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.

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.

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:

A step with $e\le 1$ 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:

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.

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:

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.

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.