Simulaciones individuales y familias de simulaciones

En los capítulos anteriores trabajamos con parámetros elegidos de antemano. Una simulación individual responde qué hace un modelo bajo una combinación concreta de ecuaciones, parámetros, condición inicial, horizonte temporal y método numérico. En las ciencias de la complejidad, sin embargo, rara vez basta con conocer una sola trayectoria. Interesa saber cómo cambia el comportamiento cuando se modifican las condiciones que definen el problema.

Un barrido paramétrico ejecuta sistemáticamente el modelo para múltiples valores de uno o más parámetros y registra resultados definidos de antemano. La misma lógica se aplica a las condiciones iniciales, a la intensidad de una perturbación o a la topología de una red, entre otras cantidades que forman parte de la especificación del experimento computacional.

La operación parece sencilla, porque basta recorrer una cuadrícula y guardar las salidas. El problema científico comienza precisamente allí. Una cuadrícula finita no revela automáticamente la estructura del espacio de parámetros. Puede pasar por alto regiones estrechas o confundir un transitorio prolongado con un estado estacionario. También puede sugerir una transición donde solo existe un criterio arbitrario de clasificación. Si dos simulaciones vecinas producen resultados diferentes, todavía hay que determinar si la diferencia corresponde al modelo, al horizonte de observación, a la condición inicial o al procedimiento numérico.

Por eso un barrido debe definirse por algo más que los valores que prueba. Hay que especificar qué parámetro se modifica y qué condiciones se mantienen constantes, además de las condiciones iniciales consideradas y del tiempo de integración. También hay que declarar qué eventos detienen la simulación y qué resolución se usa, junto con la lista de magnitudes que se almacenan. Si cambia cualquiera de esas decisiones, el mapa obtenido puede cambiar.

En este capítulo construiremos un ejemplo cuya estructura matemática conocemos de antemano. Esa ventaja permitirá comprobar qué descubre un barrido y qué puede ocultar. Después trasladaremos las conclusiones a modelos en los que la respuesta exacta no está disponible.

Un modelo de crecimiento limitado con extracción constante

Consideremos el modelo adimensional

dxdt=x(1−x)−h,\frac{dx}{dt}=x(1-x)-h,

donde x(t)x(t) es una cantidad normalizada y h≥0h\geq0, una extracción constante. Puede leerse como un modelo mínimo de una población sometida a cosecha, aunque aquí interesa sobre todo como sistema dinámico de una variable con una estructura paramétrica conocida (Strogatz, 2015).

El término x(1−x)x(1-x) corresponde a un crecimiento logístico con capacidad de carga y tasa intrínseca normalizadas a uno. La extracción hh no depende del tamaño de la población. Esa suposición es deliberadamente simplificadora y pierde plausibilidad cuando xx se acerca a cero. Por eso trataremos el cruce de x=0x=0 como un evento terminal, en el que la integración se detiene antes de que el modelo produzca poblaciones negativas sin interpretación física.

Los equilibrios satisfacen x(1−x)−h=0x(1-x)-h=0, es decir, x2−x+h=0x^2-x+h=0. Cuando h≤1/4h\leq1/4, las soluciones son

xu(h)=1−1−4h2,xs(h)=1+1−4h2.x_u(h)=\frac{1-\sqrt{1-4h}}{2}, \qquad x_s(h)=\frac{1+\sqrt{1-4h}}{2}.

La derivada del campo es fx(x;h)=1−2xf_x(x;h)=1-2x. Para h<1/4h<1/4, el equilibrio inferior, xu<1/2x_u<1/2, tiene fx>0f_x>0 y es inestable, mientras que el superior, xs>1/2x_s>1/2, tiene fx<0f_x<0 y es localmente asintóticamente estable. En h=1/4h=1/4 ambos coinciden en x=1/2x=1/2. Para h>1/4h>1/4 no hay equilibrios reales y x(1−x)−h<0x(1-x)-h<0 para todo xx, de modo que la trayectoria decrece hasta alcanzar cero.

El valor hc=1/4h_c=1/4 es, por tanto, un umbral estructural del modelo. Más adelante estudiaremos formalmente esta situación como una bifurcación silla-nodo (Kuznetsov, 2004; Strogatz, 2015). Aquí la usaremos como patrón de referencia para evaluar un barrido numérico.

