Diagramas de bifurcación con incertidumbre numérica

Qué objeto construye el computador

Un diagrama de bifurcación muestra cómo cambia la organización de las soluciones al variar un parámetro, la noción que definió el capítulo 43 y que los capítulos 53 y 54 siguieron en el mapa logístico. El capítulo 58 añadió dos dificultades, porque las ventanas periódicas pueden ocupar intervalos muy estrechos y, en sistemas multiestables, una sola condición inicial puede revelar solo parte de los atractores. El diagrama parece reunir todo eso en una figura y, por su misma capacidad de condensación, puede ocultar decisiones numéricas importantes.

Un diagrama calculado por computador no es el conjunto abstracto de todas las soluciones de una familia. Es el resultado de un experimento finito. Se elige una malla de parámetros

un conjunto de $S$ condiciones iniciales, un número $B$ de iteraciones que se descartan como transitorio, un número $M$ de iteraciones que se conservan y una precisión aritmética $p$. Conviene dar un nombre explícito al resultado, $D_{J,B,M,S,p}$, porque no hay ninguna razón matemática para suponer que dos elecciones distintas produzcan figuras indistinguibles. La estabilidad del resultado frente a cambios razonables de esas elecciones es, en sí misma, una propiedad que hay que comprobar.

La incertidumbre de la que hablamos aquí no es aleatoriedad ni un intervalo de confianza estadístico. Surge de observar una estructura asintótica e infinitamente detallada mediante un cálculo finito. Algunas fuentes de error pueden acotarse analíticamente y otras se evalúan mediante estudios de convergencia del protocolo. Este capítulo usa el mapa logístico como banco de pruebas, porque conocemos varios de sus umbrales de forma exacta.

El transitorio no es una cantidad universal

La decisión más habitual consiste en descartar las primeras 500 o 1000 iteraciones y suponer que el resto representa el régimen asintótico. La regla funciona lejos de los umbrales y no existe un número universal de iteraciones de calentamiento. Si $x^*$ es un punto fijo hiperbólico con multiplicador $\lambda$, una desviación pequeña cumple $|e_n|\approx|\lambda|^n|e_0|$ mientras valga la aproximación lineal, de modo que reducir un error inicial $d_0$ hasta una tolerancia $\varepsilon$ exige

El denominador tiende a cero cuando $|\lambda|\to1$, por lo que el número de iteraciones necesario diverge al acercarse a una pérdida de estabilidad. Es la versión discreta del enlentecimiento crítico de los capítulos 53 y 57. En el mapa logístico, el punto fijo no nulo tiene multiplicador $2-r$. Reducir una desviación de $10^{-1}$ a $10^{-8}$ requiere unas 153 iteraciones con $r=2{,}9$, unas 1604 con $r=2{,}99$ y unas 16 111 con $r=2{,}999$. Un protocolo con $B=500$ sería holgado en el primer caso e insuficiente en el tercero.

El efecto es visible en la figura 59.1. Al barrer la región de la primera duplicación, situada exactamente en $r=3$, un transitorio de 100 iteraciones produce dos bandas para $r<3$ que podrían confundirse con un ciclo de período 2. Son oscilaciones transitorias alrededor de un único punto fijo estable cuyo multiplicador está muy cerca de $-1$. Con 20 000 iteraciones descartadas, esas bandas colapsan sobre la rama fija antes del umbral. La conclusión no es que haya que usar siempre 20 000 iteraciones, porque eso sustituiría una regla arbitraria por otra. La longitud apropiada depende del régimen y de la precisión que exija la pregunta.

Resolución paramétrica

Otra fuente de incertidumbre está en el eje horizontal. Un diagrama no explora un continuo de parámetros, sino $J$ valores. Por eso una estructura más estrecha que la separación de la malla puede desaparecer por completo. El capítulo 58 situó el nacimiento del ciclo de período 3 del mapa logístico en $r_{SN}=1+\sqrt8\approx3{,}828427$ y su primera duplicación interna en $r\approx3{,}841499$. Ese primer segmento estable ocupa un ancho de $0{,}0131$.

