Modelos matemáticos del caos

Introduccion


Habr ya discutió la teoría del caos en los artículos [1, 2, 3]. Los siguientes aspectos de la teoría del caos se consideran en estos artículos: diagrama generalizado del generador Chua; modelando la dinámica del sistema Lorentz; atractores programados por circuitos integrados lógicos Lorentz, Ressler, Rikitake y Nose-Hoover.

Sin embargo, las técnicas de la teoría del caos también se utilizan para modelar sistemas biológicos, que, sin duda, son uno de los sistemas más caóticos de todos los que se pueden imaginar. Se utilizaron sistemas de ecualización dinámica para modelar todo, desde el crecimiento de la población y la epidemia hasta los latidos cardíacos arrítmicos [4].

De hecho, casi cualquier sistema caótico puede ser modelado: el mercado de valores genera curvas que pueden analizarse fácilmente utilizando atractores extraños, el proceso de soltar gotas de un grifo que gotea es aleatorio cuando se analiza con el oído desnudo, pero si se presenta como un atractor extraño, se abre Un orden sobrenatural que no podía esperarse de los medios tradicionales.

El propósito de este artículo es examinar la teoría del caos usando el ejemplo de aumentar el número de poblaciones biológicas y duplicar el ciclo en sistemas mecánicos con visualización gráfica de modelos matemáticos basados ​​en programas intuitivos simples escritos en Python.

El artículo fue escrito con el propósito de enseñar, pero permitirá que incluso un lector sin experiencia en programación use los programas anteriores para resolver de forma independiente la mayoría de los nuevos problemas educativos sobre el tema del modelado del fenómeno del caos.

Duplicar el período de ciclos y caos en el ejemplo del crecimiento del número de poblaciones biológicas.


Comencemos observando la ecuación diferencial logística que modela un aumento limitado en lugar de exponencial en el número de poblaciones biológicas:

 fracdNdt=aNbN2,(a,b>0).(1)

Es esta ecuación la que puede predecir patrones de comportamiento exóticos e inesperados de algunas poblaciones. De hecho, de acuerdo con (1), para t rightarrow+ inftyEl tamaño de la población se aproxima al límite igual a a / b .

Para la solución numérica de la solución diferencial logística, puede usar el algoritmo más simple, con el valor numérico del paso de tiempo, tomando tn+1=tn+h, entonces la solución (1) se puede obtener aplicando repetidamente la siguiente relación:

Nn+1=Nn+(aNnbNn2)h.(2)

Representamos la ecuación (2) en forma de una ecuación logística en diferencias finitas:

Nn+1=rNnsNn2. (3)

donde: r = 1 + ah y s = bh .
Sustitución en (3) Nn= fracrsxnobtenemos la fórmula iterativa:

xn+1=rxn(1xn), (4)

Calculando los valores dados por la relación (3), podemos generar una secuencia x1,x2,x3,.....
Los valores máximos de la cantidad de poblaciones que el medio ambiente soportará en determinados momentos. t1,t2,t3..

Suponemos que hay un valor límite de fracciones que expresan una parte del tamaño de la población:

x infty= limn to inftyxn, (5).

Investigaremos cómo depende x inftydel parámetro de crecimiento r en la ecuación (4). Para hacer esto, en Python, escribimos un programa que, comenzando con x1=$0.calcula los resultados en varios cientos de iteraciones ( n = 200 ) para r = 1.5; 2.0; 2.5 :

Código del programa
# -*- coding: utf8 -*- from numpy import * print(" nr=1,5 r=2,0 r=2,5 ") M=zeros([201,3]) a=[1.5,2.0,2.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,201,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,201,1): if 0<=i<=2: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) elif 2<i<=5: print(".") elif 197<=i<=200: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) 


El resultado del programa (para reducir la salida de los resultados, se dan los primeros tres y últimos cuatro valores):

  nr=1,5 r=2,0 r=2,5 0 0.5000 0.5000 0.5000 1 0.3750 0.5000 0.6250 2 0.3516 0.5000 0.5859 . . . 197 0.3333 0.5000 0.6000 198 0.3333 0.5000 0.6000 199 0.3333 0.5000 0.6000 200 0.3333 0.5000 0.6000 

