Funciones, integradores y pruebas numéricas

Contrato de una función numérica

Un integrador necesita evaluar

La función debe recibir tiempo, estado y parámetros, y devolver un vector de la misma forma que el estado. La sección oficial sobre funciones explica argumentos, valores de retorno y anotaciones. El contrato científico añade forma, tipo, dominio y comportamiento ante entradas inválidas.

from collections.abc import Callable
import numpy as np
from numpy.typing import NDArray

Vector = NDArray[np.float64]
Field = Callable[[float, Vector, object], Vector]

def evaluate_field(fun: Field, t: float, state: Vector, params) -> Vector:
    state = np.asarray(state, dtype=np.float64)
    slope = np.asarray(fun(t, state, params), dtype=np.float64)
    if slope.shape != state.shape:
        raise ValueError(
            f"El campo devolvió {slope.shape}; se esperaba {state.shape}"
        )
    if not np.isfinite(slope).all():
        raise FloatingPointError("El campo devolvió NaN o infinito")
    return slope

La anotación ayuda a leer y revisar; no valida en tiempo de ejecución. Las comprobaciones sí. El error se lanza cuando aparece la inconsistencia, no varias operaciones después. La documentación oficial sobre excepciones sirve para distinguir un valor inválido de un fallo numérico.

En evaluate_field(fun, t, state, params), fun también es una función: Python permite recibirla como argumento y llamarla dentro con fun(t, state, params). Esta composición se construye aquí desde cero. Primero se entienden y prueban evaluate_field, euler_step e integrate_fixed; solo después se llama a un integrador validado de SciPy. La comparación permite distinguir el algoritmo que el estudiante programó de la interfaz pública proporcionada por una biblioteca.

Aproximación por el método de Euler

La expansión de Taylor da

Sustituir $\dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x})$ y omitir el resto produce Euler:

Su error local es $O(h^2)$ y el error global, bajo hipótesis regulares en un intervalo fijo, es $O(h)$. La implementación de un paso conserva la entrada:

def euler_step(fun, t, state, h, params):
    if not np.isfinite(h) or h <= 0.0:
        raise ValueError("h debe ser positivo y finito")
    state = np.asarray(state, dtype=np.float64)
    slope = evaluate_field(fun, t, state, params)
    candidate = state + h * slope
    if not np.isfinite(candidate).all():
        raise FloatingPointError("Euler produjo NaN o infinito")
    return candidate

Una función de paso no decide cuántos pasos ejecutar ni dónde guardar. Esa responsabilidad pertenece al conductor:

def integrate_fixed(step, fun, t_span, y0, h, params):
    t0, tf = map(float, t_span)
    if tf <= t0:
        raise ValueError("El intervalo debe avanzar")
    count_float = (tf - t0) / h
    count = round(count_float)
    if not np.isclose(count_float, count):
        raise ValueError("El intervalo no contiene un número entero de pasos")

    times = t0 + h * np.arange(count + 1)
    states = np.empty((count + 1, np.size(y0)), dtype=np.float64)
    states[0] = y0
    for index in range(count):
        states[index + 1] = step(
            fun, times[index], states[index], h, params
        )
    return times, states

Rechazar un último paso fraccionario simplifica esta versión pedagógica. Un integrador general puede ajustar el último paso, pero debe hacerlo de manera explícita.

Práctica interactiva. Ejecuta en una celda las funciones evaluate_field, euler_step e integrate_fixed. Define después decay(t, state, rate) para $x'=-rx$ y calcula hasta $t=1$ con $h=0.2$, $0.1$ y $0.05$. Compara el último valor con $e^{-r}$, registra el error y comprueba si disminuye aproximadamente a la mitad. El ejercicio usa Python y NumPy; no abre archivos ni necesita hardware externo.

Cálculo manual de un paso

La ecuación $x'=-2x$ permite revisar la aritmética antes de confiar en un conductor completo. Con $x_0=1$ y $h=0.1$, Euler usa la pendiente inicial $-2$ y obtiene $x_1=1+0.1(-2)=0.8$. La solución exacta en ese instante es $e^{-0.2}\approx0.8187308$; el error del primer paso es, por tanto, aproximadamente $1.87\times10^{-2}$. Este valor prueba simultáneamente el signo del campo, el tamaño del paso y la actualización del estado.

Si se llega al mismo tiempo mediante dos pasos de longitud $0.05$, cada paso multiplica por $0.9$ y el resultado es $0.81$. El error baja a $8.73\times10^{-3}$, cerca de la mitad pero no exactamente. La relación de orden es asintótica: se vuelve clara al refinar una sucesión de mallas, no por una igualdad rígida en cualquier $h$. Una prueba sensata calcula varios errores y acepta un intervalo para el orden observado cuando ya se distingue el régimen de truncamiento y todavía no domina el redondeo.

