Modelos Matemáticos do Caos

1. Introdução


Habr já discutiu a teoria do caos em artigos [1,2,3]. Os seguintes aspectos da teoria do caos são considerados nesses artigos: diagrama generalizado do gerador Chua; modelar a dinâmica do sistema Lorentz; atratores programados pelos circuitos lógicos Lorentz, Ressler, Rikitake e Nose-Hoover.

No entanto, técnicas da teoria do caos também são usadas para modelar sistemas biológicos, que, sem dúvida, são um dos sistemas mais caóticos de tudo o que se pode imaginar. Sistemas dinâmicos de equalização foram usados ​​para modelar tudo, desde crescimento populacional e epidêmico até batimentos cardíacos arrítmicos [4].

De fato, quase todo sistema caótico pode ser modelado - o mercado de valores mobiliários gera curvas que podem ser facilmente analisadas usando atratores estranhos, o processo de soltar quedas de uma torneira vazando é aleatório quando analisado com a orelha nua, mas se retratado como um atrator estranho, ele abre uma ordem sobrenatural que não se poderia esperar dos meios tradicionais.

O objetivo deste artigo é considerar a teoria do caos usando o exemplo de aumentar o número de populações biológicas e dobrar o ciclo em sistemas mecânicos com visualização gráfica de modelos matemáticos com base em programas simples e intuitivos escritos em Python.

O artigo foi escrito com o objetivo de ensinar, mas permitirá que mesmo um leitor sem experiência em programação use os programas acima para resolver independentemente a maioria dos novos problemas educacionais sobre o tópico modelagem de fenômenos do caos.

Dobrar o período de ciclos e caos no exemplo do crescimento do número de populações biológicas


Vamos começar examinando a equação diferencial logística que modela um aumento limitado e não exponencial no número de populações biológicas:

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

É essa equação que pode prever padrões de comportamento exóticos e inesperados de algumas populações. De fato, de acordo com (1), para t rightarrow+ infty o tamanho da população se aproxima do limite igual a / b .

Para a solução numérica da solução diferencial logística, você pode usar o algoritmo mais simples, com o valor numérico da etapa do tempo, considerando tn+1=tn+h , a solução (1) pode ser obtida aplicando repetidamente a seguinte relação:

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

Representamos a equação (2) na forma de uma equação logística em diferenças finitas:

Nn+1=rNnsNn2 . (3)

onde: r = 1 + ah es = bh .
Substituição em (3) Nn= fracrsxn obtemos a fórmula iterativa:

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

Calculando os valores dados pela relação (3), podemos gerar uma sequência x1,x2,x3,.....
os valores máximos do número de populações que o ambiente suportará em determinados momentos t1,t2,t3. .

Assumimos que existe um valor limitador de frações que expressa uma parte do tamanho da população:

x infty= limn to inftyxn (5)

Vamos investigar como isso depende x infty do parâmetro de crescimento r na equação (4). Para fazer isso, em Python, escrevemos um programa que, começando com x1=$0, calcula os resultados em várias centenas de iterações ( n = 200 ) para r = 1,5; 2,0; 2,5 :

Código do 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])) 


O resultado do programa (para reduzir a saída dos resultados, são fornecidos os três primeiros e os quatro últimos 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 

Uma análise do modelo discreto mostra que para r = 1,5; 2,0; 2,5 com número crescente de iterações, o valor xn estabiliza e se torna quase igual ao limite x infty , que é determinado pela relação (5). Além disso, para os valores dados de r, a quantidade x infty correspondentemente igual x infty=0,3333;0,5000;$0,600 .
Aumentamos r = 3,1; 3,25; 3,5 e o número de iterações n = 1008; para isso, fazemos as seguintes alterações no programa:

Código do 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])) 