Con un paso $\Delta r=0{,}03$, la ventana es más estrecha que la separación de la malla y su detección depende del desplazamiento de la malla. Si los valores son $3{,}80$, $3{,}83$ y $3{,}86$, el segundo cae dentro del segmento. Si son $3{,}82$, $3{,}85$ y $3{,}88$, ninguno lo hace y la figura no muestra rastro de la ventana (figura 59.2). Para una malla con origen arbitrario, la probabilidad de que algún punto caiga dentro de una estructura de ancho $w$ es aproximadamente $w/\Delta r$ cuando $w<\Delta r$, es decir, un 44 % en este caso. La ausencia de una rama en la figura no demuestra que la rama no exista.

La regla de interpretación es que la resolución debe ser menor que la escala de las estructuras que se quieren detectar, aunque esa escala suele desconocerse antes del análisis. Conviene entonces trabajar de forma adaptativa, con un barrido inicial que localice regiones de cambio y refinamientos posteriores donde aparezcan transiciones o cambios rápidos de una medida auxiliar como el exponente de Lyapunov. Aumentar sin más el número de parámetros tampoco resuelve todo, porque con un transitorio insuficiente una malla muy fina solo dibuja con mayor resolución el mismo artefacto.

Cuántos puntos conservar

El número $M$ de estados conservados se interpreta según el régimen. Si el atractor es un ciclo de período $p$, bastaría conservar al menos $p$ estados, pero como $p$ suele desconocerse, un valor pequeño puede submuestrear un ciclo largo y hacerlo parecer otro objeto. Si el atractor es aperiódico, $M$ determina qué parte de su distribución empírica alcanzamos a representar.

Conviene además distinguir una nube de puntos de una densidad invariante. El diagrama dibuja visitas. En las regiones caóticas, las zonas más densas reflejan en parte la frecuencia de visita, aunque la intensidad visual también depende del tamaño del marcador, la transparencia, el número de iteraciones y la resolución de la imagen. Una banda oscura no es una estimación normalizada de probabilidad. Cuando la pregunta exige comparar distribuciones, conviene construir histogramas y estudiar su convergencia, como hicimos en el capítulo 55 con la densidad del arcoseno. Cuando la pregunta es topológica, el énfasis debe estar en la periodicidad y la estabilidad. Un ciclo de período alto puede confundirse con caos porque la figura contiene muchos niveles verticales. El diagnóstico obliga a comprobar que $|F_r^p(x_n)-x_n|$ sea pequeño después del transitorio y descartar los divisores propios de $p$.

Condiciones iniciales y continuación

Muchos diagramas del mapa logístico usan una sola condición inicial genérica por parámetro. La práctica funciona bien allí porque el atractor relevante tiene una cuenca grande dentro del intervalo. No debe convertirse, con todo, en regla general. Retomemos el mapa biestable

cuyos puntos fijos cumplen $\mu+x-x^3=0$. Los pliegues de la curva de equilibrios aparecen donde además $1-3x^2=0$, es decir, en $x=\pm1/\sqrt3$ y $|\mu_f|=2/(3\sqrt3)\approx0{,}38490$, los mismos valores del caso resuelto del capítulo 50. Con $|\mu|<|\mu_f|$ hay tres puntos fijos, las dos ramas externas atraen y la intermedia separa las cuencas.

Si para cada $\mu$ se usa siempre una semilla positiva, se verá sobre todo la rama positiva. Una semilla negativa revela la otra. Si, en cambio, se hace un barrido ascendente usando como condición inicial el último estado del parámetro anterior, la trayectoria sigue una rama hasta que desaparece en un pliegue. El barrido descendente sigue la rama opuesta. El resultado es una figura de histéresis que depende de la historia del protocolo (figura 59.3a). La figura no contiene un error. El error aparecería al interpretarla como el conjunto completo de atractores disponibles para cada $\mu$.

Un diagrama puede mezclar, por tanto, dos experimentos distintos. En un barrido reinicializado, cada parámetro se estudia desde una semilla independiente. En un barrido por continuación, el resultado anterior se transporta al parámetro siguiente e incorpora dependencia de trayectoria en el procedimiento. Ninguno es superior en general, porque cada uno responde a una pregunta distinta. Para detectar coexistencia hacen falta varias condiciones iniciales por parámetro.