Este ejemplo también descubre la diferencia entre actualizar desde el estado anterior y hacerlo desde uno ya modificado. Si una función cambia state en sitio, las pendientes intermedias de RK4 dejan de partir de los estados definidos por su fórmula. Por eso las funciones de paso devuelven un arreglo nuevo y una prueba conserva una copia de la entrada para compararla después de la llamada. La corrección de una fórmula matemática y la ausencia de efectos laterales son dos propiedades distintas; ambas requieren prueba.

Etapas internas del método RK4

El método clásico calcula:

y actualiza

def rk4_step(fun, t, state, h, params):
    state = np.asarray(state, dtype=np.float64)
    k1 = evaluate_field(fun, t, state, params)
    k2 = evaluate_field(fun, t + h/2, state + h*k1/2, params)
    k3 = evaluate_field(fun, t + h/2, state + h*k2/2, params)
    k4 = evaluate_field(fun, t + h, state + h*k3, params)
    candidate = state + h*(k1 + 2*k2 + 2*k3 + k4)/6
    if not np.isfinite(candidate).all():
        raise FloatingPointError("RK4 produjo NaN o infinito")
    return candidate

RK4 tiene orden global cuatro en problemas suaves, pero no es automáticamente adecuado para rigidez, discontinuidades o cualquier paso. El orden se verifica en el código concreto antes de usarlo como referencia.

Orden de convergencia y estabilidad

Un método puede exhibir el orden esperado y aun ser inútil con un paso grande. Para la ecuación de prueba $x'=-\lambda x$, Euler produce $x_{n+1}=(1-h\lambda)x_n$. Si $\lambda>0$, la solución exacta decae, pero la aproximación solo lo hace sin crecimiento cuando $|1-h\lambda|\leq1$. En el caso real eso exige $0\leq h\lambda\leq2$. Con $\lambda=100$ y $h=0.05$, el factor es $-4$: la secuencia alterna y su amplitud crece aunque la ecuación describa una relajación rápida.

Reducir el paso a $0.01$ estabiliza Euler en ese ejemplo, pero obliga a resolver la escala rápida aun si el interés está en una evolución mucho más lenta. Esa separación de escalas caracteriza la dificultad práctica de muchos problemas rígidos. RK4 posee una región de estabilidad mayor que Euler, no ilimitada. Un método explícito de alto orden no sustituye automáticamente a un método diseñado para rigidez.

Una prueba de orden y una prueba de estabilidad responden entonces preguntas distintas. La primera refina $h$ en un intervalo fijo y observa cómo decrece el error. La segunda varía $h\lambda$ y comprueba si una solución que debe decaer se mantiene acotada. En una aplicación científica se documentan ambas cuando las escalas del modelo puedan comprometer el paso. Declarar solo que “RK4 es más preciso” oculta esta condición.

Pruebas con una solución conocida

Para $x'=-2x$, $x(0)=1$, la solución es $x(t)=e^{-2t}$. El error final es

Si un método tiene orden $p$,

def decay(t, state, rate):
    return -rate * state

def final_error(step, h):
    _, states = integrate_fixed(step, decay, (0.0, 1.0), [1.0], h, 2.0)
    return abs(states[-1, 0] - np.exp(-2.0))

errors = np.array([final_error(euler_step, h)
                   for h in (0.2, 0.1, 0.05, 0.025)])
orders = np.log2(errors[:-1] / errors[1:])

La prueba no exige orden exacto en mallas gruesas ni cuando el redondeo domina. Define una ventana razonable y comprueba tendencia. Alligood, Sauer y Yorke usan problemas controlados para estudiar error de integración Alligood, 1996 Kathleen T. Alligood, Tim D. Sauer y James A. Yorke (1996) Chaos: An Introduction to Dynamical Systems Springer Ubicación consultada: apéndice B, esp. §B.2 . El libro aporta el experimento numérico; Python implementa la prueba.

También se prueban invariantes. Para el oscilador armónico $q'=p$, $p'=-q$, la energía exacta

es constante. Medir deriva de energía durante cien periodos descubre errores que una comparación del último estado puede ocultar. El error de fase se presenta por separado del error de amplitud.

Comparación verificable con SciPy

Después de comprender el paso fijo se usa la documentación oficial de solve_ivp. Su contrato incluye función, intervalo, estado inicial, método, tolerancias, tiempos de salida, eventos y salida densa:

solve_ivp no es una función escrita en esta unidad. Se importa desde la interfaz pública scipy.integrate. Sus tres primeros argumentos nombran el campo, el intervalo y el estado inicial. Los argumentos method, rtol y atol seleccionan el método y sus controles; no modifican las ecuaciones del modelo.

