Biblioteca numérica

01 · Integrar lo aprendido

El objetivo del proyecto

Construiremos el núcleo de una biblioteca para problemas de valor inicial:

La biblioteca debe resolver modelos escalares y vectoriales con Euler, Heun, punto medio y RK4, sin duplicar el bucle de integración.

Un proyecto numérico no termina cuando produce números. También debe explicar qué método usó, cuánto trabajo realizó, cómo trató el intervalo y con qué pruebas se verificó.

02 · Antes de programar

Requisitos funcionales y numéricos

Definimos un alcance concreto:

• estados representados siempre mediante arreglos numéricos;
• intervalos crecientes con paso máximo positivo;
• último paso acortado para terminar exactamente en tFinal;
• métodos explícitos de paso fijo;
• trayectoria completa o solo estado final;
• conteo de pasos y evaluaciones de f;
• detección de dimensiones incompatibles y estados no finitos;
• parada opcional ante un evento definido por el usuario.

Quedan como extensiones el paso adaptativo, localización precisa de eventos, métodos implícitos y matrices dispersas.

03 · Separar responsabilidades

Arquitectura de la biblioteca

• Operaciones vectoriales: combinan estados y derivadas.
• Métodos: implementan un único paso y declaran su costo.
• Integrador: administra tiempo, validaciones, almacenamiento y eventos.
• Modelos: contienen la ecuación y, si existe, su solución exacta.
• Pruebas: verifican casos conocidos, convergencia y comportamiento inválido.

Esta separación evita que una mejora en RK4 obligue a modificar los modelos o la administración de la malla.

04 · Una API predecible

Contratos de entrada y salida

El modelo recibe tiempo y estado, y devuelve una derivada con la misma dimensión:

El solucionador recibe un objeto de opciones:

const problema = {
  f: (t, [x, v]) => [v, -x],
  t0: 0,
  y0: [1, 0],
  tFinal: 2 * Math.PI,
  h: 0.1,
  metodo: "rk4",
  guardarTrayectoria: true
};

const resultadoEsperado = {
  t: "tiempo final alcanzado",
  y: "copia del estado final",
  trayectoria: "arreglo opcional de puntos",
  pasos: "cantidad de pasos aceptados",
  evaluaciones: "cantidad de llamadas a f",
  evento: "información de parada o null"
};

console.log(problema);
console.log(resultadoEsperado);

Usar un objeto evita listas largas de argumentos posicionales y permite agregar opciones sin romper llamadas existentes.

05 · Un núcleo pequeño

Operaciones vectoriales seguras

Las etapas necesitan sumar combinaciones de derivadas sin modificar el estado original:

function validarMismaDimension(...vectores) {
  const dimension = vectores[0].length;
  if (!vectores.every(vector => vector.length === dimension)) {
    throw new Error("Dimensiones incompatibles");
  }
}

function combinar(base, h, terminos) {
  validarMismaDimension(base, ...terminos.map(([, vector]) => vector));
  return base.map((valor, i) =>
    valor + h * terminos.reduce(
      (suma, [peso, vector]) => suma + peso * vector[i], 0));
}

const estado = [1, 2];
const k1 = [0.5, -1];
const k2 = [0.4, -0.8];
console.log(combinar(estado, 0.1, [[0.5, k1], [0.5, k2]]));
console.log("estado original:", estado);

Crear arreglos nuevos simplifica el razonamiento. Una optimización posterior puede reutilizar memoria, pero solo después de contar con pruebas que protejan el resultado.

06 · Métodos intercambiables

Un registro de estrategias de paso

Cada estrategia devuelve el estado siguiente y conoce cuántas evaluaciones realiza:

const sumar = (y, h, terminos) => y.map((valor, i) =>
  valor + h * terminos.reduce((s, [peso, k]) => s + peso * k[i], 0));