La estructura depende también de la condición inicial. Para h<1/4h<1/4, una trayectoria que comienza por encima del equilibrio inestable xu(h)x_u(h) se aproxima al equilibrio estable xs(h)x_s(h). Si comienza por debajo de xu(h)x_u(h), el campo apunta hacia valores menores y la trayectoria alcanza el evento x=0x=0. La frontera entre ambos resultados depende, así, del par (h,x0)(h,x_0).

Trayectorias desde x₀ = 1 para seis valores de h.
Figura 28.1. Trayectorias desde x₀ = 1 para seis valores de h. Por debajo del umbral h = 1/4, cada trayectoria se acerca a un equilibrio estable. Por encima del umbral, la trayectoria alcanza x = 0 en un tiempo finito y la integración se detiene. Ampliar figura ↗

Diseñar el barrido antes de ejecutarlo

Supongamos que ignoramos el valor analítico hc=1/4h_c=1/4 y queremos explorar el modelo mediante simulación. Una estrategia inicial podría recorrer

h=0,10; 0,12; 0,14; …; 0,40,h=0{,}10;\ 0{,}12;\ 0{,}14;\ \ldots;\ 0{,}40,

integrar desde x0=1x_0=1 hasta T=80T=80 y registrar si la trayectoria alcanza x=0x=0. La estructura computacional puede escribirse así.

for h in malla:
    solucion = integrar(modelo, x0=1.0, parametro=h, T=80)
    guardar(h, evento, tiempo_evento, estado_final)

El pseudocódigo parece trivial, pero contiene varias decisiones científicas. El horizonte T=80T=80 puede bastar para unas trayectorias y quedar corto para otras. El estado final solo tiene significado como resumen si se sabe qué ocurrió antes. La ausencia de evento no implica necesariamente convergencia a un equilibrio. Incluso la variable almacenada condiciona qué patrones podrán reconocerse después.

Un barrido reproducible separa, por esa razón, tres clases de información. Las entradas del experimento incluyen los parámetros, la condición inicial, las tolerancias, el método, el horizonte y las reglas de evento. Los resultados primarios son las cantidades definidas antes del análisis, como la ocurrencia y el tiempo de los eventos, el estado final o los valores extremos. Los metadatos del cálculo registran el modo de terminación del integrador y el número de evaluaciones, junto con las advertencias emitidas.

Guardar solo una etiqueta como «persiste» o «colapsa» elimina información que podría hacer falta para reinterpretar el resultado. Es más seguro conservar las cantidades a partir de las cuales se construyó la clasificación.

Cuadrícula gruesa y refinamiento local

Con la malla anterior, la última simulación sin evento ocurre en h=0,24h=0{,}24 y la primera con evento en h=0,26h=0{,}26. Esos resultados permiten afirmar que la transición observada por el barrido está entre ambos valores, bajo la condición inicial y el horizonte usados. Todavía no puede concluirse que el umbral sea 0,250{,}25. Esa cifra podría obtenerse por interpolación, pero la interpolación no reemplaza una nueva evaluación del modelo. Lo correcto es usar el resultado grueso para decidir dónde aumentar la resolución.

Una segunda malla puede recorrer h=0,235; 0,236; …; 0,270h=0{,}235;\ 0{,}236;\ \ldots;\ 0{,}270. La búsqueda refinada evita gastar el mismo número de simulaciones en regiones donde el comportamiento cambia poco y en regiones donde una pequeña diferencia de parámetros modifica el resultado cualitativo. Con el mismo horizonte T=80T=80, esa malla ubica el cambio entre 0,2510{,}251 y 0,2520{,}252, aunque el umbral exacto es 0,250{,}25 (figura 28.2). En h=0,251h=0{,}251 la trayectoria no alcanza x=0x=0 hasta t≈95,35t\approx95{,}35, después del horizonte, de modo que el barrido la clasifica sin evento. Más adelante veremos por qué el horizonte pesa tanto cerca del umbral.

Clasificación de las simulaciones de las mallas gruesa y fina según si alcanzan x = 0 antes de T = 80.
Figura 28.2. Clasificación de las simulaciones de las mallas gruesa y fina según si alcanzan x = 0 antes de T = 80. La malla gruesa sitúa el cambio entre 0,24 y 0,26. La malla fina lo sitúa entre 0,251 y 0,252, porque cerca del umbral exacto h = 1/4 el evento ocurre después del horizonte. Ampliar figura ↗