from scipy.integrate import solve_ivp

solution = solve_ivp(
    fun=lambda t, y: lorenz(t, y, params),
    t_span=(0.0, 80.0),
    y0=(1.0, 1.0, 1.0),
    method="DOP853",
    rtol=1e-9,
    atol=1e-11,
    dense_output=False,
)

if not solution.success:
    raise RuntimeError(solution.message)
states = solution.y.T
times = solution.t

SciPy devuelve variables por filas en solution.y; transponer una vez adapta la convención $(N,d)$ del curso. Las tolerancias controlan error local estimado, no garantizan una trayectoria exacta a largo plazo. Se repite con tolerancias menores y se comparan observables compatibles con sensibilidad.

Para detectar un cruce $x=0$ con orientación positiva:

def crossing(t, y):
    return y[0]

crossing.direction = 1
crossing.terminal = False

Pasar el evento al solucionador localiza el cruce mediante interpolación. Tomar la muestra de tiempo más cercana engrosa artificialmente una sección de Poincaré.

Tolerancias ligadas a la escala de cada variable

En un integrador adaptativo, la aceptación de un paso compara un error estimado con una escala semejante a atol + rtol * abs(y). La tolerancia relativa domina cuando una componente es grande; la absoluta establece el piso cerca de cero. Elegir rtol=1e-9 no significa obtener nueve cifras correctas en toda la trayectoria ni limita directamente el error global. Significa que el controlador local persigue un umbral definido por ambas tolerancias.

Una única tolerancia absoluta puede ser inadecuada cuando las variables usan escalas distintas. En un circuito, por ejemplo, un voltaje puede ser del orden de unidades y una corriente del orden de milésimas. solve_ivp acepta un arreglo de tolerancias absolutas; cada valor debe elegirse respecto a la resolución científicamente relevante de su componente, no para forzar que el programa termine. Reescalar las variables es otra posibilidad, pero la transformación y sus unidades tienen que permanecer explícitas.

El argumento t_eval decide dónde se reporta la solución, no obliga al método adaptativo a tomar exactamente esos pasos internos. Esto permite comparar dos corridas sobre una malla común sin confundir muestreo con control de error. Si se solicita dense_output=True, la interpolación continua construida por el método puede evaluarse después; no debe reemplazarse por una interpolación lineal improvisada sin medir su efecto en máximos o cruces.

Además del arreglo de estados se inspeccionan success, status, message y el número de evaluaciones del campo. Un aumento grande de evaluaciones al estrechar tolerancias puede señalar una región difícil, un evento frecuente o rigidez. El parámetro max_step puede ser necesario para no saltar una escala temporal relevante, pero imponerlo demasiado pequeño elimina la ventaja adaptativa. Cada ajuste se relaciona con una propiedad del modelo y queda registrado junto con el resultado.

Para Lorenz, dos tolerancias estrechas pueden coincidir al inicio y separarse más tarde por sensibilidad. Esa divergencia no invalida por sí sola al integrador. Se verifica el corto plazo con refinamiento, se comparan tiempos de eventos mientras sean estables y, para horizontes largos, se estudian observables como distribuciones, retornos o exponentes. La prueba numérica debe corresponder al tipo de afirmación que se hará con la trayectoria.

Fallos explícitos y validación de entradas

Se prueban entradas inválidas: paso no positivo, intervalo invertido, estado vacío, forma incorrecta y campo no finito. También se decide qué hacer si el integrador no llega al tiempo final. Atrapar toda excepción y continuar crea datos parciales que parecen completos. Solo se captura un error cuando existe una recuperación definida y registrada.

Una prueba de regresión congela una propiedad justificada, no un número elegido porque “era lo que salía”. Puede guardar el resultado de una solución exacta, el orden observado dentro de tolerancia o una identidad como la suma de exponentes. Si se actualiza el valor esperado, el cambio se explica.

Todavía no se ha ejecutado el programa.

Pruebas del integrador

1. Implementa RK4 sin modificar el estado de entrada.
2. Recupera los órdenes observados de Euler y RK4 en una ventana adecuada.
3. Mide energía y fase del oscilador durante cien periodos.
4. Provoca un error de forma y comprueba que el mensaje indique la causa.
5. Compara dos tolerancias de SciPy y registra éxito, evaluaciones y mensaje.

El conjunto se considera verificable si los casos exactos recuperan su orden, las entradas inválidas fallan cerca de la causa y las corridas adaptativas conservan su estado y diagnóstico completos. Solo entonces una dinámica sensible añade dificultad científica, no incertidumbre básica sobre el código.