Un análisis del modelo discreto muestra que para r = 1.5; 2.0; 2.5 con un aumento en el número de iteraciones, el valor xnse estabiliza y se vuelve casi igual al límite x infty, que está determinado por la relación (5). Además, para los valores dados de r, la cantidad x inftycorrespondientemente igual x infty=0.3333;0.5000;$0.600.
Aumentamos r = 3.1; 3.25; 3.5 y el número de iteraciones n = 1008, para esto hacemos los siguientes cambios en el programa:

Código del programa
 # -*- coding: utf8 -*- from numpy import * print(" nr=3,1 r=3,25 r=3,5 ") M=zeros([1008,3]) a= [3.1,3.25,3.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,1008,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,1008,1): if 0<=i<=3: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) elif 4<i<=7: print(".") elif 1000<=i<=1007: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) 


El resultado del programa (para reducir la salida de los resultados, se dan los primeros cuatro y últimos ocho valores):

  nr=3,1 r=3,25 r=3,5 0 0.5000 0.5000 0.5000 1 0.7750 0.8125 0.8750 2 0.5406 0.4951 0.3828 3 0.7699 0.8124 0.8269 . . . 1000 0.5580 0.4953 0.5009 1001 0.7646 0.8124 0.8750 1002 0.5580 0.4953 0.3828 1003 0.7646 0.8124 0.8269 1004 0.5580 0.4953 0.5009 1005 0.7646 0.8124 0.8750 1006 0.5580 0.4953 0.3828 1007 0.7646 0.8124 0.8269 

Como se desprende de los datos anteriores, en lugar de estabilizarse cerca de una sola población limitante, la parte fraccional de la población fluctúa entre dos fracciones a medida que cambia el tiempo. En comparación con r = 3.1 , el período del ciclo para r = 3.25 se duplica y para r = 3.5 cuatro veces.

Programa para mostrar gráficamente los ciclos de crecimiento de la población.
 # -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * M=zeros([1008,3]) a= [3.1,3.25,3.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,1008,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) x=arange(987,999,1) y=M[987:999,0] y1=M[987:999,1] y2=M[987:999,2] plt.title('     r=3,1;3,25;3,5') plt.plot(x,y, label="T=1,ymax=%s,ymin=%s"%(round(max(y),3),round(min(y),3))) plt.plot(x,y1, label="T=2,ymax=%s,ymin=%s"%(round(max(y1),3),round(min(y1),3))) plt.plot(x,y2, label="T=4,ymax=%s,ymin=%s"%(round(max(y2),3),round(min(y2),3))) plt.grid() plt.legend(loc="best") plt.ylabel("x(n)") plt.xlabel("n") plt.show() 


El resultado del programa:



Debido a duplicar el período de iteración, xn+1=rxn(1xn)se hizo ampliamente conocido. Cuando el valor de la tasa de crecimiento excede r = 3.56 , la duplicación del período se acelera y el caos extremo ya surge en el punto r = 3.57 . Para mostrar el inicio del caos, utilizamos el siguiente programa:

Código del programa
 # -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * print(" nr=3,57 ") M=zeros([1041,1]) a= [3.57] for j in arange(0,1,1): M[0,j]=0.5 for j in arange(0,1,1): for i in arange(1,1041,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,1041,1): if 1000<=i<=1015: print(" {0: 7.0f} {1: 7.4f}" . format(i,M[i,0])) x=arange(999,1040,1) y=M[999:1040,0] plt.title('     r=3,57') plt.plot(x,y) plt.grid() plt.ylabel("x(n)") plt.xlabel("n") plt.show() 


El resultado del programa:
  nr=3,57 1000 0.4751 1001 0.8903 1002 0.3487 1003 0.8108 1004 0.5477 1005 0.8844 1006 0.3650 1007 0.8275 1008 0.5096 1009 0.8922 1010 0.3434 1011 0.8050 1012 0.5604 1013 0.8795 1014 0.3784 1015 0.8397 

imagen