El refinamiento tampoco convierte el barrido en una demostración de bifurcación. Un barrido ordinario observa únicamente las soluciones que alcanza desde las condiciones iniciales especificadas. Puede omitir ramas inestables o estados con cuencas pequeñas. Los métodos de continuación numérica están diseñados para seguir familias de soluciones y localizar cambios estructurales con herramientas más apropiadas que una cuadrícula (Govaerts, 2000; Uecker, 2022).

Un barrido es un experimento exploratorio sobre el comportamiento calculado. Una rama de equilibrios y su estabilidad son objetos matemáticos que requieren análisis adicional.

Eventos como condiciones matemáticas

En muchos modelos interesa el momento en que ocurre una condición determinada. Un evento puede definirse mediante una función g(t,x)=0g(t,x)=0, que en nuestro ejemplo es g(t,x)=xg(t,x)=x. El evento ocurre cuando la trayectoria cruza x=0x=0. Como solo interesa el cruce descendente, se especifica una dirección negativa. Como continuar más allá carece de interpretación en el modelo poblacional, el evento se declara terminal. En solve_ivp de SciPy, una implementación mínima es la siguiente.

from scipy.integrate import solve_ivp


def rhs(t, y, h):
    x = y[0]
    return [x * (1 - x) - h]


def extincion(t, y, h):
    return y[0]


extincion.terminal = True
extincion.direction = -1

sol = solve_ivp(
    rhs,
    (0, 100),
    [1.0],
    args=(0.30,),
    events=extincion,
    rtol=1e-9,
    atol=1e-11,
)

t_evento = sol.t_events[0][0]

El integrador no necesita que uno de sus puntos almacenados coincida exactamente con el cero. Detecta un cambio de signo de la función de evento dentro de un paso y usa un procedimiento de búsqueda de raíces sobre la solución interpolada para estimar el instante correspondiente (Hairer et al., 1993). La documentación de SciPy advierte, sin embargo, que el mecanismo se apoya en cambios de signo entre pasos, de modo que pueden perderse eventos si la función cruza cero varias veces dentro de un solo paso (SciPy Developers, s. f.).

Definir un evento con precisión matemática no elimina la dependencia respecto del procedimiento numérico. Hay que comprobar que la escala temporal del integrador permite observar la estructura relevante. En ello influyen las tolerancias y el paso máximo permitido, max_step.

Tampoco todo umbral científico debe implementarse como un evento. Si un criterio clínico o experimental se define mediante una medición ruidosa, tratarlo como una superficie exacta g=0g=0 puede ser una idealización demasiado fuerte. El formalismo es correcto solo respecto de la definición introducida en el modelo.

Verificar un evento mediante una solución exacta

El ejemplo tiene una ventaja poco frecuente, porque para h>1/4h>1/4 puede calcularse exactamente el tiempo que tarda una trayectoria en ir desde x0=1x_0=1 hasta cero. Definamos δ=h−1/4>0\delta=h-1/4>0 y escribamos

x(1−x)−h=−[(x−12)2+δ].x(1-x)-h=-\left[\left(x-\frac12\right)^2+\delta\right].

Entonces dt=−dx/[(x−1/2)2+δ]dt=-dx/\left[(x-1/2)^2+\delta\right] y, al integrar desde x=1x=1 hasta x=0x=0,

t0(h)=∫01dx(x−1/2)2+δ=2δarctan⁡(12δ).t_0(h)=\int_0^1 \frac{dx}{(x-1/2)^2+\delta} =\frac{2}{\sqrt{\delta}} \arctan\left(\frac{1}{2\sqrt{\delta}}\right).

Para h=0,30h=0{,}30, por ejemplo, δ=0,05\delta=0{,}05 y t0≈10,288256t_0\approx10{,}288256. El bloque anterior, ejecutado con SciPy 1.17.1, da 10,2882560310{,}28825603, con un error absoluto de unos 1×10−81\times10^{-8}. Para h=0,251h=0{,}251, muy cerca del umbral, el tiempo sube hasta unos 95,351295{,}3512.

Cuando hh se acerca a 1/41/4 desde valores mayores, arctan⁡(1/(2δ))→π/2\arctan\bigl(1/(2\sqrt\delta)\bigr)\to\pi/2, de modo que

t0(h)∼πh−1/4.t_0(h)\sim\frac{\pi}{\sqrt{h-1/4}}.

