33. Simular un sistema masa-resorte con integración numérica

33.1 Qué vamos a simular

Consideramos una masa que se mueve en una dimensión, unida a un resorte lineal y a un amortiguador viscoso. La coordenada x se mide desde el equilibrio.

El caso libre usa Fext = 0. Mantener la entrada como función opcional permitirá reutilizar el mismo núcleo para una excitación conocida.

33.2 Del continuo a pasos discretos

El modelo físico define x(t) y v(t) en todo instante. El programa sólo conserva aproximaciones en una secuencia:

Integrar numéricamente significa estimar el estado siguiente a partir del actual y de las derivadas definidas por la ecuación.

33.3 Configuración y estado

Aplicamos la separación del tema anterior:

El paso temporal no es una propiedad física del resorte: pertenece al método con el que aproximamos su evolución.

33.4 Aceleración como función pura

function aceleracion(parametros, estado) {
  const { masaKg, constanteNpm, amortiguamientoNsPm } = parametros;
  const fuerzaExternaN = parametros.fuerzaExternaN?.(estado.tiempoS) ?? 0;
  return (fuerzaExternaN
    - constanteNpm * estado.posicionM
    - amortiguamientoNsPm * estado.velocidadMps) / masaKg;
}

const parametros = {
  masaKg: 1, constanteNpm: 25, amortiguamientoNsPm: 1,
  fuerzaExternaN: tiempoS => 2 * Math.cos(3 * tiempoS)
};
const estado = { tiempoS: 0, posicionM: 0.1, velocidadMps: 0 };
console.log(aceleracion(parametros, estado), "m/s²");

Para este estado, la fuerza externa vale +2 N, la elástica −2,5 N y la aceleración −0,5 m/s².

33.5 Euler explícito

El método más directo usa las derivadas del comienzo del paso:

Es sencillo y útil para entender la discretización, pero en un oscilador ideal tiende a agregar energía. Reducir Δt demora el crecimiento, aunque no cambia esa tendencia cualitativa.

33.6 Implementar Euler explícito

function pasoEulerExplicito(parametros, estado, dtS) {
  const fuerzaExternaN = parametros.fuerzaExternaN?.(estado.tiempoS) ?? 0;
  const aMps2 = (fuerzaExternaN
    - parametros.constanteNpm * estado.posicionM
    - parametros.amortiguamientoNsPm * estado.velocidadMps) / parametros.masaKg;
  return {
    tiempoS: estado.tiempoS + dtS,
    posicionM: estado.posicionM + estado.velocidadMps * dtS,
    velocidadMps: estado.velocidadMps + aMps2 * dtS
  };
}

const p = { masaKg: 1, constanteNpm: 25, amortiguamientoNsPm: 0 };
const e0 = { tiempoS: 0, posicionM: 0.1, velocidadMps: 0 };
console.log(pasoEulerExplicito(p, e0, 0.01));

33.7 Euler semimplícito

Primero actualizamos la velocidad y luego usamos esa velocidad nueva para la posición:

El cambio parece mínimo, pero su comportamiento en osciladores es muy diferente. En el caso conservativo mantiene la energía acotada dentro de un error oscilante para pasos suficientemente pequeños.

33.8 Implementar Euler semimplícito

function pasoEulerSemimplicito(parametros, estado, dtS) {
  const fuerzaExternaN = parametros.fuerzaExternaN?.(estado.tiempoS) ?? 0;
  const aMps2 = (fuerzaExternaN
    - parametros.constanteNpm * estado.posicionM
    - parametros.amortiguamientoNsPm * estado.velocidadMps) / parametros.masaKg;
  const velocidadMps = estado.velocidadMps + aMps2 * dtS;
  return {
    tiempoS: estado.tiempoS + dtS,
    posicionM: estado.posicionM + velocidadMps * dtS,
    velocidadMps
  };
}

const p = { masaKg: 1, constanteNpm: 25, amortiguamientoNsPm: 0 };
const e0 = { tiempoS: 0, posicionM: 0.1, velocidadMps: 0 };
console.log(pasoEulerSemimplicito(p, e0, 0.01));

33.9 Runge-Kutta de cuarto orden

RK4 evalúa la derivada cuatro veces dentro del paso y combina esas pendientes. Para el vector de estado y = [x, v]:

Su error global decrece mucho más rápidamente al reducir Δt que el de los métodos de Euler. A cambio requiere más evaluaciones y no conserva exactamente las propiedades geométricas de un sistema conservativo.