Escribiremos un programa para visualizar la dependencia del comportamiento de iteración en el parámetro de crecimiento r . Para cada valor de r en el intervalo a leqslantr leqslantbSe realizan 1000 iteraciones para lograr la estabilidad. Luego, cada 250 valores obtenidos como resultado de iteraciones se trazan a lo largo del eje vertical, formando puntos ( r, x ):

Código del programa
 # -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import* N=1000 y=[] y.append(0.5) for r in arange(2.8,4.0,0.0001): for n in arange(1,N,1): y.append(round(r*y[n-1]*(1-y[n-1]),4)) y=y[N-250:N] x=[r ]*250 plt.plot( x,y, color='black', linestyle=' ', marker='.', markersize=1) plt.title("   2,8<= r<=4,0,0<=x<=1") plt.xlabel("r") plt.ylabel("y") plt.axvline(x=3.01,color='black',linestyle='--') plt.axvline(x=3.45, color='black',linestyle='--') plt.axvline(x=3.6,color='black',linestyle='--') plt.axvline(x=3.7,color='black',linestyle='--') plt.axvline(x=3.8,color='black',linestyle='--') plt.axvline(x=3.9,color='black',linestyle='--') plt.show() 


El resultado en forma de diagrama:



El gráfico resultante se denomina "diagrama de ramificación" , que le permite determinar si un valor r dado corresponde a un ciclo o caos. El único valor del tamaño de la población se determina al valor r aproximadamente3luego un ciclo con un período de 2 a r aproximadamente3.4, luego un ciclo con un período de 4, luego un ciclo con un período de 8 en adelante con un acercamiento rápido al caos.

Cabe señalar que las áreas verticales de espacio sin llenar en el gráfico son las regiones r = 3.6 y r = 3.7 , entre r = 3.7 y r = 3.8 , entre r = 3.8 y r = 3.9 donde el orden cíclico regresa del caos anterior.
Considerar la aparición de un ciclo con un período múltiplo de 3 en la región. 3.8 leqslantr leqslant$3.hacer cambios al programa anterior:

Código del programa
 # -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import* N=1000 y=[] y.append(0.5) for r in arange(3.8,3.9,0.0001): for n in arange(1,N,1): y.append(round(r*y[n-1]*(1-y[n-1]),4)) y=y[N-250:N] x=[r ]*250 plt.plot( x,y, color='black', linestyle=' ', marker='.', markersize=1) plt.title("   3,8<= r<=3,9,0<=x<=1") plt.xlabel("r") plt.ylabel("y") plt.axvline(x=3.83,color='black',linestyle='--') plt.show() 


El resultado del programa:



El ciclo del período 3 aparece cerca del punto r = 3.83 , y luego se divide secuencialmente en los ciclos 6,12,24. La existencia de un ciclo con período 3 implica la presencia de ciclos de cualquier otro período finito, así como ciclos caóticos sin ningún período en absoluto.

El diagrama de ramificación le permite seguir el desarrollo del sistema con un cambio suave en el parámetro. Con un valor fijo del parámetro, el diagrama de araña (diagrama de Lamera) permite trazar las órbitas de los puntos.

La construcción de un diagrama de araña le permite identificar varios efectos que son invisibles en el diagrama de rama. Escribamos un programa:

Código del programa
 # -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * a=2.7 x1=0.62 def ff(x): return a*x*(1-x) b=a*x1*(1-x1)/x1 def fl(x): return b*x x=0.1 y=0 Y=[] X=[] for i in arange(1,30,1): X.append(x) Y.append(y) y=ff(x) X.append(x) Y.append(y) x=y/b plt.title('   \n  λx(1-x)  λ = 2.7') plt.plot(X,Y,'r') x1=arange(0,1,0.001) y1=[ff(x) for x in x1] y2=[fl(x) for x in x1] plt.plot(x1,y1,'b') plt.plot(x1,y2,'g') plt.grid(True) plt.show() 


Diagrama de Lamera:

Diagrama de lameria

Período de duplicación en sistemas mecánicos


Considere una ecuación diferencial que modela las vibraciones amortiguadas libres de un punto material de una masa dada en un resorte no lineal, en el que la amortiguación está determinada por la velocidad.