El tiempo del evento diverge al acercarse al umbral desde arriba (figura 28.3). Esa propiedad explica por qué un horizonte fijo puede ocultar la transición. Una trayectoria puede terminar la simulación sin haber alcanzado el evento porque todavía atraviesa una región de evolución extremadamente lenta. Con T=80T=80, el evento queda fuera del horizonte para todo hh menor que 0,25140{,}2514, el valor que resuelve t0(h)=80t_0(h)=80.

Tiempo exacto hasta el evento x = 0 para h > 1/4, junto con los tiempos que detecta la malla gruesa.
Figura 28.3. Tiempo exacto hasta el evento x = 0 para h > 1/4, junto con los tiempos que detecta la malla gruesa. La divergencia cerca del umbral muestra por qué un horizonte finito, aquí T = 80, puede ocultar un resultado que ocurre más tarde. Ampliar figura ↗

El horizonte temporal también es un parámetro del experimento computacional

El mismo fenómeno aparece desde el lado estable. Cuando h<1/4h<1/4, la tasa local de recuperación alrededor de xs(h)x_s(h) es fx(xs)=1−2xs=−1−4hf_x(x_s)=1-2x_s=-\sqrt{1-4h}, de modo que una escala temporal local característica es

τ(h)=11−4h.\tau(h)=\frac{1}{\sqrt{1-4h}}.

La escala crece cuando hh se acerca a 1/41/4 desde valores menores. Para h=0,2499h=0{,}2499, τ≈50\tau\approx50. Desde x0=1x_0=1 el equilibrio estable es xs=0,51x_s=0{,}51, pero en T=20T=20 la trayectoria todavía está en unos 0,5460{,}546. En T=80T=80 sigue en unos 0,5150{,}515. Acercarse a menos de 0,0010{,}001 del equilibrio exige llegar a t≈150t\approx150.

Por eso, clasificar el «estado estacionario» con el último valor calculado puede ser engañoso. Cerca de una transición, el tiempo necesario para acercarse al régimen asintótico puede aumentar mucho.

Una práctica habitual consiste en descartar un transitorio y calcular estadísticas solo después de cierto tiempo. Esa decisión es útil, aunque requiere justificación. Un descarte fijo de 100 unidades puede ser excesivo en una región del espacio de parámetros e insuficiente en otra. En modelos oscilatorios o caóticos, además, decidir qué parte es transitoria puede ser más difícil que en este ejemplo de una variable.

Integrar «hasta que parezca estable» no resuelve el problema. Hace falta un criterio operacional, como la distancia a un equilibrio conocido, el tamaño de la derivada, la estabilidad de una estadística en ventanas sucesivas o la detección de un evento. Cada criterio responde a una pregunta diferente.

Barrer parámetros sin barrer condiciones iniciales puede ocultar la estructura

Hasta ahora usamos x0=1x_0=1. Si el experimento se repite con distintas condiciones iniciales, aparece otra dimensión del problema. Para h<1/4h<1/4, el equilibrio inferior xu(h)x_u(h) es inestable y separa dos resultados en este sistema de una variable. Por encima de xu(h)x_u(h) la trayectoria se dirige hacia xs(h)x_s(h), mientras que por debajo alcanza cero.

Una afirmación como «con h=0,20h=0{,}20 la población persiste» está, por tanto, incompleta. Con ese parámetro hay condiciones iniciales que persisten y otras que no, porque el resultado depende del estado desde el cual comienza la trayectoria. Con h=0,20h=0{,}20, la frontera está en xu≈0,2764x_u\approx0{,}2764.

La figura 28.4 muestra un barrido en dos dimensiones, con hh en el eje horizontal y x0x_0 en el vertical. La curva analítica xu(h)x_u(h) separa las regiones de los dos comportamientos cuando h<1/4h<1/4. Por encima de h=1/4h=1/4 el modelo no tiene equilibrio positivo y toda condición inicial positiva termina alcanzando cero.

Mapa del evento de extinción según extracción y condición inicial.
Figura 28.4. Barrido conjunto de la extracción h y la condición inicial x₀, con horizonte T = 400. Las celdas oscuras corresponden a trayectorias que alcanzan x = 0. Para h < 1/4, la frontera entre resultados sigue al equilibrio inestable x_u(h), con el escalonamiento que impone la resolución de la malla. Ampliar figura ↗

En dimensiones mayores, las fronteras entre comportamientos pueden formar superficies complejas para las que normalmente no se dispone de una fórmula analítica. Un barrido de condiciones iniciales puede dar entonces una aproximación empírica de las cuencas de atracción, aunque su resolución vuelve a importar. Una cuenca estrecha puede desaparecer visualmente si la cuadrícula es demasiado gruesa.

