26. Euler semimplícito y Verlet: alternativas de integración

26.1 Por qué considerar otro integrador

Euler explícito aproxima el siguiente estado usando solamente las tasas conocidas al comienzo del paso. Es simple y útil, pero puede acumular error rápidamente. En un movimiento oscilatorio puede aumentar la amplitud sin que la regla cinemática indique ese crecimiento.

Un integrador no cambia la ecuación que queremos resolver: cambia la receta numérica utilizada para avanzar entre estados discretos. Compararemos recetas que requieren cantidades de información y operaciones diferentes.

26.2 Nuestro punto de partida: Euler explícito

Posición y velocidad usan valores del estado n. Si la aceleración depende de la posición, primero evaluamos aₙ=a(xₙ,vₙ,tₙ). Su costo por paso es pequeño y normalmente requiere una evaluación de la aceleración.

La posición ignora cómo cambia la velocidad durante el intervalo. El método es de primer orden: bajo condiciones adecuadas, su error global disminuye aproximadamente en proporción a Δt.

26.3 Euler semimplícito

Euler semimplícito —también llamado Euler simpléctico en ciertos sistemas— actualiza primero la velocidad y usa esa velocidad nueva para la posición:

La única diferencia visible es el subíndice de v en la segunda ecuación, pero el comportamiento puede cambiar notablemente. Sigue siendo un método de primer orden y no es “Euler explícito escrito en otro orden”.

26.4 Implementar Euler semimplícito

function pasoEulerSemiImplicito(estado, calcularAceleracion, dt) {
  const aceleracionMps2 = calcularAceleracion(estado);
  const nuevaVelocidadMps = estado.velocidadMps
    + aceleracionMps2 * dt;
  return {
    tiempoS: estado.tiempoS + dt,
    posicionM: estado.posicionM + nuevaVelocidadMps * dt,
    velocidadMps: nuevaVelocidadMps
  };
}

const aceleracion = estado => -4 * estado.posicionM;
console.log(pasoEulerSemiImplicito(
  { tiempoS: 0, posicionM: 1, velocidadMps: 0 },
  aceleracion,
  0.1
));

Para este oscilador, la aceleración inicial es −4 m/s². La nueva velocidad vale −0,4 m/s y la posición avanza con ella hasta 0,96 m.

26.5 Comparación de un paso

Con x₀=0 m, v₀=2 m/s, a=3 m/s² y Δt=0,5 s:

En este caso uno subestima y el otro sobreestima la posición por la misma cantidad. Eso no implica que promediarlos sea siempre un integrador adecuado ni que ambos se comporten simétricamente para cualquier aceleración.

26.6 Sistemas oscilatorios y energía

Para comparar integradores usaremos una aceleración dada por a(x)=−ω²x. Su solución exacta oscila. Si elegimos masa y constante de referencia iguales a uno, una magnitud útil para diagnosticar el cálculo es:

En la solución ideal E permanece constante. Euler explícito tiende a hacerla crecer. Euler semimplícito suele mantenerla oscilando alrededor de un valor acotado si el paso es apropiado. La energía se usa aquí como prueba numérica, no para desarrollar las causas del movimiento.

26.7 La idea de Verlet de posición

La expansión de Taylor de la posición hacia adelante y hacia atrás es:

Al sumar ambas expresiones, los términos de velocidad se cancelan. Despejando la posición futura obtenemos Verlet de posición:

El método usa la posición actual y la anterior. Su error global en posición es de segundo orden bajo condiciones habituales: al reducir Δt a la mitad, el error tiende a disminuir aproximadamente por un factor de cuatro.

26.8 Inicializar la posición anterior

El estado inicial suele proporcionar x₀ y v₀, pero Verlet también necesita x₋₁. Una aproximación coherente de segundo orden se obtiene retrocediendo mediante Taylor:

function iniciarVerlet(posicionM, velocidadMps, aceleracionMps2, dt) {
  return {
    tiempoS: 0,
    posicionAnteriorM: posicionM - velocidadMps * dt
      + 0.5 * aceleracionMps2 * dt ** 2,
    posicionM
  };
}

function pasoVerlet(estado, aceleracionMps2, dt) {
  const siguienteM = 2 * estado.posicionM
    - estado.posicionAnteriorM
    + aceleracionMps2 * dt ** 2;
  return {
    tiempoS: estado.tiempoS + dt,
    posicionAnteriorM: estado.posicionM,
    posicionM: siguienteM
  };
}

const estado = iniciarVerlet(1, 0, -4, 0.1);
console.log(estado);
console.log(pasoVerlet(estado, -4, 0.1));

26.9 Recuperar una velocidad con Verlet

Verlet de posición no guarda la velocidad como variable principal. Podemos estimarla de forma centrada cuando disponemos de la posición anterior y la siguiente:

Esta velocidad corresponde al instante central n, no exactamente al extremo nuevo. Esa diferencia importa si una regla depende de la velocidad, si debemos informar estados sincronizados o si queremos calcular colisiones con rozamiento.

26.10 Velocity Verlet

Velocity Verlet conserva posición y velocidad sincronizadas. Primero calcula la nueva posición usando aₙ, evalúa la aceleración en el estado nuevo y después actualiza la velocidad con el promedio de ambas aceleraciones:

Si la aceleración depende únicamente de posición y tiempo, el esquema es directo. Si depende también de la velocidad, evaluar aₙ₊₁ puede requerir una adaptación o un método distinto.

26.11 Implementar Velocity Verlet

function pasoVelocityVerlet(estado, aceleracionEn, dt) {
  const a0 = aceleracionEn(estado.posicionM, estado.tiempoS);
  const nuevoTiempoS = estado.tiempoS + dt;
  const nuevaPosicionM = estado.posicionM
    + estado.velocidadMps * dt
    + 0.5 * a0 * dt ** 2;
  const a1 = aceleracionEn(nuevaPosicionM, nuevoTiempoS);
  const nuevaVelocidadMps = estado.velocidadMps
    + 0.5 * (a0 + a1) * dt;

  return {
    tiempoS: nuevoTiempoS,
    posicionM: nuevaPosicionM,
    velocidadMps: nuevaVelocidadMps,
    aceleracionMps2: a1
  };
}

const aceleracion = posicionM => -4 * posicionM;
console.log(pasoVelocityVerlet(
  { tiempoS: 0, posicionM: 1, velocidadMps: 0 },
  aceleracion,
  0.1
));

Este paso evalúa la aceleración dos veces, aunque en un bucle podemos reutilizar aₙ₊₁ como aₙ del paso siguiente y reducir el costo efectivo.

26.12 Aplicación vectorial

Las ecuaciones se aplican componente a componente. En dos dimensiones, Velocity Verlet utiliza vectores para posición, velocidad y aceleración:

const sumar = (a, b) => ({ x: a.x + b.x, y: a.y + b.y });
const escalar = (v, k) => ({ x: v.x * k, y: v.y * k });

function pasoVelocityVerlet2D(estado, aceleracionEn, dt) {
  const a0 = aceleracionEn(estado.posicionM, estado.tiempoS);
  const tiempoS = estado.tiempoS + dt;
  const posicionM = sumar(
    sumar(estado.posicionM, escalar(estado.velocidadMps, dt)),
    escalar(a0, 0.5 * dt ** 2)
  );
  const a1 = aceleracionEn(posicionM, tiempoS);
  const velocidadMps = sumar(
    estado.velocidadMps,
    escalar(sumar(a0, a1), 0.5 * dt)
  );
  return { tiempoS, posicionM, velocidadMps, aceleracionMps2: a1 };
}

const gravedad = () => ({ x: 0, y: -9.8 });
console.log(pasoVelocityVerlet2D({
  tiempoS: 0,
  posicionM: { x: 0, y: 2 },
  velocidadMps: { x: 6, y: 8 }
}, gravedad, 0.1));

26.13 Aceleración constante como caso especial

Cuando la aceleración es constante, aₙ₊₁=aₙ=a. Velocity Verlet se reduce a:

Son exactamente las ecuaciones del movimiento uniformemente acelerado. Por eso este método reproduce ese caso salvo errores de coma flotante. No debemos concluir que será exacto cuando la aceleración cambie con posición, velocidad o tiempo.

26.14 Actividad interactiva: comparar integradores

El oscilador comienza en x=1 m y v=0, con aceleración a=−ω²x. Elegí un método y un paso. La curva azul es la solución exacta; la naranja es la secuencia numérica.

26.15 Comparación de propiedades

“Orden 2” describe cómo disminuye el error al reducir Δt bajo supuestos adecuados. No garantiza por sí solo menor error para cualquier paso, problema o implementación.

26.16 Costo por paso y reutilización

Si calcular la aceleración es costoso, debemos contar evaluaciones:

• Euler explícito y semimplícito necesitan normalmente una evaluación por paso.
• Verlet de posición necesita una evaluación en xₙ.
• Velocity Verlet parece necesitar aₙ y aₙ₊₁, pero puede conservar aₙ del estado y evaluar solo aₙ₊₁.

function pasoVelocityVerletConAceleracion(estado, aceleracionEn, dt) {
  const posicionM = estado.posicionM + estado.velocidadMps * dt
    + 0.5 * estado.aceleracionMps2 * dt ** 2;
  const tiempoS = estado.tiempoS + dt;
  const aceleracionMps2 = aceleracionEn(posicionM, tiempoS);
  const velocidadMps = estado.velocidadMps
    + 0.5 * (estado.aceleracionMps2 + aceleracionMps2) * dt;
  return { tiempoS, posicionM, velocidadMps, aceleracionMps2 };
}