Precisión finita

Los errores anteriores podrían reducirse aumentando $B$, $M$, $J$ o $S$. La aritmética finita introduce una dificultad distinta. Los números de punto flotante forman un conjunto finito, por lo que un programa que itera de manera determinista sobre un tipo numérico fijo está iterando una función sobre un conjunto finito de estados. Por el principio del palomar, toda órbita computacional suficientemente larga repite algún estado y desde entonces es exactamente periódica, sin que el mapa real tenga por qué tener allí una órbita periódica atractora (Higham, 2002).

El efecto aparece mucho antes de una repetición exacta. Con $r=4$, donde el exponente de Lyapunov vale $\ln2$ para las órbitas típicas, los errores pequeños se duplican en promedio en cada iteración. Partiendo del mismo decimal, $0{,}4$, escrito en simple y en doble precisión, la diferencia entre ambas órbitas supera $10^{-4}$ en la iteración 12, $10^{-2}$ en la 19 y $10^{-1}$ en la 22 (figura 59.3b). Después, las dos secuencias ya no deben leerse como aproximaciones punto a punto de una misma trayectoria real. Con media precisión, el efecto es todavía más brusco. La órbita de $0{,}4$ repite exactamente en la iteración 44 un estado ya visitado en la iteración 4, con lo que entra en un ciclo de longitud 40 que no pertenece al mapa matemático.

Eso no vuelve inútil toda estadística calculada a largo plazo. La pregunta correcta es qué propiedad se quiere aproximar y con qué justificación. Una órbita numérica puede perder la correspondencia punto a punto con la órbita exacta y aun así producir estimaciones robustas de ciertas propiedades. La teoría de sombreado estudia cuándo una pseudoórbita numérica permanece cerca de alguna órbita verdadera, que no es la que compartiría su estado inicial idealizado. Hammel y colaboradores (1987) mostraron, para órbitas concretas del mapa logístico, que una pseudoórbita calculada con errores de redondeo pequeños permanece cerca de una órbita verdadera durante un número muy grande de iteraciones. El resultado se obtuvo para esos casos y no puede darse por supuesto en cualquier mapa. Conviene separar, por tanto, la predicción de una trayectoria de la estimación de propiedades del atractor, una distinción que el capítulo 60 usará para construir un horizonte de predicción.

Auditar el protocolo

La expresión incertidumbre numérica puede inducir a dibujar barras de error verticales alrededor de cada punto, lo que sería conceptualmente incorrecto, porque los puntos de un diagrama no son estimadores independientes de una media observacional. La incertidumbre pertinente es estructural respecto del protocolo y se expresa mejor mediante comparaciones de convergencia. Aumentar $B$ y comprobar si desaparecen ramas transitorias, refinar $\Delta\mu$ y comprobar si aparecen estructuras antes omitidas, aumentar $M$ y evaluar si cambia la periodicidad inferida, repetir con varias semillas cuando la multiestabilidad sea plausible y repetir con mayor precisión cuando la conclusión dependa de detalles finos.

El bloque siguiente ejecuta la primera de esas comparaciones. Clasifica el período del atractor en 1101 valores de $r$ con dos longitudes de transitorio y cuenta en cuántos parámetros cambia la clasificación.

import numpy as np

def periodo(r, transitorio, maximo=16, tol=1e-9):
    """Período detectado del atractor tras descartar un transitorio."""
    x = 0.3
    for _ in range(transitorio):
        x = r * x * (1 - x)
    orbita = [x]
    for _ in range(maximo):
        x = r * x * (1 - x)
        orbita.append(x)
    for p in range(1, maximo + 1):
        if abs(orbita[p] - orbita[0]) < tol:
            return p
    return 0                                   # ningún período detectado hasta el máximo

malla = np.linspace(2.9, 4.0, 1101)
corto = np.array([periodo(r, 1_000) for r in malla])
largo = np.array([periodo(r, 50_000) for r in malla])
cambian = corto != largo
print(f"parámetros examinados: {len(malla)}")
print(f"clasificación distinta con 1000 y con 50 000 iteraciones descartadas: {cambian.sum()}")
print(f"fracción: {cambian.mean():.3%}")
print("ejemplos de parámetros que cambian:",
      ", ".join(f"{r:.4f} ({a} -> {b})" for r, a, b in
                zip(malla[cambian][:4], corto[cambian][:4], largo[cambian][:4])))