La clasificación de un punto muy cercano a una frontera puede ser, además, sensible al error numérico, lo que conecta con el capítulo 27. Cuando una frontera dinámica separa dos resultados, una perturbación de redondeo o de integración puede desplazar una condición inicial a uno u otro lado. Con h=0,20h=0{,}20, las condiciones xu−10−6x_u-10^{-6} y xu+10−6x_u+10^{-6} terminan en resultados opuestos. Desde xu−10−6x_u-10^{-6} la trayectoria alcanza x=0x=0 cerca de t=27t=27, mientras que desde xu+10−6x_u+10^{-6} se acerca al equilibrio xs≈0,7236x_s\approx0{,}7236. En sistemas caóticos la cuestión será todavía más delicada y no se resolverá simplemente aumentando la precisión.

La resolución puede crear la apariencia de una transición o esconderla

Supongamos que se guarda solo el resultado de una malla gruesa y se dibuja una figura con líneas que unen los puntos. La interpolación gráfica puede sugerir una frontera continua donde el modelo nunca se ejecutó. Si además las regiones se colorean según una clasificación binaria, el ojo puede tomar el borde entre celdas como una estimación precisa del umbral. Esa precisión la añade la figura, porque el modelo no se evaluó entre los puntos de la malla.

También ocurre el problema inverso. Una malla demasiado gruesa puede saltarse por completo una ventana estrecha de comportamiento distinto. Que dos puntos vecinos de la cuadrícula tengan el mismo resultado no implica que todos los valores intermedios lo compartan.

Por eso una estrategia razonable de exploración es jerárquica. Comienza con un barrido amplio y relativamente grueso, que identifica regiones de variación rápida o de sensibilidad elevada. En esas regiones se refina la resolución y se comprueba que la conclusión se mantiene ante un refinamiento adicional. Cuando corresponde, el barrido se reemplaza por métodos diseñados para seguir la estructura matemática relevante.

Esa última etapa es importante. En problemas de bifurcación, la continuación numérica puede seguir ramas de soluciones y detectar puntos especiales de una manera que un barrido desde condiciones iniciales fijas no puede reproducir (Govaerts, 2000; Kuznetsov, 2004; Uecker, 2022). Un barrido paramétrico y una continuación numérica son herramientas relacionadas, pero no intercambiables.

Qué guardar durante un barrido

Almacenar la trayectoria completa de cada simulación puede producir volúmenes enormes de datos. Guardar demasiado tampoco garantiza un mejor análisis. Hay que distinguir entre los datos necesarios para la pregunta principal y los datos suficientes para auditar el cálculo. Para este ejemplo, el cuadro 28.1 propone los campos de una fila por simulación.

Cuadro 28.1. Campos de una fila por simulación en el barrido del capítulo.

CampoSignificado
hparámetro explorado
x0condición inicial
Thorizonte máximo
eventosi ocurrió el cruce de x = 0
t_eventotiempo del evento, si existe
x_finalúltimo estado calculado
statusmodo de terminación del integrador
nfevevaluaciones de la función dinámica
rtol, atoltolerancias numéricas

Si la pregunta requiere oscilaciones, máximos locales o estadísticas temporales, esos resultados deben añadirse explícitamente. Si más tarde se descubre que faltaba una cantidad, habrá que repetir la simulación o haber conservado una representación de la trayectoria.

La estructura de almacenamiento debe evitar una ambigüedad frecuente. Un NaN en t_evento no significa que el evento sea imposible. Solo indica que no se detectó en las condiciones de ejecución especificadas. El evento pudo ocurrir después de TT o pasar inadvertido para el procedimiento de detección. También es posible que la integración haya fallado antes.

Ejemplo reproducible de barrido y refinamiento

Una implementación compacta del barrido y de su refinamiento puede escribirse así.

import numpy as np
from scipy.integrate import solve_ivp


def rhs(t, y, h):
    x = y[0]
    return [x * (1 - x) - h]


def evento_cero(t, y, h):
    return y[0]


evento_cero.terminal = True
evento_cero.direction = -1


def simular(h, x0=1.0, T=80.0):
    sol = solve_ivp(
        rhs,
        (0, T),
        [x0],
        args=(h,),
        events=evento_cero,
        rtol=1e-9,
        atol=1e-11,
        max_step=0.2,
    )
    hubo_evento = len(sol.t_events[0]) > 0
    t_evento = sol.t_events[0][0] if hubo_evento else np.nan
    return {
        "h": h,
        "x0": x0,
        "T": T,
        "evento": hubo_evento,
        "t_evento": t_evento,
        "x_final": sol.y[0, -1],
        "status": sol.status,
        "nfev": sol.nfev,
    }