33.10 Implementar un paso RK4

function pasoRK4(parametros, estado, dtS) {
  const derivada = e => ({
    dx: e.velocidadMps,
    dv: ((parametros.fuerzaExternaN?.(e.tiempoS) ?? 0)
      - parametros.constanteNpm * e.posicionM
      - parametros.amortiguamientoNsPm * e.velocidadMps) / parametros.masaKg
  });
  const mover = (e, k, factor) => ({
    tiempoS: e.tiempoS + factor * dtS,
    posicionM: e.posicionM + factor * dtS * k.dx,
    velocidadMps: e.velocidadMps + factor * dtS * k.dv
  });

  const k1 = derivada(estado);
  const k2 = derivada(mover(estado, k1, 0.5));
  const k3 = derivada(mover(estado, k2, 0.5));
  const k4 = derivada(mover(estado, k3, 1));
  return {
    tiempoS: estado.tiempoS + dtS,
    posicionM: estado.posicionM + dtS * (k1.dx + 2*k2.dx + 2*k3.dx + k4.dx) / 6,
    velocidadMps: estado.velocidadMps + dtS * (k1.dv + 2*k2.dv + 2*k3.dv + k4.dv) / 6
  };
}

const p = { masaKg: 1, constanteNpm: 25, amortiguamientoNsPm: 0 };
const e0 = { tiempoS: 0, posicionM: 0.1, velocidadMps: 0 };
console.log(pasoRK4(p, e0, 0.01));

33.11 Elegir el integrador

• Euler explícito: educativo, barato y generalmente inadecuado para evoluciones oscilatorias largas.
• Euler semimplícito: simple, robusto para muchas simulaciones interactivas y mejor comportamiento energético.
• RK4: alta precisión por paso para trayectorias suaves, con mayor costo.

No existe un integrador mejor para todo propósito. Importan la precisión requerida, la duración, el costo y las propiedades que deseamos preservar.

33.12 Amortiguamiento y métodos de Verlet

Verlet clásico resulta especialmente natural cuando la aceleración depende de la posición, pero el rozamiento viscoso introduce dependencia de la velocidad.

Es posible adaptar el método, pero no conviene insertar sin más −bv en una fórmula derivada para aceleraciones independientes de v. En este tema usamos semimplícito y RK4, que aceptan directamente a(t, x, v).

33.13 Elegir Δt según la escala física

El período natural proporciona una primera escala:

En vez de pensar sólo “Δt = 0,01”, observamos cuántos pasos hay por período:

Un valor como 100 pasos por período puede ser un punto de partida, no una garantía. La adecuación se comprueba reduciendo el paso y comparando resultados.

33.14 Llegar exactamente al tiempo final

Si la duración no es múltiplo de Δt, el último paso debe acortarse:

Así evitamos terminar después del instante solicitado. También usamos una tolerancia para que los redondeos no produzcan un paso residual diminuto.

33.15 Integrar y muestrear no son lo mismo

El integrador puede avanzar cada 0,001 s y la gráfica guardar un punto cada 0,02 s. Separar ambas frecuencias reduce memoria sin degradar el cálculo interno.

Si sólo guardamos cuando el tiempo supera la próxima marca, conviene registrar también el estado inicial y el final.

33.16 Guardar instantáneas independientes

Cada muestra debe ser una copia. Si agregamos repetidamente el mismo objeto mutable, el arreglo contendrá muchas referencias al último estado.

Para millones de muestras pueden utilizarse arreglos tipados separados para tiempo, posición y velocidad.

33.17 Un simulador reutilizable

function aceleracionLibre(parametros, estado) {
  return (-parametros.constanteNpm * estado.posicionM
    - parametros.amortiguamientoNsPm * estado.velocidadMps) / parametros.masaKg;
}

function pasoEulerSemimplicito(parametros, estado, dtS) {
  const velocidadMps = estado.velocidadMps
    + aceleracionLibre(parametros, estado) * dtS;
  return {
    tiempoS: estado.tiempoS + dtS,
    posicionM: estado.posicionM + velocidadMps * dtS,
    velocidadMps
  };
}

function simular({ parametros, estadoInicial, dtS, duracionS, integrador }) {
  if (dtS <= 0 || duracionS < 0) throw new RangeError("Tiempos inválidos");
  let estado = { ...estadoInicial };
  const muestras = [{ ...estado }];

  while (estado.tiempoS < duracionS - 1e-12) {
    const paso = Math.min(dtS, duracionS - estado.tiempoS);
    estado = integrador(parametros, estado, paso);
    if (![estado.tiempoS, estado.posicionM, estado.velocidadMps].every(Number.isFinite)) {
      throw new Error("La integración produjo un estado no finito");
    }
    muestras.push({ ...estado });
  }
  return muestras;
}