Al ejecutarlo se obtiene esta salida.

parámetros examinados: 1101
clasificación distinta con 1000 y con 50 000 iteraciones descartadas: 40
fracción: 3.633%
ejemplos de parámetros que cambian: 2.9830 (2 -> 1), 2.9840 (2 -> 1), 2.9850 (2 -> 1), 2.9860 (2 -> 1)

El 3,6 % de los parámetros cambia de clasificación. Los primeros ejemplos están justo antes de $r=3$, donde el transitorio corto convierte una oscilación amortiguada en un falso ciclo de período 2. Un resultado así no invalida el diagrama, pero acota qué afirmaciones resisten el cambio de protocolo.

El resultado de estas comparaciones puede resumirse como una envolvente de conclusiones, del tipo «la ventana de período 3 permanece localizada entre estos parámetros con tres resoluciones, dos longitudes de transitorio y dos semillas genéricas». Esa afirmación es más informativa que declarar que el diagrama es preciso porque contiene millones de puntos. Cuando se busca un valor crítico concreto, además, el gráfico no es el instrumento final y debe reemplazarse por una condición matemática. Para una duplicación se busca la órbita y el cruce de su multiplicador por $-1$. Para una nodo-silla se buscan a la vez el punto fijo de la iterada y el multiplicador $+1$. El diagrama localiza y la ecuación confirma (Parker y Chua, 1989).

La palabra convergencia exige aquí una precisión adicional, porque no hay una distancia universal entre dos diagramas. En una región periódica, dos cálculos se comparan por el período detectado, la posición de las ramas y sus multiplicadores. En una región caótica, exigir coincidencia punto a punto carece de sentido después del horizonte de predictibilidad. Resulta más apropiado comparar el rango visitado, el exponente de Lyapunov o los histogramas de ocupación. Si hay multiestabilidad, la comparación debe incluir qué atractores se recuperan y con qué fracción de semillas.

Una práctica útil consiste en trabajar con mallas anidadas. Si una primera exploración usa paso $\Delta\mu$, la siguiente puede usar $\Delta\mu/2$ conservando los puntos anteriores, lo que permite distinguir las estructuras nuevas producidas por mayor resolución de las diferencias causadas por muestrear posiciones distintas. Conviene también definir un criterio de detención antes de mirar la figura final, por ejemplo exigir que duplicar el transitorio no cambie la clasificación de más de una fracción prefijada de los parámetros y que refinar la malla no desplace un umbral estimado más que una tolerancia declarada. Estos criterios no convierten una simulación en demostración, pero hacen explícito cuándo se considera estable una conclusión numérica. La auditoría debe conservar, por último, los resultados negativos. Si un umbral no converge con la resolución utilizada o si aparecen atractores nuevos al ampliar el conjunto de semillas, la respuesta correcta es mantener la conclusión abierta.

Cuatro conclusiones con distinto soporte

Cuatro afirmaciones obtenidas al explorar estos dos mapas ilustran el punto. Una dice que el punto fijo del mapa logístico pierde estabilidad en $r=3$. El diagrama es ilustrativo y la conclusión está respaldada analíticamente, porque el multiplicador vale $2-r$ y cruza $-1$ allí, de modo que la incertidumbre gráfica no afecta al valor exacto.

Otra dice que no existe ninguna ventana periódica entre $r=3{,}82$ y $r=3{,}88$. Si proviene de una malla de paso $0{,}03$, carece de soporte, porque la figura 59.2 muestra que la malla puede saltarse el intervalo periódico. La limitación es de resolución paramétrica.

Una tercera afirma que, con $\mu=0$, el mapa biestable tiene un único atractor positivo. Si el experimento usó solo semillas positivas, tampoco está justificada, porque las semillas negativas convergen a la rama negativa. La limitación es de cobertura del espacio de estados.