malla_gruesa = np.arange(0.10, 0.401, 0.02)
resultados = [simular(h) for h in malla_gruesa]

# Refinamiento alrededor del cambio observado entre 0,24 y 0,26
malla_fina = np.arange(0.235, 0.271, 0.001)
resultados_finos = [simular(h) for h in malla_fina]

sin_evento = [r["h"] for r in resultados_finos if not r["evento"]]
con_evento = [r["h"] for r in resultados_finos if r["evento"]]
print(f"{max(sin_evento):.3f} {min(con_evento):.3f}")

El código no pretende ser un sistema general de experimentación. Su propósito es hacer visibles las decisiones mínimas, que son la malla, la condición inicial, el horizonte, el evento, las tolerancias y las cantidades almacenadas. La última línea imprime 0.251 0.252, el cambio aparente que muestra la figura 28.2.

Un estudio real añadiría control de versiones de la configuración, semillas si hay aleatoriedad, manejo de fallos, exportación estructurada y pruebas automatizadas. El capítulo 29 estará dedicado a la verificación del código mediante soluciones conocidas e invariantes, entre otros recursos.

Transición observada y mecanismo

Nuestro ejemplo es deliberadamente favorable, porque conocemos las ecuaciones y podemos demostrar que los equilibrios colisionan en h=1/4h=1/4. En un sistema empírico, observar que una variable cambia bruscamente al modificar un parámetro no basta para afirmar que existe una bifurcación silla-nodo. Otras explicaciones posibles incluyen un umbral impuesto por la medición, una regla de decisión discontinua, una transición estocástica, una mezcla de poblaciones, un error numérico o una variable no modelada que cambia al mismo tiempo.

El barrido muestra qué hace el modelo especificado. Para usar ese resultado como explicación de un sistema real hay que justificar además la representación y contrastar predicciones discriminantes.

La distinción será especialmente importante al aplicar estas herramientas a la psiquiatría. Una trayectoria clínica puede mostrar una recaída abrupta y un modelo dinámico puede generar cambios abruptos mediante una bifurcación, pero ambas afirmaciones juntas no demuestran que esa recaída concreta se haya producido por una bifurcación. El modelo debe generar consecuencias adicionales que puedan compararse con observaciones.

Desde el punto de vista filosófico, un mapa de parámetros es una propiedad de una familia de modelos. Su relación con el fenómeno empírico pasa por la manera en que se definieron y estimaron los parámetros, aunque tengan nombres clínicos o biológicos.

Límites de la exploración exhaustiva

El número de simulaciones crece rápidamente con el número de dimensiones. Cien valores de un parámetro exigen cien ejecuciones. Con cien valores en cada uno de cinco parámetros, una cuadrícula completa tendría 1005=1010100^5=10^{10} combinaciones. Si cada simulación tardara un milisegundo, el barrido completo duraría unos 116 días.

El problema es una forma de la llamada maldición de la dimensionalidad. Una cuadrícula uniforme desperdicia cálculos en regiones poco informativas y puede seguir siendo demasiado pobre donde el comportamiento cambia rápidamente.

Más adelante estudiaremos métodos de muestreo, análisis de sensibilidad, diseños adaptativos y reducción dimensional. Por ahora basta una regla metodológica. La resolución debe elegirse según la estructura que se quiere detectar. Un número arbitrario de puntos por eje carece de justificación.

El mismo principio se aplica a las condiciones iniciales. En un sistema de diez variables, barrer cien valores por variable es inviable. Se requieren hipótesis sobre las regiones relevantes o procedimientos de muestreo más eficientes.

La limitación computacional tiene también una consecuencia metodológica, porque obliga a declarar qué parte del espacio se considera relevante. Esa decisión forma parte de la inferencia.

Consolidación

Un barrido paramétrico es un experimento computacional sobre una familia de configuraciones. Su calidad depende de que los parámetros, las condiciones iniciales, el horizonte, los eventos, la resolución y los resultados almacenados estén alineados con la pregunta. El número bruto de simulaciones importa mucho menos.

