La ecuación logística es uno de los sistemas dinámicos más conocidos y estudiados de las matemáticas. Y es que aunque un sencillo programa en BASIC para un MSX de 8 bits basta para calcular su evolución, su comportamiento es tan complejo que literalmente se han escrito libros enteros tratando de desentrañarlo.
La ecuación logística se desarrolló inicialmente para describir la dinámica del crecimiento de poblaciones. Imaginemos por ejemplo que deseamos analizar el tamaño de una determinada colonia de conejos. En una primera aproximación, podemos considerar que el crecimiento de la colonia será proporcional a su tamaño, pues cuanto más grande sea la colonia, más parejas reproductoras de conejos se formarán y más crías habrá en la siguiente generación. Esta idea da lugar a la siguiente ecuación:
Pi+1 = k·Pi
Donde Pi es el tamaño de la colonia en la generación i, Pi+1 el tamaño en la siguiente generación i+1 y k un parámetro (en general, mayor que 1) que describe la rapidez del crecimiento de la colonia. Como es fácil de comprobar, este modelo resulta en un crecimiento indefinido (de hecho, exponencial) de la población hasta el infinito, por lo que si bien puede ser útil para describir la evolución de la colonia en una primera etapa de su existencia, no sirve para analizar su comportamiento a largo plazo, porque los ecosistemas no pueden soportar poblaciones de tamaño indefinidamente grande.
De ahí surge la modificación de este modelo que se conoce como ecuación logística, que es el objeto de esta entrada:
Si+1 = μ · Si · (1-Si)
En esta ecuación, Si se define como el ratio de la población respecto de una población máxima que hipotéticamente podría soportar su ecosistema, Si = (Población)/(Población máxima). De este modo, si la población alcanza su valor máximo, Si = 1 y el término (1-Si) de la ecuación hace que la población se extinga por completo en la siguiente generación; y, sin llegar a estos extremos, a medida que la población crece y Si se aproxima a 1, el término (1-Si) hace que el crecimiento de la población se ralentice.
El programa en MSX BASIC que se muestra en el vídeo presenta los cálculos que se obtienen con esta ecuación. Como puede verse, basta con un programa muy sencillo de apenas quince líneas (de las que de hecho la mayoría se dedican a representar los resultados) en un antiguo ordenador de 8 bits para hacer este cálculo de una forma bastante rápida y cómoda:
El programa calcula la evolución de la población para un tamaño inicial fijo S1 = 0.1 y para diferentes valores del parámetro μ que van desde 0 hasta 4; cada vez que la pantalla se borra, se comienza un nuevo cálculo con un valor del parámetro μ un poco más alto. A lo largo del eje horizontal se representa el avance en las generaciones, con cada punto correspondiendo al tamaño de una generación dada. Ese tamaño se representa en el eje vertical. El sonido que produce el programa también es proporcional al tamaño de la población: cuanto más agudo es, mayor es el tamaño.
La evolución de la dinámica de la ecuación a valores bajos de μ es bastante sencilla: con μ muy bajos, la población cae rápidamente a S = 0; es decir, para esos valores de μ, la población inicial de S1 = 0.1 no basta para sostener la colonia y esta se extingue. En cambio, con valores algo más altos de μ, la población acaba estabilizándose en un determinado tamaño mayor que cero.
Este comportamiento no es difícil de justificar y de calcular analíticamente. El tamaño en el que se estabiliza la población es lo que se denomina un punto fijo de la ecuación: para que se obtenga esa población estable, la solución proporcionada por la ecuación en dos generaciones sucesivas ha de ser constante; es decir:
Si+1 = μ · Si · (1-Si) = Si
Ecuación que se puede resolver fácilmente para la variable Si, obteniéndose dos soluciones: o bien S = 0 (una población que tiene inicialmente tamaño 0 permanece indefinidamente en ese valor), o bien S = 1 – 1/μ. El valor hacia el que tiende a estabilizarse la población en ese segundo grupo de simulaciones es precisamente este segundo.
Hasta aquí llega el comportamiento sencillo de la ecuación, el que se ajusta a lo que se pretendía con ella: un modelo dinámico que predice el crecimiento de la colonia hasta un determinado tamaño estable, no infinito, o bien su extinción, en determinadas condiciones. Lo interesante aparece, como suele ocurrir, cuando el modelo deja de hacer lo que se esperaba de él.
Así, tras algunas simulaciones en las que se observa que el tamaño de la población ya no tiende hacia su valor estacionario de una forma directa, sino de una forma oscilatoria, para μ = 3.25 se tiene un nuevo cambio cualitativo: la ecuación ya no tiende a ningún valor constante, sino que oscila indefinidamente entre dos valores alternativos; se tiene lo que se denomina una órbita periódica, en este caso de periodo 2, con lo que se quiere expresar que se tiene una solución formada por dos valores, cada uno de los cuales se repite cada dos generaciones (en este caso concreto, estos valores que se repiten son, aproximadamente: 0.4953, 0.8124, 0.4953, 0.8124, 0.4953…). Este comportamiento, aunque puede resultar inesperado, de nuevo no resulta demasiado difícil de justificar analíticamente. Se tiene que en esta órbita periódica de periodo 2, la solución ha de ser igual a la que se tenía dos generaciones antes, es decir:
Si+2 = μ · Si+1 · (1-Si+1)=μ2 · Si · (1-Si) · (1 – μ ·Si · (1-Si)) = Si
Ecuación que con un breve ejercicio de aritmética y de paciencia se puede resolver, obteniéndose que además de los dos puntos fijos ya mencionados ( S = 0 y S = 1 – 1/μ), se tienen los siguientes puntos que conforman la mencionada órbita periódica de periodo 2:
Sin embargo, la cosa no acaba ahí. Para μ = 3.5 se tiene otra órbita periódica, pero ya no de periodo 2, sino 4: 0.5009, 0.8750, 0.3828, 0.8270 y vuelta a empezar con 0.5009. Esa solución se puede analizar con el mismo procedimiento propuesto para la anterior órbita periódica de orden 2, aunque la cosa se va complicando, porque la resolución de órbitas periódicas de periodos crecientes implica trabajar con polinomios de orden también creciente. Y las complicaciones no acaban ahí, porque para μ = 3.75 y μ = 4 se obtienen respuestas aún más complejas, en las que etapas de una curiosa regularidad, como el momento hacia la mitad de la penúltima simulación (la correspondiente a μ = 3.75) en el que la solución se estabiliza en un valor casi constante durante varias iteraciones, se alternan con otras en las que parece imposible distinguir ningún orden.
El comportamiento es tan complejo que entra en lo que se define, de un modo matemático riguroso, como comportamiento caótico. La ecuación logística jugó de hecho un papel histórico prominente en el desarrollo de la teoría del caos, y a día de hoy sigue siendo, como se ha dicho más arriba, el objeto de libros enteros (para una introducción más breve, pero bastante completa, se pueden leer los capítulos 10 y 11 del libro Ordinary Differential Equations and Dynamical Systems de Gerald Teschl).
La definición técnica de lo que constituye un sistema caótico que se presenta en el libro de Teschl pretende capturar diferentes aspectos, de los que uno de los fundamentales es la sensibilidad. En efecto, como se puede apreciar en el vídeo, la ecuación logística es extremadamente sensible a cambios en el valor del parámetro μ, en especial cuando el parámetro crece y se aproxima a 4. De hecho, por encima de μ = 3.25, cambios en el tercer o cuarto decimal de μ pueden llevar ya no solo a cambios cuantitativos significativos en la respuesta, sino también a cambios cualitativos, que alteran completamente su forma, como el paso de una órbita periódica de periodo 2 a una de periodo 4, o a soluciones de formas más complejas; si no se tiene un MSX a mano, esta afirmación se puede verificar fácilmente construyéndose, por ejemplo, una hoja de Excel. Pero esta sensibilidad se extiende a otros factores, como por ejemplo el valor inicial de S.
Para ilustrarlo, en el siguiente vídeo se compara la evolución para varios valores iniciales de S (variando desde 0.1 hasta 0.9; es decir, valores que cubren prácticamente todo el rango de posibles valores de S con sentido físico), tomando μ = 2. La representación de los datos es similar a la del vídeo anterior, con la excepción de que para facilitar la comparación entre las curvas obtenidas, estas se representan juntas, sin borrar la pantalla entre cada una de ellas, y en diferentes colores. Como puede apreciarse, la sección inicial de la curva naturalmente se modifica cuando se cambia la condición inicial, pues la curva arranca desde un punto distinto, pero en todos los casos se tiene el mismo comportamiento a largo plazo: la estabilización en torno al punto fijo S = 1 – 1/μ = 1/2.
La verdad es que el vídeo resulta bastante aburrido y no puede culparse a nadie si no quiere verlo entero: salvo unos pocos puntos al principio de las curvas, en todos los casos rápidamente se obtiene la misma línea horizontal. Compárese esta situación con la que se obtiene al tomar μ = 4 y hacer cambios en la condición inicial tan pequeños como una unidad en el tercer decimal, es decir, tomar por ejemplo S = 0.499 en vez de 0.498 o 0.500:
Como puede observarse, ya desde las primeras generaciones la evolución dinámica se vuelve radicalmente diferente. Este fenómeno de la sensibilidad respecto de lo parámetros y las condiciones iniciales es extremadamente importante en la práctica, porque, en el mundo real, las condiciones iniciales y los parámetros nunca se pueden medir experimentalmente de forma absolutamente exacta: siempre van acompañados de una incertidumbre, por pequeña que sea. Lo que caracteriza a los sistemas caóticos no es que sean azarosos; al contrario, son completamente deterministas: en el caso de la ecuación logística, toda su evolución pasada o futura está determinada por una pequeña expresión matemática que ocupa apenas una línea. Los sistemas caóticos son perfectamente deterministas, pero lo son de tal modo que amplifican las inevitables incertidumbres en los parámetros de entrada, hasta extremos que vuelven imposible predecir su evolución a largo plazo.
Esta sensibilidad es tan grande que, de hecho, esta entrada comenzó con una pequeña (gran) mentira: la evolución de la ecuación logística no se puede calcular con un MSX. El motivo es que la sensibilidad de esta ecuación a los parámetros es tan alta que los errores de redondeo en los cálculos provocados por la lógica de 8 bits que emplea el MSX para almacenar los números son suficientes para corromper los resultados. Pero, siendo justos con el MSX, este problema no es solo suyo: en el MSX, con su pequeño procesador Z80, el problema es especialmente grave (es decir, se empezará a manifestar tras un número de iteraciones relativamente pequeño), pero, tarde o temprano, se acabará llegando a la misma situación en cualquier otro ordenador digital, sea cual sea su potencia; es decir, incluso aunque no se tuviese en cuenta la necesaria incertidumbre empírica en los parámetros del modelo y se los supusiese dados con absoluta exactitud, se acabaría chocando contra el mismo muro: la dinámica a un plazo suficientemente largo de la ecuación logística (o de cualquier otro sistema caótico) no se puede ni se podrá calcular nunca con ningún ordenador digital habido o por haber. Y si esta limitación no le parece particularmente importante, piense que, por ejemplo, la dinámica del sistema solar también es caótica.
Esta entrada forma parte de la serie de entradas dedicada a la ecuación logística.