La última sostiene que la órbita calculada con $r=4$ predice el estado exacto cien iteraciones después. La comparación entre precisiones la refuta, ya que la separación alcanza el tamaño del intervalo mucho antes. De ahí no se sigue que el exponente de Lyapunov o la densidad invariante sean incalculables, sino que la pregunta debe pasar de la trayectoria individual a la propiedad estadística (Alligood et al., 1996).

Qué queda establecido

Un diagrama de bifurcación comprime la dependencia paramétrica y el comportamiento asintótico. Su poder visual no elimina la necesidad de análisis. Cerca de las bifurcaciones los transitorios se alargan, las estructuras estrechas pueden quedar entre los puntos de la malla, la multiestabilidad obliga a explorar más de una condición inicial y la aritmética finita convierte el continuo en un conjunto finito de estados representables. La respuesta adecuada no es desconfiar de toda simulación, sino especificar qué se calculó y comprobar si la conclusión resiste al refinar el procedimiento. Cuando existe una condición analítica para la bifurcación, debe usarse. Cuando la propiedad es estadística, debe comprobarse la convergencia estadística. Cuando la dinámica es multiestable, debe explorarse el espacio de estados.

El capítulo 60 cerrará esta parte con la última cuestión. Partiendo de dos condiciones iniciales muy próximas, seguirá su separación y calculará hasta cuándo una predicción puede considerarse informativa. El diagrama habrá mostrado cómo cambia la dinámica al variar los parámetros. El caso final mostrará qué ocurre cuando los parámetros se mantienen fijos y la información sobre el estado inicial es inevitablemente finita.

Ejercicios

Ejercicio 59.1. (a) Calcule cuántas iteraciones hacen falta para reducir una desviación de $10^{-1}$ a $10^{-8}$ en el punto fijo del mapa logístico con $r=2{,}95$. (b) Repita con $r=3{,}05$, usando el multiplicador del ciclo de período 2. (c) ¿Qué ocurre exactamente en $r=3$?

Ejercicio 59.2. (a) Calcule el ancho del segmento de período 3 estable del mapa logístico. (b) ¿Qué paso de malla garantiza al menos dos valores dentro de él? (c) Estime la probabilidad de detectarlo con un paso de $0{,}02$ y origen arbitrario.

Ejercicio 59.3. Un diagrama del mapa logístico con $B=200$ muestra dos ramas en $r=2{,}99$. (a) ¿Qué comprobación decidiría si el atractor es un punto fijo o un ciclo de período 2? (b) ¿Qué valor de $B$ sería razonable? (c) ¿Cómo se informaría el resultado si ambas clasificaciones aparecen según el protocolo?

Ejercicio 59.4. En el mapa biestable con $h=0{,}2$, (a) calcule los pliegues. (b) Describa qué ocurre en un barrido ascendente desde $\mu=-0{,}6$ hasta $\mu=0{,}6$. (c) ¿Qué protocolo revelaría los dos atractores en $\mu=0$?

Ejercicio 59.5. (a) Explique por qué toda órbita calculada en punto flotante termina siendo periódica. (b) ¿Qué distingue esa periodicidad de la del mapa matemático? (c) ¿Qué precaución sugiere para los diagramas calculados en media precisión?

Ejercicio 59.6. Proponga un protocolo de auditoría para la afirmación «entre $r=3{,}6$ y $r=3{,}7$ el mapa logístico no tiene ventanas periódicas de período menor que 10». Indique qué variaciones del procedimiento haría y qué resultado obligaría a retirar la afirmación.

Referencias

Alligood, K. T., Sauer, T. D. y Yorke, J. A. (1996). Chaos: An introduction to dynamical systems. Springer. https://doi.org/10.1007/b97589

Hammel, S. M., Yorke, J. A. y Grebogi, C. (1987). Do numerical orbits of chaotic dynamical processes represent true orbits? Journal of Complexity, 3(2), 136–145. https://doi.org/10.1016/0885-064X(87)90024-0

Higham, N. J. (2002). Accuracy and stability of numerical algorithms (2.ª ed.). SIAM. https://doi.org/10.1137/1.9780898718027

Parker, T. S. y Chua, L. O. (1989). Practical numerical algorithms for chaotic systems. Springer. https://doi.org/10.1007/978-1-4612-3486-9