const p = { masaKg: 1, constanteNpm: 25, amortiguamientoNsPm: 1 };
const e0 = { tiempoS: 0, posicionM: 0.1, velocidadMps: 0 };
const datos = simular({ parametros: p, estadoInicial: e0,
  dtS: 0.002, duracionS: 5, integrador: pasoEulerSemimplicito });
console.log({ cantidad: datos.length, estadoFinal: datos.at(-1) });

33.18 Tiempo real con acumulador

La frecuencia de dibujo del navegador es variable. No debemos usar directamente la duración de cada cuadro como gran paso físico.

Este patrón conserva un paso físico fijo aunque los cuadros visuales lleguen con intervalos irregulares.

33.19 Interpolación para dibujar

Después de consumir pasos completos suele quedar una fracción en el acumulador. Podemos interpolar entre el estado anterior y el actual para representar una posición visual suave.

La interpolación sólo afecta la vista. No debe reemplazar el estado físico ni retroalimentarse en el integrador.

33.20 Detener una simulación amortiguada

Una exponencial no llega exactamente a cero. Podemos considerar reposo cuando durante cierto intervalo se cumplen simultáneamente:

Las tolerancias deben expresarse en unidades físicas y elegirse según el propósito. Comprobar sólo posición puede detener el sistema mientras cruza rápidamente el equilibrio.

33.21 Detectar fallos durante la ejecución

Después de cada paso conviene verificar valores finitos y límites razonables. Un estado enorme puede indicar un paso inadecuado antes de llegar a Infinity.

Ante un fallo, el programa debería detener la evolución y mostrar parámetros, integrador, tiempo y último estado válido. Ocultar el error mediante un reinicio automático dificulta el diagnóstico.

33.22 Laboratorio de integración

Las curvas comparan la misma ecuación y el mismo paso durante diez períodos naturales. La referencia naranja es analítica; las restantes pertenecen a los tres métodos numéricos.

33.23 Rendimiento sin alterar la física

• No crear objetos temporales innecesarios dentro de millones de pasos si el rendimiento es crítico.
• Separar la frecuencia de integración de la cantidad de muestras dibujadas.
• Limitar los pasos por cuadro para evitar una espiral de atraso en tiempo real.
• No aumentar Δt sólo para obtener más cuadros por segundo sin evaluar el error.
• Medir antes de optimizar: en simulaciones pequeñas, la claridad suele ser más valiosa.

33.24 Errores frecuentes y ejercicio

• Actualizar x antes de v y llamarlo semimplícito: ese orden corresponde a Euler explícito.
• Usar el tiempo de cuadro como Δt: vuelve variable la física.
• Modificar el mismo estado al calcular las etapas de RK4: contamina las pendientes.
• Guardar referencias repetidas: destruye el historial.
• Comparar métodos con pasos diferentes: mezcla dos causas de error.
• Confundir amortiguamiento físico con pérdida numérica: deben analizarse por separado.

Con m = 1 kg, k = 25 N/m, b = 1 N·s/m, x₀ = 0,10 m, v₀ = 0 y Δt = 0,02 s, realizá un paso de Euler explícito y uno semimplícito.

const m = 1, k = 25, b = 1, x0 = 0.1, v0 = 0, dt = 0.02;
const a0 = (-k * x0 - b * v0) / m;
const velocidadNueva = v0 + a0 * dt;
const explicito = { x: x0 + v0 * dt, v: velocidadNueva };
const semimplicito = { x: x0 + velocidadNueva * dt, v: velocidadNueva };

console.log({ a0, explicito, semimplicito });

33.25 Ideas para recordar

• Integrar convierte derivadas continuas en actualizaciones discretas.
• La aceleración debe ser una función comprobable de parámetros, estado y tiempo.
• Euler explícito tiende a agregar energía a un oscilador ideal.
• Euler semimplícito mejora mucho el comportamiento con un cambio de orden.
• RK4 usa cuatro evaluaciones para obtener mayor precisión por paso.
• Δt debe compararse con las escalas temporales del sistema.
• La integración, el muestreo y el dibujo son procesos diferentes.
• En tiempo real, un acumulador permite mantener el paso físico fijo.