mx+cx+kx+ betax3=0(6)

En la ecuación (6), el término kx representa la fuerza de un resorte lineal aplicado a un punto material de una masa dada, y el término  betax3representa la no linealidad real de la primavera.

Si una fuerza actúa sobre el sistema de oscilaciones libres (6), la ecuación diferencial de Duffing para oscilaciones forzadas describe el desplazamiento del punto material de la masa a la que se aplica esta fuerza:

mx+cx+kx+ betax3=F0cos omegat(7)

La ecuación (7) para la mayoría de los parámetros incluidos se resuelve numéricamente. El sistema mecánico para el modelo matemático de acuerdo con la ecuación (7) se muestra en la figura:

imagen

Una característica del sistema dado es que, en lugar de un resorte, se usa un hilo de metal flexible, que oscila en un plano vertical, para el cual la constante de gancho k es negativa. En este esquema, los puntos de equilibrio estable (a) y (c), y el punto de equilibrio inestable (b).

Cuando un punto material se desplaza de la posición (b), la fuerza que actúa sobre él es repulsiva. Si la fuerza periódica, por ejemplo, creada por el campo magnético oscilante está parcialmente amortiguada por la resistencia del aire. Entonces, la ecuación (7) es un modelo matemático aceptable para el desplazamiento horizontal x (t) de un punto material con los siguientes rangos de parámetros k<0,c>0, beta>0.

Para estudiar el comportamiento de dicho sistema no lineal, tomamos k=1,m=c= beta= omega=1entonces la ecuación diferencial (7) toma la forma:

x+xx+x3=F0cos(t), (8)

Escribimos un programa para la integración numérica de la ecuación (8) bajo las condiciones iniciales x(0)=1,x(0)=0en el campo 100 leqslantt leqslant200y para cada uno de los siguientes valores de amplitud F0=0.6;0.7;0.75;0.8, y en cada caso, grafica las soluciones para los planos x(t),x(t)y t,x(t):

Código del programa
 # -*- coding: utf8 -*- from numpy import * from scipy. integrate import odeint import matplotlib.pyplot as plt for F in [0.6,0.7,0.75,0.8]: def f(y,t): y1,y2=y return [y2,-y2-y1**3+y1+F*cos(t)] t=arange(100,200,0.001) y0=[1.0,0.0] [y1,y2]=odeint(f, y0, t, full_output=False,rtol=1e-12).T if F==0.6: plt.subplot(221) plt.title('  F=0.6,T=2'r'$\pi$') plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) plt.subplot(222) plt.title(' x(t): F=0.6,T=2'r'$\pi$') plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) elif F==0.7: plt.subplot(223) plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1, label='  \n F=0.7,T=4'r'$\pi$') plt.legend(loc='upper left') plt.grid(True) plt.subplot(224) plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1, label=' x(t): F=0.7,T=4'r'$\pi$') plt.legend(loc='upper left') plt.grid(True) plt.show() if F==0.75: plt.subplot(221) plt.title('  F=0.75,T=8'r'$\pi$') plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) plt.subplot(222) plt.title(' x(t): F=0.75,T=8'r'$\pi$') plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) elif F==0.8: plt.subplot(223) plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1, label=' \n F=0.8,') plt.legend(loc='upper left') plt.grid(True) plt.subplot(224) plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1, label=' x(t): F=0.8,') plt.legend(loc='upper left') plt.grid(True) plt.show() 


Gráficos como resultado del programa

imagen



Esta transición del período doble al caos muestra el comportamiento general de un sistema mecánico no lineal en respuesta a un cambio en el parámetro físico correspondiente, por ejemplo: k,m,c, beta, omega,F0. Tales fenómenos no ocurren en sistemas mecánicos lineales.

Atractor Lorenz


La sustitución en la ecuación de Duffing por oscilaciones forzadas (7) conduce a un sistema bidimensional no lineal de ecuaciones diferenciales, que se mostró en la lista anterior. E.N. consideró un sistema tridimensional no lineal de ecuaciones diferenciales aplicado a problemas meteorológicos. Lorenz:

 fracdxdt=sx+sy,
 fracdydt=xz+rxy,(9)
 fracdzdt=xydz