const aceleracionEn = x => -4 * x;
console.log(pasoVelocityVerletConAceleracion({
  tiempoS: 0,
  posicionM: 1,
  velocidadMps: 0,
  aceleracionMps2: -4
}, aceleracionEn, 0.1));

26.17 Colisiones y cambios discontinuos

Los integradores desarrollados suponen una evolución suave durante cada paso. Una colisión, un rebote o una restricción puede cambiar la velocidad de manera instantánea. En ese evento debemos aplicar explícitamente la regla correspondiente.

Después de modificar posición o velocidad, también hay que restaurar la coherencia interna del integrador. Velocity Verlet necesita recalcular la aceleración; Verlet de posición puede requerir reconstruir la posición anterior para que la nueva velocidad implícita sea la deseada.

26.18 Elegir un método

La elección debe verificarse con el sistema real, distintos pasos y magnitudes relevantes. No hay un integrador universalmente mejor.

26.19 Pruebas comparativas

Una misma interfaz permite ejecutar varios métodos sobre el mismo caso:

function energiaOscilador(estado, omega) {
  return 0.5 * estado.velocidadMps ** 2
    + 0.5 * omega ** 2 * estado.posicionM ** 2;
}

function pasoSemiImplicito(estado, omega, dt) {
  const a = -(omega ** 2) * estado.posicionM;
  const v = estado.velocidadMps + a * dt;
  return {
    posicionM: estado.posicionM + v * dt,
    velocidadMps: v
  };
}

let estado = { posicionM: 1, velocidadMps: 0 };
const energiaInicial = energiaOscilador(estado, 2);
for (let n = 0; n < 100; n++) {
  estado = pasoSemiImplicito(estado, 2, 0.05);
}

console.log(energiaInicial);
console.log(energiaOscilador(estado, 2));

Además de posición final, conviene registrar error máximo, energía mínima y máxima, amplitud y costo. Un único instante puede ocultar un mal comportamiento ocurrido antes.

26.20 Errores frecuentes

• Llamar Euler explícito a cualquier orden de asignación: el subíndice de la velocidad utilizada define el esquema.
• Usar Velocity Verlet con a dependiente de v sin adaptación: aₙ₊₁ no puede evaluarse directamente si aún falta vₙ₊₁.
• Inicializar Verlet con x₋₁=x₀: eso impone una velocidad implícita que puede no coincidir con v₀.
• Calcular la velocidad centrada en el instante equivocado: su marca temporal es la del punto central.
• Confundir orden con estabilidad absoluta: un método de mayor orden también puede fallar con un paso inadecuado.
• Evaluar dos veces sin necesidad: conservar la aceleración nueva para el paso siguiente.
• No reparar el estado tras una colisión: las variables internas pueden quedar incompatibles.
• Comparar métodos con distintos Δt o condiciones: mantener el caso de prueba controlado.

26.21 Ejercicios propuestos

1. Aplicá un paso de Euler semimplícito con x=2 m, v=1 m/s, a=−3 m/s² y Δt=0,2 s.
2. Compará ese resultado con Euler explícito y con la solución de aceleración constante.
3. Calculá x₋₁ para iniciar Verlet con x₀=5 m, v₀=2 m/s, a₀=−1 m/s² y Δt=0,1 s.
4. Usá Verlet para hallar x₁ si x₋₁=4,795 m, x₀=5 m, a₀=−1 m/s² y Δt=0,1 s.
5. Aplicá un paso de Velocity Verlet al mismo estado y aceleración constante.
6. Explicá por qué Velocity Verlet requiere conocer aₙ₊₁.
7. ¿Qué método elegirías si necesitás informar velocidad en cada estado: Verlet de posición o Velocity Verlet?
8. Diseñá una comparación de convergencia entre Euler semimplícito y Velocity Verlet.

26.22 Ideas para recordar

• El orden de actualización forma parte de la definición del integrador.
• Euler semimplícito actualiza primero v y luego x con la velocidad nueva.
• Verlet de posición usa xₙ y xₙ₋₁ y no conserva directamente v.
• Velocity Verlet mantiene x y v mediante aceleraciones de ambos extremos.
• Con aceleración constante, Velocity Verlet reproduce las ecuaciones analíticas.
• Un orden mayor suele reducir más rápido el error, pero no garantiza estabilidad universal.
• Energía, amplitud y error máximo revelan problemas que la posición final puede ocultar.
• El método, el paso temporal y el modelo deben verificarse juntos.

En el próximo tema estudiaremos error acumulado, estabilidad y elección del paso temporal.