O resultado do programa (para reduzir a saída dos resultados, são fornecidos os quatro primeiros e os últimos oito 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 segue a partir dos dados acima, em vez de estabilizar perto de uma única população limitante, a parte fracionária da população flutua entre duas frações à medida que o tempo muda. Comparado com r = 3,1 , o período do ciclo para r = 3,25 duplica e para r = 3,5 quatro vezes.

Programa para exibição gráfica de ciclos de crescimento populacional
 # -*- 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() 


O resultado do programa:



Devido à duplicação do período de iteração, xn+1=rxn(1xn) tornou-se amplamente conhecido. Quando o valor da taxa de crescimento excede r = 3,56 , a duplicação do período acelera e o caos extremo já surge no ponto r = 3,57 . Para exibir o início do caos, usamos o seguinte programa:

Código do 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() 


O resultado do 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 

imagem

Escreveremos um programa para visualizar a dependência do comportamento da iteração no parâmetro de crescimento r . Para cada valor de r no intervalo a leqslantr leqslantb 1000 iterações são executadas para obter estabilidade. Em seguida, todos os 250 valores obtidos como resultado de iterações são plotados ao longo do eixo vertical, formando pontos ( r, x ):

Código do 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() 


O resultado na forma de um diagrama:



O gráfico resultante é chamado de "diagrama de ramificação" , que permite determinar se um determinado valor de r corresponde a um ciclo ou caos. O único valor do tamanho da população é determinado pelo valor r aproximadamente3 depois um ciclo com um período de 2 a r aproximadamente3,4 , depois um ciclo com um período de 4, depois um ciclo com um período de 8 em diante, com uma abordagem rápida ao caos.

Deve-se notar que as áreas verticais do espaço não preenchido no gráfico são as regiões r = 3,6 er = 3,7 , entre r = 3,7 er = 3,8 , entre r = 3,8 er = 3,9 onde a ordem cíclica retorna do caos anterior.
Considerar a aparência de um ciclo com um período múltiplo de 3 na região 3,8 leqslantr leqslant$3, faça alterações no programa anterior:

Código do 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() 


O resultado do programa:



O ciclo do período 3 aparece nas proximidades do ponto r = 3,83 e, em seguida, é dividido seqüencialmente em ciclos 6,12,24. A existência de um ciclo com o período 3 implica a presença de ciclos de qualquer outro período finito, bem como ciclos caóticos sem nenhum período.

O diagrama de ramificação permite acompanhar o desenvolvimento do sistema com uma mudança suave no parâmetro. Com um valor fixo do parâmetro, o diagrama de aranha (diagrama de Lamera) permite rastrear as órbitas dos pontos.

A construção de um diagrama de aranha permite identificar vários efeitos que são invisíveis no diagrama de ramificação. Vamos escrever um programa:

Código do 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 Lamera:

Diagrama de Lameria

Período duplicado em sistemas mecânicos


Considere uma equação diferencial que modela as vibrações amortecidas livres de um ponto material de uma determinada massa em uma mola não linear, na qual o amortecimento é determinado pela velocidade.

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

Na equação (6), o termo kx representa a força de uma mola linear aplicada a um ponto material de uma dada massa e o termo  betax3 representa a não linearidade real da primavera.

Se uma força atua no sistema de oscilações livres (6), o deslocamento do ponto do material da massa à qual essa força é aplicada é descrito pela equação diferencial de Duffing para oscilações forçadas:

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

A equação (7) para a maioria dos parâmetros incluídos é resolvida numericamente. O sistema mecânico para o modelo matemático de acordo com a equação (7) é mostrado na figura:

imagem

Uma característica do sistema fornecido é que, em vez de uma mola, é usada uma rosca de metal flexível, que oscila em um plano vertical, para o qual a constante de gancho k é negativa. Nesse esquema, os pontos de equilíbrio estável (a) e (c) e o ponto de equilíbrio instável (b).

Quando um ponto material é deslocado da posição (b), a força que atua sobre ele é repulsiva. Se a força periódica, por exemplo, criada pelo campo magnético oscilante for parcialmente amortecida pela resistência do ar. Então, a equação (7) é um modelo matemático aceitável para o deslocamento horizontal x (t) de um ponto de material com os seguintes intervalos de parâmetros k<0,c>0, beta>0 .

Para estudar o comportamento de um sistema não-linear, tomamos k=1,m=c= beta= ômega=1 então a equação diferencial (7) assume a forma:

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

Escrevemos um programa para a integração numérica da equação (8) nas condições iniciais x(0)=1,x(0)=0 no campo 100 leqslantt leqslant200 e para cada um dos seguintes valores de amplitude F0=0,6;0,7;0,75;0,8 e, em cada caso, plote as soluções para os planos x(t),x(t) e t,x(t) :

Código do 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 do programa

imagem



Essa transição do período de duplicação para o caos mostra o comportamento geral de um sistema mecânico não linear em resposta a uma alteração no parâmetro físico correspondente, por exemplo: k,m,c, beta, ômega,F0 . Tais fenômenos não ocorrem em sistemas mecânicos lineares.

Atrator Lorenz


A substituição na equação de Duffing por oscilações forçadas (7) leva a um sistema bidimensional não linear de equações diferenciais, que foi mostrado na listagem anterior. Um sistema tridimensional não linear de equações diferenciais aplicado a problemas meteorológicos foi considerado por E.N. Lorenz:

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

A solução do sistema (9) é melhor visualizada na projeção em um dos três planos. Escrevemos um programa de integração numérica para os valores dos parâmetros b = \ frac {8} {3}, s = 10, r = 28 e as condições iniciais x (0) = - 8, y (0) = 8, z (0) = 27 :

Código do 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() 


O resultado do programa:

imagem

Considerando a imagem no gráfico ao longo do tempo, pode-se supor que o ponto P (x (t), y (t), z (t)) faça um número aleatório de oscilações, direita ou esquerda. Para uma aplicação meteorológica do sistema Lorenz, após um número aleatório de dias claros, segue-se um número aleatório de dias chuvosos.

Considere um programa para mapear o atrator Lorentz no plano xyz para condições iniciais ligeiramente diferentes:

Código do 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() 


Os resultados do programa são mostrados nos seguintes gráficos:

imagem

imagem

A partir dos gráficos acima, conclui-se que uma mudança na condição inicial de 1,0 para 1.0001 altera drasticamente a natureza da mudança no atrator de Lorentz.

Sistema Rossler


Este é um sistema tridimensional não linear muito estudado:
 fracdxdt=yz,
 fracdydt=x alphay, (10)
 fracdzdt=b+x(xc).

Escreveremos um programa para a integração numérica do sistema (10) para os seguintes parâmetros a = 0,39, b = 2, c = 4 nas condições iniciais x (0) = 0, y (0) = 0, z (0) = 0 :

Código do 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() 


O resultado do programa:

Graph

No avião, a fita de Rossler parece um laço, mas no espaço acaba sendo torcida como uma fita de Moebius.

Conclusões


Para demonstrar o fenômeno do caos, são apresentados programas simples e intuitivos em uma linguagem de programação Python de alto nível, fáceis de atualizar para novos projetos sobre esse tópico. O artigo tem foco educacional e metodológico e pode ser usado no processo de aprendizagem.

Referências


  1. Um pouco sobre o caos e como criá-lo
  2. Um olhar crítico sobre o atrator Lorenz
  3. Geradores de caos FPGA
  4. Equações diferenciais e problemas de valor-limite: modelagem e cálculo usando o Mathematica, Maple e MATLAB. 3ª edição.: Per. do inglês - M .: LLC “I.D. Williams, 2008. - 1104 p .: III. Paral. tit. Inglês

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


All Articles