Integración y el paso de tiempo

Integrar es convertir una aceleración en un movimiento, y hay tres formas de hacerlo que se diferencian en el orden de dos líneas de código. Esa diferencia trivial decide si un sistema de partículas orbitando un atractor se mantiene estable durante horas o se desintegra en veinte segundos, y la explicación no está en la precisión sino en si el método conserva o no una cantidad relacionada con la energía. Es de los pocos sitios de la programación gráfica donde la matemática correcta y el código corto coinciden exactamente.

Los tres métodos

Todos parten de lo mismo: una aceleración a calculada a partir del estado, y un paso h. La diferencia está en qué velocidad se usa para mover la posición.

Euler explícito mueve la posición con la velocidad antigua y después actualiza la velocidad:

let a = aceleracion(pos, vel);
pos = pos + vel * h;          // velocidad vieja
vel = vel + a * h;

Euler semi-implícito, también llamado simpléctico o de Cromer, actualiza primero la velocidad y mueve la posición con la nueva:

let a = aceleracion(pos, vel);
vel = vel + a * h;
pos = pos + vel * h;          // velocidad nueva

Verlet en velocidad usa la media de las aceleraciones antes y después del paso:

let a0 = aceleracion(pos, vel);
pos = pos + vel * h + a0 * (0.5 * h * h);
let a1 = aceleracion(pos, vel);
vel = vel + (a0 + a1) * (0.5 * h);

Los dos primeros son de primer orden: el error local por paso es proporcional a h². Verlet es de segundo: el error local es proporcional a h³. Uno esperaría que la diferencia relevante fuera esa, y no lo es.

Por qué el explícito explota

Toma un oscilador armónico —un muelle, una órbita circular, cualquier fuerza restauradora— y sigue la energía total del sistema paso a paso. Con Euler explícito, la energía crece en cada paso, en una cantidad proporcional a h². No es un error de redondeo que se pueda ignorar: es un sesgo sistemático en una dirección. Un planeta en órbita se aleja en espiral, un muelle amplifica su oscilación, un sistema de partículas atrapadas en un pozo de potencial acaba escapando.

Con Euler semi-implícito la energía no se conserva exactamente, pero oscila alrededor del valor correcto sin deriva. La órbita no es exacta —precesa levemente— pero es cerrada y estable indefinidamente. Esa propiedad es lo que significa que un integrador sea simpléctico: conserva una cantidad muy próxima a la energía real, en vez de la energía real, y esa cantidad no deriva.

La consecuencia práctica es contundente: Euler semi-implícito cuesta exactamente lo mismo que el explícito y es cualitativamente mejor. No hay ninguna razón para escribir el explícito jamás. Y sin embargo se escribe todo el tiempo, porque el orden natural en que uno piensa el problema es “muevo la partícula y luego actualizo la velocidad”.

Verlet en velocidad también es simpléctico y además de segundo orden, pero cuesta dos evaluaciones de la aceleración por paso. Merece la pena cuando la aceleración es barata y la precisión importa —trayectorias balísticas, órbitas de las que se espera exactitud— y no merece la pena cuando la aceleración implica una búsqueda de vecinos, porque duplicar eso duplica el coste del frame entero.

Para un sistema de partículas visual, la respuesta correcta casi siempre es Euler semi-implícito.

El kernel

struct Params {
  conteo: u32,
  h:      f32,          // paso de simulacion, en segundos
  tiempo: f32,
  _r:     f32,
};
@group(0) @binding(0) var<uniform> p: Params;

struct Particula { pos: vec3f, vida: f32, vel: vec3f, masa: f32 };
@group(0) @binding(1) var<storage, read_write> parts: array<Particula>;

const G: vec3f = vec3f(0.0, -9.81, 0.0);
const ARRASTRE: f32 = 0.12;
const VEL_MAX: f32 = 60.0;