const metodos = {
  euler: {
    evaluaciones: 1,
    paso(f, t, y, h) {
      return sumar(y, h, [[1, f(t, y)]]);
    }
  },
  heun: {
    evaluaciones: 2,
    paso(f, t, y, h) {
      const k1 = f(t, y);
      const k2 = f(t + h, sumar(y, h, [[1, k1]]));
      return sumar(y, h, [[1 / 2, k1], [1 / 2, k2]]);
    }
  },
  puntoMedio: {
    evaluaciones: 2,
    paso(f, t, y, h) {
      const k1 = f(t, y);
      const k2 = f(t + h / 2, sumar(y, h, [[1 / 2, k1]]));
      return sumar(y, h, [[1, k2]]);
    }
  },
  rk4: {
    evaluaciones: 4,
    paso(f, t, y, h) {
      const k1 = f(t, y);
      const k2 = f(t + h / 2, sumar(y, h, [[1 / 2, k1]]));
      const k3 = f(t + h / 2, sumar(y, h, [[1 / 2, k2]]));
      const k4 = f(t + h, sumar(y, h, [[1, k3]]));
      return sumar(y, h, [[1 / 6, k1], [1 / 3, k2], [1 / 3, k3], [1 / 6, k4]]);
    }
  }
};

console.log(Object.keys(metodos));
console.log("RK4 usa", metodos.rk4.evaluaciones, "evaluaciones por paso");

07 · Orquestar la simulación

El bucle de integración

El integrador no conoce las fórmulas internas del método. Solo valida, solicita pasos y reúne metadatos:

function integrar({ f, t0, y0, tFinal, h, metodo, guardar = true }) {
  if (typeof f !== "function") throw new TypeError("f debe ser una función");
  if (!Array.isArray(y0) || y0.length === 0) throw new TypeError("y0 inválido");
  if (!(h > 0) || tFinal < t0) throw new RangeError("Intervalo o paso inválido");
  if (!metodo?.paso || !(metodo.evaluaciones > 0)) throw new Error("Método inválido");

  let t = t0;
  let y = [...y0];
  let pasos = 0;
  const trayectoria = guardar ? [{ t, y: [...y] }] : undefined;

  while (t < tFinal) {
    const pasoActual = Math.min(h, tFinal - t);
    const siguiente = metodo.paso(f, t, y, pasoActual);
    if (siguiente.length !== y.length || !siguiente.every(Number.isFinite)) {
      throw new Error("El método produjo un estado inválido");
    }
    t += pasoActual;
    y = siguiente;
    pasos++;
    if (guardar) trayectoria.push({ t, y: [...y] });
  }
  return {
    t, y: [...y], trayectoria, pasos,
    evaluaciones: pasos * metodo.evaluaciones
  };
}

const euler = {
  evaluaciones: 1,
  paso: (f, t, y, h) => {
    const k = f(t, y);
    return y.map((valor, i) => valor + h * k[i]);
  }
};
const resultado = integrar({
  f: (t, [y]) => [y], t0: 0, y0: [1], tFinal: 1,
  h: 0.25, metodo: euler
});
console.log(resultado);

El último paso acortado evita sobrepasar el extremo. Copiar estados impide que el usuario modifique accidentalmente el historial al conservar una referencia.

08 · Fallar de forma informativa

Validaciones y errores

Las entradas inválidas deben detener el cálculo cerca de su origen:

• h≤0, intervalo invertido o valores no finitos;
• estado vacío o componentes no numéricas;
• derivada con dimensión diferente;
• nombre de método inexistente;
• límite máximo de pasos excedido;
• estado que se vuelve NaN o infinito.

Conviene incluir en el error el tiempo, número de paso y método. Esa información transforma un «resultado roto» en un problema diagnosticable.

09 · Detenerse por una condición

Eventos y salidas

Una simulación puede terminar cuando una cantidad cruza cero, alcanza un umbral o sale de un dominio. Una función de evento devuelve un valor g(t,𝐲); un cambio de signo indica un cruce.

La versión mínima detiene la integración en el primer nodo posterior al cruce. Una biblioteca más precisa localiza el evento dentro del paso mediante interpolación y búsqueda de raíces.

Separá tres decisiones:

• pasos internos usados para calcular;
• tiempos de salida solicitados por el usuario;
• eventos que pueden cortar o modificar la simulación.

Guardar más puntos no aumenta la precisión de una trayectoria ya calculada.

10 · Verificación básica

Problemas con solución conocida

Una batería pequeña debe incluir comportamientos diferentes:

