Integration: Forces into Motion

in reviewThis lesson is complete and being validated before release - details may still change.

Integration

Loading interactive demo…

One particle on a spring to the anchor, under drag and gravity if they are on, stepped by the chosen integrator. Each step measures the forces where the particle is, kicks the velocity by the acceleration times dt, and moves along a velocity for dt; the method decides which happens first and which velocity the move uses. The exact motion is known in closed form, so every step is audited against it.Take the tour, or press and explore. Legend:

The panel beside the stage is the real ParticleIntegrator class from the published library executing in your browser, with the current line highlighted.

The particle, its path so far, and the ghost of where the run ends.
The exact motion, dashed, and the hollow point where the particle should be now.
The velocity.
The spring's pull, back toward the anchor.
Gravity's pull, when it is on.
Drag, against the motion, when it is on; the spring itself, dashed to the anchor.
The sum of the pulls: the acceleration.

time is the seconds simulated, energy the particle's kinetic plus potential energy, drift how far that has strayed from the exact motion's energy as a share of the energy the run began with, error the distance from the particle to its exact place, and speed its speed.

The code panel shows bare executable code by default. The bar above it names the file (hover the name for the full published path); options holds the toggles: comments, the tracer instrumentation that produces this visualization, and verbose names for the terse identifiers (display only: the default is the byte-exact published file). The view switch adds plain-English narration: both interleaves it with the code, pseudo shows only the pseudocode. hide code collapses the panel; the rail it leaves behind brings it back.

Keys: Space play/pause · ←/→ step · Home/End jump · R reset · F fullscreen · G guided tour

The problem. Newton says $a = F / m$ : the forces tell you the acceleration, never the position. A computer can only take steps. With a particle at $p$ moving at $v$ , and $a(p, v)$ the sum of the spring $-k p$ , the drag $-c v$ , and gravity $-g$ , a step of $\Delta t$ must produce a new $p$ and a new $v$ . The question is in which order, and with which velocity.

Explicit Euler. Move first, along the velocity you have, then kick: $p' = p + v\,\Delta t$ , then $v' = v + a\,\Delta t$ . On a spring with $\omega = \sqrt{k}$ the step multiplies the amplitude by $\sqrt{1 + \omega^2 \Delta t^2}$ , which is always more than one. Every step adds energy, so the orbit spirals outward whatever $\Delta t$ you choose; a smaller step only slows the spiral. The drift readout is the energy invented so far, as a share of what the run began with.

Semi-implicit Euler. Kick first, then move along the new velocity: $v' = v + a\,\Delta t$ , then $p' = p + v'\,\Delta t$ . The same arithmetic in the other order, and the amplitude factor becomes exactly one: the step is symplectic, it preserves area in the $(p, v)$ plane, and the energy oscillates about its true value without ever creeping. The orbit is a little off in phase and a little squashed, but it is an orbit, for as long as you care to run it, provided $\omega \Delta t < 2$ .