La solución al sistema (9) se ve mejor en proyección en uno de los tres planos. Escribimos un programa de integración numérica para los valores de los parámetros b = \ frac {8} {3}, s = 10, r = 28 y las condiciones iniciales x (0) = - 8, y (0) = 8, z (0) = 27 :

Código del programa
 # -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt s,r,b=10,28,8/3 def f(y, t): y1, y2, y3 = y return [s*(y2-y1), -y2+(r-y3)*y1, -b*y3+y1*y2] t = np.linspace(0,20,2001) y0 = [-8, 8, 27] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T plt.plot(y1,y3, color='black', linestyle=' ', marker='.', markersize=2) plt.xlabel('x') plt.ylabel('z') plt.grid(True) plt.title("     xz") plt.show() 


El resultado del programa:

imagen

Teniendo en cuenta la imagen en el gráfico a lo largo del tiempo, se puede suponer que el punto P (x (t), y (t), z (t)) genera un número aleatorio de oscilaciones, derecha o izquierda. Para una aplicación meteorológica del sistema Lorenz, después de un número aleatorio de días despejados, sigue un número aleatorio de días lluviosos.

Considere un programa para mapear el atractor de Lorentz en el plano xyz para condiciones iniciales ligeramente diferentes:

Código del programa
 # -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #     . s,r,b=10,25,3 def f(y, t): y1, y2, y3 = y return [s*(y2-y1), -y2+(r-y3)*y1, -b*y3+y1*y2] #        t = np.linspace(0,20,2001) y0 = [1, -1, 10] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white') ax=Axes3D(fig) ax.plot(y1,y2,y3,linewidth=2) plt.xlabel('y1') plt.ylabel('y2') plt.title(" : y0 = [1, -1, 10]") y0 = [1.0001, -1, 10] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white') ax=Axes3D(fig) ax.plot(y1,y2,y3,linewidth=2) plt.xlabel('y1') plt.ylabel('y2') plt.title(" : y0 = [1.0001, -1, 10]") plt.show() 


Los resultados del programa se muestran en los siguientes gráficos:

imagen

imagen

De los gráficos anteriores se deduce que un cambio en la condición inicial de 1.0 a 1.0001 cambia drásticamente la naturaleza del cambio en el atractor de Lorentz.

Sistema de Rossler


Este es un sistema tridimensional no lineal muy estudiado:
 fracdxdt=yz,
 fracdydt=x alphay,(10)
 fracdzdt=b+x(xc).

Escribiremos un programa para la integración numérica del sistema (10) para los siguientes parámetros a = 0.39, b = 2, c = 4 bajo las condiciones iniciales x (0) = 0, y (0) = 0, z (0) = 0 :

Código del programa
 # -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt a,b,c=0.398,2.0,4.0 def f(y, t): y1, y2, y3 = y return [(-y2-y3), y1+a*y2, b+y3*(y1-c)] t = np.linspace(0,50,5001) y0 = [0,0, 0] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=2) plt.xlabel('x') plt.ylabel('y') plt.grid(True) plt.title("     xy") plt.show() 


El resultado del programa:

Graph

En el avión, la cinta de Rossler parece un lazo, pero en el espacio parece estar torcida como una cinta de Moebius.

Conclusiones


Para demostrar los fenómenos del caos, se presentan programas simples e intuitivos en un lenguaje de programación Python de alto nivel, que son fáciles de actualizar a nuevos proyectos sobre este tema. El artículo tiene un enfoque educativo y metodológico y puede usarse en el proceso de aprendizaje.

Referencias


  1. Un poco sobre el caos y cómo crearlo
  2. Una mirada crítica al atractor de Lorenz
  3. FPGA Chaos Generators
  4. Ecuaciones diferenciales y problemas de valor límite: modelado y cálculo usando Mathematica, Maple y MATLAB. 3ª edición.: Per. del ingles - M .: LLC "I.D. Williams, 2008. - 1104 p .: Ill. - Paral. tit. Ingles

Source: https://habr.com/ru/post/es436014/


All Articles