Además de comparar el estado final, verificá el tiempo alcanzado, cantidad de pasos, dimensión, inmutabilidad del estado inicial y costo declarado.

11 · Probar el orden

Una prueba de refinamiento automatizada

Si E(h)≈Chᵖ, entonces:

function rk4Escalar(f, t, y, h) {
  const k1 = f(t, y);
  const k2 = f(t + h / 2, y + h * k1 / 2);
  const k3 = f(t + h / 2, y + h * k2 / 2);
  const k4 = f(t + h, y + h * k3);
  return y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
}

function resolver(pasos) {
  const h = 1 / pasos;
  let t = 0;
  let y = 1;
  for (let n = 0; n < pasos; n++) {
    y = rk4Escalar((tiempo, estado) => estado, t, y, h);
    t += h;
  }
  return y;
}

const errores = [4, 8, 16, 32].map(pasos =>
  Math.abs(Math.E - resolver(pasos)));
for (let i = 1; i < errores.length; i++) {
  const orden = Math.log2(errores[i - 1] / errores[i]);
  console.log(`refinamiento ${i}: orden observado ${orden.toFixed(4)}`);
}

La prueba no debería exigir exactamente cuatro. Una tolerancia razonable admite que pasos gruesos aún no estén en el régimen asintótico y que pasos diminutos sufran redondeo.

12 · Pruebas cualitativas

Estabilidad e invariantes

Las pruebas exactas no cubren todo. También conviene verificar propiedades:

• un decaimiento estable no debe crecer por elección incorrecta del paso;
• poblaciones y concentraciones deberían respetar positividad cuando corresponda;
• la energía de un oscilador permite medir deriva a largo plazo;
• un problema rígido debe advertir la limitación de métodos explícitos;
• resultados con h y h/2 deben ser coherentes con el orden.

Una biblioteca general no puede garantizar todas las propiedades físicas. Debe exponer suficientes métricas para que el usuario pueda evaluarlas.

13 · Optimizar con evidencia

Costo, memoria y reproducibilidad

Las métricas mínimas son pasos aceptados, pasos rechazados, evaluaciones de f y tiempo total. Si se agrega un método implícito, también cuentan iteraciones y resoluciones lineales.

• Permití desactivar el historial cuando solo interesa el estado final.
• Guardá resultados en tiempos seleccionados para trayectorias largas.
• Reutilizá arreglos solo después de medir y probar aliasing.
• No compares métodos únicamente por cantidad de pasos.
• Registrá configuración, versión y tolerancias para repetir experimentos.

Una optimización que cambia silenciosamente el resultado necesita una justificación numérica, no solo una mejora de velocidad.

14 · Banco de pruebas interactivo

Laboratorio: una API, varios modelos y métodos

El laboratorio utiliza estados vectoriales y el mismo integrador para todos los casos. La curva dorada es exacta y la turquesa es la salida del método seleccionado.

Preguntas para explorar

1. Duplicá los pasos y compará cómo cambia el error de cada método.
2. Buscá métodos que coincidan en crecimiento pero difieran en el modelo logístico.
3. En el oscilador, observá la amplitud después de una vuelta completa.
4. En dos escalas, reducí los pasos hasta que un método explícito se vuelva inestable.

15 · Comprensión de diseño

Confusiones frecuentes

16 · Extender el proyecto

Ejercicios y mejoras

Ejercicio 1: estado final sin historial

Modificá el integrador para no crear trayectoria cuando guardar=false y escribí una prueba que lo confirme.

Ejercicio 2: evento terminal

Agregá un evento que detenga una caída vertical cuando la altura cruce cero.

Ejercicio 3: prueba parametrizada de orden

Creá una función que reciba método, orden esperado y una lista de refinamientos.

Ejercicio 4: salida CSV

Convertí una trayectoria vectorial a columnas t,y0,y1,….

Ejercicio 5: paso adaptativo

Diseñá una extensión que compare un paso de tamaño h con dos de h/2.

17 · Cierre del recorrido

De fórmulas a software numérico confiable

El resultado del proyecto no es solo un solucionador: es una base extensible para experimentar con interpolación, cuadratura, paso adaptativo, métodos implícitos y nuevos modelos con criterios numéricos verificables.