Velocity Verlet. Take half the kick, move, measure the force at the new place, and finish the kick with the mean of the two accelerations: $p' = p + v\,\Delta t + \tfrac{1}{2} a\,\Delta t^2$ , then $v' = v + \tfrac{1}{2}(a + a')\,\Delta t$ . It is second order, so halving $\Delta t$ quarters the error instead of halving it, and it is symmetric in time, which is why the energy stays put. It costs a second force evaluation per step, which is why the trace shows two force events.

Order of accuracy. After a fixed simulated time, explicit and semi-implicit Euler are off by an amount proportional to $\Delta t$ ; Verlet by an amount proportional to $\Delta t^2$ . Set $\Delta t$ to $0.05$ and then $0.5$ and compare the error readout for each method: the Euler pair grow tenfold, Verlet a hundredfold, and at the coarse step all three lose the curve, Verlet least.

Stability. A method is stable when a step cannot grow an oscillation. Explicit Euler never is. The symplectic pair hold while $\omega \Delta t < 2$ : with stiffness $10$ , $\omega \approx 3.16$ , and $\Delta t = 0.5$ sits just inside the limit, so the orbit holds though it is a poor shape. Stiffer springs need smaller steps, which is why a cloth with stiff constraints is such hard work for an explicit method.

The exact motion. The forces here are linear, so each axis obeys $u'' + c u' + k u = f$ and the solution is a damped cosine about the rest point $f / k$ : the origin across, and $-g / k$ down once gravity is on. The stage draws it dashed, and the error readout is the distance from the particle to where that motion puts it now. With drag on, the exact energy decays too; drift is still measured against it, so a method that loses energy at the right rate reads as zero.

Cost. Every method is $O(1)$ per step; Verlet pays two force evaluations where the others pay one. In a game the forces are the expensive part, so engines favour the semi-implicit step: one force, one kick, one move, and an orbit that never explodes.

Choices worth reading closely in the implementation (the code panel's comments and tracer checkboxes reveal the full published file):

The three methods share two private moves, push for the kick and move for the step along a velocity, so each tag has one call site and the trace shows the method's order rather than three copies of the arithmetic. Explicit Euler calls move with the velocity the step began with; semi-implicit calls it with the velocity push just returned.
Verlet calls accelerationAt twice, once at the old place and once at the new, then hands push the mean of the two: that is why a Verlet step traces five events and the others four.
The forces are a static, forces(state, options), so the stage draws exactly the pulls the step summed; energy is a static for the same reason.
exact is the closed form, a damped cosine per axis measured from a reference state, and retune moves that reference to the state it finds: a live run keeps a truth to compare against after every knob change.
The drag is evaluated at the half-step velocity in Verlet's second force measurement; with drag off, which is the default, the two measurements differ only by the spring at the new place.
Validation rejects a drag heavier than the spring can damp, so the underdamped closed form is always the right one.

Euler gave the method that bears his name in 1768, in his Institutionum calculi integralis: when an equation says how fast a quantity changes, step forward a little along that rate and repeat. Newton had drawn the picture in the Principia, building an orbit from a chain of straight pushes, and the scheme Verlet made standard for molecular dynamics in 1967 had served Delambre for comets in the 1790s and Stormer for the aurora in 1907. Games met the question when every frame had to move a world: which cheap step keeps a spring, a ragdoll, or an orbit from exploding?

1687Isaac Newton - In the Principia's first proposition, builds an orbit under a central force from a chain of straight-line pushes, the picture every numerical integrator still draws.source
1768Leonhard Euler - Publishes the Institutionum calculi integralis, with the method of stepping a differential equation forward along its current rate of change.source
1907Carl Stormer - Computes the paths of charged particles in the Earth's magnetic field with the central-difference scheme later named for Verlet, to explain the aurora.source
1967Loup Verlet - Uses the scheme for molecular dynamics of argon in Computer Experiments on Classical Fluids, and it becomes the integrator of choice for particle simulation.source
1981Alan Cromer - Publishes Stable Solutions Using the Euler Approximation, crediting a student with the observation that updating the velocity first keeps an oscillator's energy bounded: the semi-implicit, or Euler-Cromer, method.source
2001Thomas Jakobsen - Presents Advanced Character Physics at the Game Developers Conference: Verlet integration with constraints as the ragdoll and cloth engine of Hitman: Codename 47.source

Used today

Game physics engines - Rigid-body engines step velocities first and then positions, the semi-implicit order, because a game cannot afford a spring that gains energy every frame.
Ragdolls and cloth - Position-based Verlet with distance constraints moves character limbs and cloth, since a constraint is just a position to nudge back.
Molecular dynamics - Velocity Verlet steps millions of atoms for nanoseconds of simulated time, its time symmetry keeping the total energy honest over billions of steps.
Orbits and spacecraft - Symplectic integrators let the solar system be run for millions of years without the planets drifting outward through invented energy.

method

dt0.10 s

stiffness4.0

drag0.0

gravity

steps

show

time: 5.0 s
energy: 184.77
drift: +610.7%
error: 4.85
speed: 11.83

particle-integrator.ts

import { Vec2 } from "../linear/vec";import type { Tracer } from "../trace";/** How a step turns the acceleration into a new velocity and position. */export type IntegrationMethod = "euler" | "semi-implicit" | "verlet";/** Where the particle is and how it is moving, in world units and units per second. */export interface ParticleState {  position: Vec2;  velocity: Vec2;}export interface IntegratorOptions {  method: IntegrationMethod;  /** How many seconds one step advances the clock. */  dt: number;  timeStep: number;  /** The spring to the origin: a pull of -k p. */  stiffness: number;  /** Linear drag: a pull of -c v against the motion (default 0, none). */  drag?: number;  /** The downward pull, in units per second squared (default 0, none). */  gravity?: number;}/** The pulls on the particle at a state, per unit mass, so each is an acceleration. */export interface ParticleForces {  spring: Vec2;  drag: Vec2;  gravity: Vec2;  total: Vec2;}export interface ParticleEnergy {  kinetic: number;  potential: number;  total: number;}/** * One particle of unit mass on a spring to the origin, under optional * drag and gravity, stepped by one of three integrators. Explicit Euler * moves along the old velocity and then updates it, so every step adds * energy and an orbit spirals out. Semi-implicit (symplectic) Euler * updates the velocity first and moves along the new one: the same * arithmetic in the other order, and the energy stays bounded. * Velocity Verlet takes half the kick, moves, measures the force at the * new place, and finishes the kick with the mean of the two: second * order, and symmetric in time. The forces are linear, so the exact * motion is a damped oscillation about the rest point on each axis, and * `exact(t)` gives it in closed form for the error and drift readouts; * `retune` measures it from the state the retune found, so a live run * keeps a truth to compare against. Every step traces the force, the * velocity update, the move, and an energy audit, in the order the * method does them. * * Time: O(1) per step. Space: O(1). */export class ParticleIntegrator {  private current: ParticleState;  private clock = 0;  private tick = 0;  private optionsValue: IntegratorOptions;  /** The state and time the exact solution is measured from: the start, or the last retune. */  private reference: { state: ParticleState; time: number };  constructor(    initial: ParticleState,    options: IntegratorOptions,    private readonly tracer?: Tracer,  ) {    ParticleIntegrator.validate(options);    this.optionsValue = options;    this.current = { position: initial.position, velocity: initial.velocity };    this.reference = { state: this.current, time: 0 };  }  /** The method, the step, the spring, the drag, and the gravity in force. */  get options(): IntegratorOptions {    return this.optionsValue;  }  /** The particle as it stands, after the steps so far. */  get state(): ParticleState {    return this.current;  }  /** Seconds simulated so far: the steps times their dt. */  get time(): number {    return this.clock;  }  /** How many steps have been taken. */  get ticks(): number {    return this.tick;  }  /**   * Changes the method or the forces without moving anything: the next   * step obeys the new options from the same state, and the exact   * solution is measured afresh from here, so the error readouts stay   * meaningful. Invalid options throw and leave the old ones in force.   */  retune(options: IntegratorOptions): void {    ParticleIntegrator.validate(options);    this.optionsValue = options;    this.reference = { state: this.current, time: this.clock };  }  private static validate(o: IntegratorOptions): void {    if (!(o.dt > 0)) throw new RangeError(`dt ${o.dt} must be positive`);    if (!(o.timeStep > 0)) throw new RangeError(`dt ${o.dt} must be positive`);    if (!(o.stiffness > 0)) throw new RangeError(`stiffness ${o.stiffness} must be positive`);    if ((o.drag ?? 0) < 0 || (o.gravity ?? 0) < 0) throw new RangeError("drag and gravity must not be negative");    // The closed form below is the underdamped one; heavier drag would need the other two cases.    if ((o.drag ?? 0) ** 2 >= 4 * o.stiffness) throw new RangeError(`drag ${o.drag} is too heavy for stiffness ${o.stiffness}`);    if (!["euler", "semi-implicit", "verlet"].includes(o.method)) throw new RangeError(`unknown method ${o.method}`);  }  /** The pulls at a state: the spring toward the origin, drag against the velocity, gravity down, and their sum. */  static forces(state: ParticleState, o: IntegratorOptions): ParticleForces {    const spring = state.position.scale(-o.stiffness);    const drag = state.velocity.scale(-(o.drag ?? 0));    const gravity = new Vec2(0, -(o.gravity ?? 0));    return { spring, drag, gravity, total: spring.add(drag).add(gravity) };  }  /** Kinetic energy plus the spring's and gravity's potential, per unit mass. */  static energy(state: ParticleState, o: IntegratorOptions): ParticleEnergy {    const kinetic = state.velocity.dot(state.velocity) / 2;    const potential = (o.stiffness * state.position.dot(state.position)) / 2 + (o.gravity ?? 0) * state.position.y;    return { kinetic, potential, total: kinetic + potential };  }  /**   * The exact state at time t, measured from the reference: on each axis   * a damped oscillation about the rest point (the origin across, and   * -g/k down, where gravity and the spring balance).   */  exact(t: number): ParticleState {    const o = this.options;    const k = o.stiffness;    const c = o.drag ?? 0;    const omega = Math.sqrt(k - (c * c) / 4);    const tau = t - this.reference.time;    if (tau === 0) return this.reference.state;    const decay = Math.exp((-c * tau) / 2);    const cos = Math.cos(omega * tau);    const sin = Math.sin(omega * tau);    const axis = (u0: number, du0: number, rest: number): [number, number] => {      const w0 = u0 - rest;      const b = (du0 + (c * w0) / 2) / omega;      const w = decay * (w0 * cos + b * sin);      const dw = decay * ((-c / 2) * (w0 * cos + b * sin) + omega * (b * cos - w0 * sin));      return [rest + w, dw];    };    const { position, velocity } = this.reference.state;    const [x, vx] = axis(position.x, velocity.x, 0);    const [y, vy] = axis(position.y, velocity.y, -(o.gravity ?? 0) / k);    return { position: new Vec2(x, y), velocity: new Vec2(vx, vy) };  }  /** One step of the chosen method, then an audit of the energy against the exact motion. */  step(): ParticleState {    this.tick += 1;    const dt = this.options.dt;    const timeStep = this.options.timeStep;    const p = this.current.position;    const particlePosition = this.current.position;    const v = this.current.velocity;    const particleVelocity = this.current.velocity;    const a = this.accelerationAt(p, v);    const acceleration = this.accelerationAt(particlePosition, particleVelocity);    let next: ParticleState;    if (this.options.method === "euler") {      // Explicit: move along the velocity the step began with, then update it.      const moved = this.move(p, v, dt);      const moved = this.move(particlePosition, particleVelocity, timeStep);      const pushed = this.push(v, a, dt);      const pushed = this.push(particleVelocity, acceleration, timeStep);      next = { position: moved, velocity: pushed };    } else if (this.options.method === "semi-implicit") {      // Symplectic: update the velocity first, then move along the new one.      const pushed = this.push(v, a, dt);      const pushed = this.push(particleVelocity, acceleration, timeStep);      const moved = this.move(p, pushed, dt);      const moved = this.move(particlePosition, pushed, timeStep);      next = { position: moved, velocity: pushed };    } else {      // Velocity Verlet: half the kick, the move, the force at the new place, the mean of both kicks.      const half = v.add(a.scale(dt / 2));      const half = particleVelocity.add(acceleration.scale(timeStep / 2));      const moved = this.move(p, half, dt);      const moved = this.move(particlePosition, half, timeStep);      const a2 = this.accelerationAt(moved, half);      const pushed = this.push(v, a.add(a2).scale(0.5), dt);      const pushed = this.push(particleVelocity, acceleration.add(a2).scale(0.5), timeStep);      next = { position: moved, velocity: pushed };    }    this.current = next;    this.clock += dt;    this.clock += timeStep;    this.audit();    return next;  }  /** Steps `ticks` times and returns the history, the starting state first. */  run(ticks: number): ParticleState[] {    if (!Number.isInteger(ticks) || ticks < 0) throw new RangeError(`ticks ${ticks} must be a whole number`);    const history: ParticleState[] = [this.current];    for (let t = 0; t < ticks; t++) history.push(this.step());    return history;  }  /** The acceleration at a place and velocity: the forces summed, per unit mass. */  private accelerationAt(p: Vec2, v: Vec2): Vec2 {  private accelerationAt(particlePosition: Vec2, particleVelocity: Vec2): Vec2 {    const f = ParticleIntegrator.forces({ position: p, velocity: v }, this.options);    const f = ParticleIntegrator.forces({ position: particlePosition, velocity: particleVelocity }, this.options);    // ▸ sum the spring, the drag, and gravity at this place: the acceleration    this.tracer?.emit("force", { tick: this.tick, x: p.x, y: p.y, sx: f.spring.x, sy: f.spring.y, dx: f.drag.x, dy: f.drag.y, gx: f.gravity.x, gy: f.gravity.y, ax: f.total.x, ay: f.total.y });    this.tracer?.emit("force", { tick: this.tick, x: particlePosition.x, y: particlePosition.y, sx: f.spring.x, sy: f.spring.y, dx: f.drag.x, dy: f.drag.y, gx: f.gravity.x, gy: f.gravity.y, ax: f.total.x, ay: f.total.y });    return f.total;  }  /** The velocity after a kick of a for dt seconds. */  private push(v: Vec2, a: Vec2, dt: number): Vec2 {  private push(particleVelocity: Vec2, acceleration: Vec2, timeStep: number): Vec2 {    const pushed = v.add(a.scale(dt));    const pushed = particleVelocity.add(acceleration.scale(timeStep));    // ▸ kick the velocity by the acceleration times dt    this.tracer?.emit("accelerate", { tick: this.tick, ax: a.x, ay: a.y, vx: pushed.x, vy: pushed.y });    this.tracer?.emit("accelerate", { tick: this.tick, ax: acceleration.x, ay: acceleration.y, vx: pushed.x, vy: pushed.y });    return pushed;  }  /** The position after moving along v for dt seconds. */  private move(p: Vec2, v: Vec2, dt: number): Vec2 {  private move(particlePosition: Vec2, particleVelocity: Vec2, timeStep: number): Vec2 {    const moved = p.add(v.scale(dt));    const moved = particlePosition.add(particleVelocity.scale(timeStep));    // ▸ move along a velocity for dt: the old one, or the one just kicked    this.tracer?.emit("move", { tick: this.tick, vx: v.x, vy: v.y, x: moved.x, y: moved.y });    this.tracer?.emit("move", { tick: this.tick, vx: particleVelocity.x, vy: particleVelocity.y, x: moved.x, y: moved.y });    return moved;  }  /** The energy now against the exact motion's, and how far the particle is from where it should be. */  private audit(): void {    const energy = ParticleIntegrator.energy(this.current, this.options);    const truth = this.exact(this.clock);    const exact = ParticleIntegrator.energy(truth, this.options).total;    const error = this.current.position.sub(truth.position).length();    // ▸ audit: kinetic plus potential, against the exact motion    this.tracer?.emit("energy", { tick: this.tick, t: this.clock, kinetic: energy.kinetic, potential: energy.potential, total: energy.total, exact, drift: energy.total - exact, error });  }}

steps =: 50
time =: 5.00
drift =: +158.77
error =: 4.85

50 steps, 5.00 s simulated: the energy drifted +158.77 from the exact motion's, and the particle ended 4.85 units from its exact place.

import { Vec2 } from "../linear/vec";import type { Tracer } from "../trace";/** How a step turns the acceleration into a new velocity and position. */export type IntegrationMethod = "euler" | "semi-implicit" | "verlet";/** Where the particle is and how it is moving, in world units and units per second. */export interface ParticleState { position: Vec2; velocity: Vec2;}export interface IntegratorOptions { method: IntegrationMethod; /** How many seconds one step advances the clock. */ dt: number; timeStep: number; /** The spring to the origin: a pull of -k p. */ stiffness: number; /** Linear drag: a pull of -c v against the motion (default 0, none). */ drag?: number; /** The downward pull, in units per second squared (default 0, none). */ gravity?: number;}/** The pulls on the particle at a state, per unit mass, so each is an acceleration. */export interface ParticleForces { spring: Vec2; drag: Vec2; gravity: Vec2; total: Vec2;}export interface ParticleEnergy { kinetic: number; potential: number; total: number;}/** * One particle of unit mass on a spring to the origin, under optional * drag and gravity, stepped by one of three integrators. Explicit Euler * moves along the old velocity and then updates it, so every step adds * energy and an orbit spirals out. Semi-implicit (symplectic) Euler * updates the velocity first and moves along the new one: the same * arithmetic in the other order, and the energy stays bounded. * Velocity Verlet takes half the kick, moves, measures the force at the * new place, and finishes the kick with the mean of the two: second * order, and symmetric in time. The forces are linear, so the exact * motion is a damped oscillation about the rest point on each axis, and * `exact(t)` gives it in closed form for the error and drift readouts; * `retune` measures it from the state the retune found, so a live run * keeps a truth to compare against. Every step traces the force, the * velocity update, the move, and an energy audit, in the order the * method does them. * * Time: O(1) per step. Space: O(1). */export class ParticleIntegrator { private current: ParticleState; private clock = 0; private tick = 0; private optionsValue: IntegratorOptions; /** The state and time the exact solution is measured from: the start, or the last retune. */ private reference: { state: ParticleState; time: number }; constructor( initial: ParticleState, options: IntegratorOptions, private readonly tracer?: Tracer, ) { ParticleIntegrator.validate(options); this.optionsValue = options; this.current = { position: initial.position, velocity: initial.velocity }; this.reference = { state: this.current, time: 0 }; } /** The method, the step, the spring, the drag, and the gravity in force. */ get options(): IntegratorOptions { return this.optionsValue; } /** The particle as it stands, after the steps so far. */ get state(): ParticleState { return this.current; } /** Seconds simulated so far: the steps times their dt. */ get time(): number { return this.clock; } /** How many steps have been taken. */ get ticks(): number { return this.tick; } /** * Changes the method or the forces without moving anything: the next * step obeys the new options from the same state, and the exact * solution is measured afresh from here, so the error readouts stay * meaningful. Invalid options throw and leave the old ones in force. */ retune(options: IntegratorOptions): void { ParticleIntegrator.validate(options); this.optionsValue = options; this.reference = { state: this.current, time: this.clock }; } private static validate(o: IntegratorOptions): void { if (!(o.dt > 0)) throw new RangeError(`dt ${o.dt} must be positive`); if (!(o.timeStep > 0)) throw new RangeError(`dt ${o.dt} must be positive`); if (!(o.stiffness > 0)) throw new RangeError(`stiffness ${o.stiffness} must be positive`); if ((o.drag ?? 0) < 0 || (o.gravity ?? 0) < 0) throw new RangeError("drag and gravity must not be negative"); // The closed form below is the underdamped one; heavier drag would need the other two cases. if ((o.drag ?? 0) ** 2 >= 4 * o.stiffness) throw new RangeError(`drag ${o.drag} is too heavy for stiffness ${o.stiffness}`); if (!["euler", "semi-implicit", "verlet"].includes(o.method)) throw new RangeError(`unknown method ${o.method}`); } /** The pulls at a state: the spring toward the origin, drag against the velocity, gravity down, and their sum. */ static forces(state: ParticleState, o: IntegratorOptions): ParticleForces { const spring = state.position.scale(-o.stiffness); const drag = state.velocity.scale(-(o.drag ?? 0)); const gravity = new Vec2(0, -(o.gravity ?? 0)); return { spring, drag, gravity, total: spring.add(drag).add(gravity) }; } /** Kinetic energy plus the spring's and gravity's potential, per unit mass. */ static energy(state: ParticleState, o: IntegratorOptions): ParticleEnergy { const kinetic = state.velocity.dot(state.velocity) / 2; const potential = (o.stiffness * state.position.dot(state.position)) / 2 + (o.gravity ?? 0) * state.position.y; return { kinetic, potential, total: kinetic + potential }; } /** * The exact state at time t, measured from the reference: on each axis * a damped oscillation about the rest point (the origin across, and * -g/k down, where gravity and the spring balance). */ exact(t: number): ParticleState { const o = this.options; const k = o.stiffness; const c = o.drag ?? 0; const omega = Math.sqrt(k - (c * c) / 4); const tau = t - this.reference.time; if (tau === 0) return this.reference.state; const decay = Math.exp((-c * tau) / 2); const cos = Math.cos(omega * tau); const sin = Math.sin(omega * tau); const axis = (u0: number, du0: number, rest: number): [number, number] => { const w0 = u0 - rest; const b = (du0 + (c * w0) / 2) / omega; const w = decay * (w0 * cos + b * sin); const dw = decay * ((-c / 2) * (w0 * cos + b * sin) + omega * (b * cos - w0 * sin)); return [rest + w, dw]; }; const { position, velocity } = this.reference.state; const [x, vx] = axis(position.x, velocity.x, 0); const [y, vy] = axis(position.y, velocity.y, -(o.gravity ?? 0) / k); return { position: new Vec2(x, y), velocity: new Vec2(vx, vy) }; } /** One step of the chosen method, then an audit of the energy against the exact motion. */ step(): ParticleState { this.tick += 1; const dt = this.options.dt; const timeStep = this.options.timeStep; const p = this.current.position; const particlePosition = this.current.position; const v = this.current.velocity; const particleVelocity = this.current.velocity; const a = this.accelerationAt(p, v); const acceleration = this.accelerationAt(particlePosition, particleVelocity); let next: ParticleState; if (this.options.method === "euler") { // Explicit: move along the velocity the step began with, then update it. const moved = this.move(p, v, dt); const moved = this.move(particlePosition, particleVelocity, timeStep); const pushed = this.push(v, a, dt); const pushed = this.push(particleVelocity, acceleration, timeStep); next = { position: moved, velocity: pushed }; } else if (this.options.method === "semi-implicit") { // Symplectic: update the velocity first, then move along the new one. const pushed = this.push(v, a, dt); const pushed = this.push(particleVelocity, acceleration, timeStep); const moved = this.move(p, pushed, dt); const moved = this.move(particlePosition, pushed, timeStep); next = { position: moved, velocity: pushed }; } else { // Velocity Verlet: half the kick, the move, the force at the new place, the mean of both kicks. const half = v.add(a.scale(dt / 2)); const half = particleVelocity.add(acceleration.scale(timeStep / 2)); const moved = this.move(p, half, dt); const moved = this.move(particlePosition, half, timeStep); const a2 = this.accelerationAt(moved, half); const pushed = this.push(v, a.add(a2).scale(0.5), dt); const pushed = this.push(particleVelocity, acceleration.add(a2).scale(0.5), timeStep); next = { position: moved, velocity: pushed }; } this.current = next; this.clock += dt; this.clock += timeStep; this.audit(); return next; } /** Steps `ticks` times and returns the history, the starting state first. */ run(ticks: number): ParticleState[] { if (!Number.isInteger(ticks) || ticks < 0) throw new RangeError(`ticks ${ticks} must be a whole number`); const history: ParticleState[] = [this.current]; for (let t = 0; t < ticks; t++) history.push(this.step()); return history; } /** The acceleration at a place and velocity: the forces summed, per unit mass. */ private accelerationAt(p: Vec2, v: Vec2): Vec2 { private accelerationAt(particlePosition: Vec2, particleVelocity: Vec2): Vec2 { const f = ParticleIntegrator.forces({ position: p, velocity: v }, this.options); const f = ParticleIntegrator.forces({ position: particlePosition, velocity: particleVelocity }, this.options); // ▸ sum the spring, the drag, and gravity at this place: the acceleration this.tracer?.emit("force", { tick: this.tick, x: p.x, y: p.y, sx: f.spring.x, sy: f.spring.y, dx: f.drag.x, dy: f.drag.y, gx: f.gravity.x, gy: f.gravity.y, ax: f.total.x, ay: f.total.y }); this.tracer?.emit("force", { tick: this.tick, x: particlePosition.x, y: particlePosition.y, sx: f.spring.x, sy: f.spring.y, dx: f.drag.x, dy: f.drag.y, gx: f.gravity.x, gy: f.gravity.y, ax: f.total.x, ay: f.total.y }); return f.total; } /** The velocity after a kick of a for dt seconds. */ private push(v: Vec2, a: Vec2, dt: number): Vec2 { private push(particleVelocity: Vec2, acceleration: Vec2, timeStep: number): Vec2 { const pushed = v.add(a.scale(dt)); const pushed = particleVelocity.add(acceleration.scale(timeStep)); // ▸ kick the velocity by the acceleration times dt this.tracer?.emit("accelerate", { tick: this.tick, ax: a.x, ay: a.y, vx: pushed.x, vy: pushed.y }); this.tracer?.emit("accelerate", { tick: this.tick, ax: acceleration.x, ay: acceleration.y, vx: pushed.x, vy: pushed.y }); return pushed; } /** The position after moving along v for dt seconds. */ private move(p: Vec2, v: Vec2, dt: number): Vec2 { private move(particlePosition: Vec2, particleVelocity: Vec2, timeStep: number): Vec2 { const moved = p.add(v.scale(dt)); const moved = particlePosition.add(particleVelocity.scale(timeStep)); // ▸ move along a velocity for dt: the old one, or the one just kicked this.tracer?.emit("move", { tick: this.tick, vx: v.x, vy: v.y, x: moved.x, y: moved.y }); this.tracer?.emit("move", { tick: this.tick, vx: particleVelocity.x, vy: particleVelocity.y, x: moved.x, y: moved.y }); return moved; } /** The energy now against the exact motion's, and how far the particle is from where it should be. */ private audit(): void { const energy = ParticleIntegrator.energy(this.current, this.options); const truth = this.exact(this.clock); const exact = ParticleIntegrator.energy(truth, this.options).total; const error = this.current.position.sub(truth.position).length(); // ▸ audit: kinetic plus potential, against the exact motion this.tracer?.emit("energy", { tick: this.tick, t: this.clock, kinetic: energy.kinetic, potential: energy.potential, total: energy.total, exact, drift: energy.total - exact, error }); }}