fn aceleracion(q: Particula) -> vec3f {
  // Gravedad mas arrastre lineal. El arrastre estabiliza el integrador.
  return G - q.vel * ARRASTRE;
}

@compute @workgroup_size(64)
fn integrar(@builtin(global_invocation_id) gid: vec3u) {
  let i = gid.x;
  if (i >= p.conteo) { return; }

  var q = parts[i];
  if (q.vida <= 0.0) { return; }

  // Euler semi-implicito: primero la velocidad, luego la posicion.
  q.vel = q.vel + aceleracion(q) * p.h;

  // Saturar la velocidad evita que un paso grande mande la particula
  // al infinito y con ella el resto del sistema si hay interaccion.
  let v = length(q.vel);
  if (v > VEL_MAX) { q.vel = q.vel * (VEL_MAX / v); }

  q.pos = q.pos + q.vel * p.h;

  // Rebote contra el suelo con perdida de energia.
  if (q.pos.y < 0.0) {
    q.pos.y = -q.pos.y * 0.4;
    q.vel.y = -q.vel.y * 0.4;
    q.vel = vec3f(q.vel.x * 0.85, q.vel.y, q.vel.z * 0.85);
  }

  q.vida = q.vida - p.h;
  parts[i] = q;
}

Los return tempranos son legales porque no hay ninguna barrera en el kernel. La saturación de velocidad no es un adorno: sin ella, un paso de tiempo anómalo produce una velocidad enorme, la partícula sale del dominio, y si hay una rejilla espacial de por medio su índice de celda se desborda y corrompe la estructura de todos los demás.

El acumulador de tiempo

El intervalo entre frames varía. Si se lo pasas directamente al integrador, la simulación depende del framerate: un usuario a 144 Hz ve un comportamiento distinto que uno a 60, y una caída puntual de rendimiento produce un salto visible. La solución estándar es fijar el paso de simulación y acumular el tiempo real.

const H = 1 / 120;              // paso fijo de simulacion
const MAX_PASOS = 8;            // techo por frame
let acumulado = 0;
let anterior = performance.now();

function frame(ahora) {
  let dt = (ahora - anterior) / 1000;
  anterior = ahora;

  // Techo al delta: si la pestana estuvo en segundo plano, dt puede
  // valer treinta segundos y no queremos simularlos.
  dt = Math.min(dt, 0.25);
  acumulado += dt;

  let pasos = 0;
  const enc = device.createCommandEncoder();
  const pass = enc.beginComputePass();
  pass.setPipeline(pipeIntegrar);
  pass.setBindGroup(0, bgSim);
  while (acumulado >= H && pasos < MAX_PASOS) {
    pass.dispatchWorkgroups(Math.ceil(N / 64));
    acumulado -= H;
    pasos++;
  }
  pass.end();
  // ... render ...
  device.queue.submit([enc.finish()]);
  requestAnimationFrame(frame);
}

Fíjate en que los pasos se encolan como dispatches consecutivos en el mismo pass. Cada uno ve el resultado del anterior por la garantía entre dispatches, así que no hace falta nada más, y el coste de CPU de encolar ocho dispatches es despreciable.

El MAX_PASOS evita la espiral de la muerte: si simular un paso tarda más que el paso mismo, el acumulado crece cada frame, el bucle hace más pasos, tarda más todavía, y el programa se cuelga. Con el techo, la simulación se ralentiza respecto al tiempo real en vez de bloquear la página. Es una degradación honesta, y con partículas visuales nadie la nota.

El uniform con p.h se escribe una vez y no cambia, que es otra ventaja de tener el paso fijo: no hay que actualizar el buffer entre dispatches, cosa que dentro de un pass no se puede hacer.

Si quieres suavidad visual con un paso de simulación bajo, la técnica complementaria es interpolar en el shader de render entre el estado anterior y el actual usando acumulado / H como factor. Eso exige guardar la posición anterior, cuatro bytes más por partícula, y para partículas pequeñas y rápidas rara vez compensa.