En el modelo x˙=x(1−x)−h\dot x=x(1-x)-h, el análisis exacto muestra un cambio estructural en h=1/4h=1/4. Una malla gruesa lo localiza solo dentro de un intervalo. Refinarla aumenta la resolución empírica del barrido, aunque no reemplaza el análisis matemático ni los métodos de continuación. Con un horizonte fijo, además, la malla fina puede desplazar la transición aparente, como ocurrió al situarla entre 0,2510{,}251 y 0,2520{,}252.

Los eventos permiten registrar condiciones relevantes en el momento en que ocurren. En este modelo, el tiempo hasta x=0x=0 puede calcularse exactamente cuando h>1/4h>1/4, lo que da una referencia independiente para verificar el detector numérico. Ese tiempo crece sin límite al acercarse al umbral, de modo que un horizonte finito puede modificar la clasificación observada.

La condición inicial también es parte del experimento. Para h<1/4h<1/4, el equilibrio inestable xu(h)x_u(h) separa resultados distintos. Barrer solo el parámetro habría ocultado esa dependencia. La resolución, finalmente, puede esconder estructura o fabricar su apariencia, porque la cuadrícula y la clasificación son decisiones de representación que deben documentarse y someterse a pruebas de refinamiento.

Ejercicios

Ejercicio 28.1. Con h=0,16h=0{,}16, (a) calcule los equilibrios xux_u y xsx_s. (b) Evalúe fxf_x en cada uno y clasifique su estabilidad. (c) Calcule la escala temporal τ\tau del equilibrio estable.

Ejercicio 28.2. (a) Calcule t0(h)t_0(h) para h=0,26h=0{,}26 y para h=0,30h=0{,}30. (b) Muestre que t0(h)=π/δ−4+O(δ)t_0(h)=\pi/\sqrt\delta-4+O(\delta) cuando δ→0\delta\to0. (c) Compare esa aproximación con el valor exacto para h=0,251h=0{,}251.

Ejercicio 28.3. (a) Determine numéricamente el menor hh para el cual el evento ocurre antes de T=80T=80. (b) Explique con ese valor por qué la malla fina del capítulo sitúa el cambio entre 0,2510{,}251 y 0,2520{,}252. (c) ¿Qué horizonte haría falta para detectar el evento con h=0,2501h=0{,}2501?

Ejercicio 28.4. Con h=0,2499h=0{,}2499 y x0=1x_0=1, la variable u=x−1/2u=x-1/2 cumple u˙=a2−u2\dot u=a^2-u^2 con a=0,01a=0{,}01. (a) Verifique que u(t)=acoth⁡(at+C)u(t)=a\coth(at+C) es solución y calcule CC. (b) Calcule x(20)x(20) y x(80)x(80). (c) ¿Cuánto tarda la trayectoria en quedar a menos de 0,0010{,}001 de xsx_s?

Ejercicio 28.5. Con h=0,20h=0{,}20, (a) calcule xux_u. (b) Una malla de condiciones iniciales con paso 0,050{,}05 comienza en 0,050{,}05. ¿Entre qué valores sitúa la frontera? (c) Explique por qué las condiciones xu±10−6x_u\pm10^{-6} terminan en resultados opuestos y por qué el redondeo podría decidir la clasificación de un punto todavía más cercano.

Ejercicio 28.6. (a) ¿Cuántas simulaciones exige una cuadrícula de 100 valores en cada uno de cinco parámetros y cuánto tardaría si cada una dura un milisegundo? (b) Con un presupuesto de 10510^5 simulaciones, ¿cuántos valores por eje admite una cuadrícula completa en cinco dimensiones? (c) ¿Qué alternativas a la cuadrícula completa menciona el capítulo?

Consultar las soluciones de los ejercicios

Soluciones de los ejercicios

Solución 28.1. (a) Como 1−0,64=0,6\sqrt{1-0{,}64}=0{,}6, los equilibrios son xu=0,2x_u=0{,}2 y xs=0,8x_s=0{,}8. (b) Se obtiene fx(0,2)=0,6>0f_x(0{,}2)=0{,}6>0, así que xux_u es inestable. En xsx_s se obtiene fx(0,8)=−0,6<0f_x(0{,}8)=-0{,}6<0, de modo que xsx_s es estable. (c) La escala es τ=1/0,6≈1,67\tau=1/0{,}6\approx1{,}67.

Solución 28.2. (a) Con h=0,26h=0{,}26, δ=0,01\delta=0{,}01 y t0=20arctan⁡5≈27,468t_0=20\arctan5\approx27{,}468. Con h=0,30h=0{,}30, t0≈10,288t_0\approx10{,}288. (b) Para v>0v>0 se cumple arctan⁡(1/v)=π/2−arctan⁡v\arctan(1/v)=\pi/2-\arctan v. Con v=2δv=2\sqrt\delta resulta t0=π/δ−(2/δ)arctan⁡(2δ)t_0=\pi/\sqrt\delta-(2/\sqrt\delta)\arctan(2\sqrt\delta). Como arctan⁡(2δ)=2δ+O(δ3/2)\arctan(2\sqrt\delta)=2\sqrt\delta+O(\delta^{3/2}), el segundo término vale 4+O(δ)4+O(\delta). (c) La aproximación da π/0,001−4≈95,346\pi/\sqrt{0{,}001}-4\approx95{,}346, frente al valor exacto 95,35195{,}351. La diferencia, de 0,0050{,}005, coincide con el término siguiente del desarrollo, 16δ/316\delta/3.

Solución 28.3. (a) La ecuación t0(h)=80t_0(h)=80 tiene solución h≈0,2514h\approx0{,}2514. (b) El valor 0,2510{,}251 queda por debajo de ese límite y su evento ocurre después de T=80T=80, mientras que 0,2520{,}252 queda por encima y su evento ocurre antes. (c) Con h=0,2501h=0{,}2501, δ=10−4\delta=10^{-4} y t0=200arctan⁡50≈310,2t_0=200\arctan50\approx310{,}2. Haría falta un horizonte mayor que 310.

Solución 28.4. (a) La derivada de acoth⁡(at+C)a\coth(at+C) es −a2/sinh⁡2(at+C)=a2(1−coth⁡2(at+C))=a2−u2-a^2/\sinh^2(at+C)=a^2\bigl(1-\coth^2(at+C)\bigr)=a^2-u^2. La condición u(0)=0,5u(0)=0{,}5 exige coth⁡C=50\coth C=50, es decir, C=12ln⁡(51/49)≈0,0200C=\tfrac12\ln(51/49)\approx0{,}0200. (b) Se obtiene x(20)≈0,5462x(20)\approx0{,}5462 y x(80)≈0,5148x(80)\approx0{,}5148. (c) La distancia a xsx_s es a(coth⁡(at+C)−1)a\bigl(\coth(at+C)-1\bigr) y queda bajo 0,0010{,}001 cuando coth⁡(at+C)<1,1\coth(at+C)<1{,}1. Eso ocurre si at+C>12ln⁡21≈1,522at+C>\tfrac12\ln21\approx1{,}522, es decir, para t>150,2t>150{,}2.

Solución 28.5. (a) Se obtiene xu=(1−0,2)/2≈0,2764x_u=(1-\sqrt{0{,}2})/2\approx0{,}2764. (b) La malla sitúa la frontera entre 0,250{,}25 y 0,300{,}30. (c) En xux_u la derivada del campo vale 0,2≈0,45>0\sqrt{0{,}2}\approx0{,}45>0, de modo que cualquier desviación inicial crece y su signo decide el resultado. Para un punto a una distancia de xux_u comparable con el error de redondeo del valor inicial, ese signo puede quedar decidido por el redondeo o por el error de integración.

Solución 28.6. (a) Exige 101010^{10} simulaciones, que a un milisegundo cada una tardarían unos 116 días. (b) Admite 10 valores por eje, porque 10510^5 combinaciones corresponden a 10 valores en cada uno de cinco ejes. (c) El capítulo menciona el refinamiento jerárquico y la continuación numérica, además de los métodos de muestreo, el análisis de sensibilidad y los diseños adaptativos. Todos exigen declarar qué región del espacio se considera relevante.

Referencias

Govaerts, W. J. F. (2000). Numerical methods for bifurcations of dynamical equilibria. Society for Industrial and Applied Mathematics.

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer.

Kuznetsov, Y. A. (2004). Elements of applied bifurcation theory (3.ª ed.). Springer.

SciPy Developers. (s. f.). scipy.integrate.solve_ivp. Recuperado el 19 de septiembre de 2026, de https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Strogatz, S. H. (2015). Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering (2.ª ed.). Westview Press.

Uecker, H. (2022). Continuation and bifurcation in nonlinear PDEs: Algorithms, applications, and experiments. Jahresbericht der Deutschen Mathematiker-Vereinigung, 124(1), 43–80. https://doi.org/10.1365/s13291-021-00241-5