Русский
Паблик ВК
Телеграм Канал
Магазин
Обратная связь
18+
SE7EN.ws

Математические модели хаоса

Научно-популярное
Прокомментировать
179

Введение

На Habr уже обсуждалась теория хаоса в статьях [1,2,3]. В этих статьях рассмотрены следующие аспекты теории хаоса: обобщённая схема генератора Чуа; моделирование динамики системы Лоренца; программируемые логическими интегральными схемами аттракторы Лоренца, Ресслера, Рикитаке и Нозе-Гувера.

Однако, техники теории хаоса используются и для моделирования биологических систем, которые, бесспорно, являются одними из наиболее хаотических систем из всех, что можно себе представить. Системы динамических равенств использовались для моделирования всего — от роста популяций и эпидемий, до аритмических сердцебиений [4].

В действительности, почти любая хаотическая система может быть смоделирована — рынок ценных бумаг порождает кривые, которые можно легко анализировать при помощи странных аттракторов, процесс падения капель из протекающего водопроводного крана кажется случайным при анализе невооруженным ухом, но, если его изобразить как странный аттрактор, открывается сверхъестественный порядок, которого нельзя было бы ожидать от традиционных средств.

Целью настоящей статьи является рассмотрение теории хаоса на примере роста численности биологических популяций и удвоения цикла в механических системах с графической визуализацией математических моделей основанной на простых интуитивно понятных программах, написанных на Python.

Статья написана с целью обучения, но позволит, даже не имеющему опыта программирования читателю, используя приведенные программы, самостоятельно решить большинство новых учебных задач по теме моделирования явлений хаоса.

Удвоение периода циклов и хаос на примере роста численности биологических популяций

Начнём с рассмотрения логистического дифференциального уравнения, которое моделирует ограниченный, а не экспоненциальный рост численности биологических популяций:

Математические модели хаоса численность популяции приближается к граничной равной a/b.

Для численного решения логистического дифференциального решения можно использовать самый простой алгоритм, с численным значением шага времени, приняв $t_{n+1}=t_{n}+h$, тогда решение (1) можно получить путём многократного применения следующего соотношения:

$N_{n+1}=N_{n}+(aN_n-bN{_{n}}^{2})h.$ (2)

Представим уравнение (2) в виде логистического уравнения в конечных разностях:

$N_{n+1}=rN_{n}-sN{_{n}}^{2}$. (3),

где: r=1+ah и s=bh.
Подстановкой в (3) $N_{n}=frac{r}{s}x_{n}$, получим итерационную формулу:

$x_{n+1}=rx_{n}(1-x_{n})$, (4)

Рассчитывая значения, заданные соотношением (3), можно сгенерировать последовательность $x_{1},x_{2},x_{3},.....$
максимальных значений численности популяций, которые будет поддерживать среда в заданные моменты времени $t_{1},t_{2},t_{3}.$.

Предполагаем, что существует предельное значение дробей, выражающих часть численности популяций:

$x_{infty }=lim_{n to infty }x_{n}$, (5).

Будем исследовать, как зависит $x_{infty }$ от параметра роста r в уравнении (4). Для этого на Python напишем программу, которая, начиная с $x_{1}=0,5$, вычисляет результаты при нескольких сотен итераций (n=200) для r=1,5;2,0;2,5:

Код программы

# -*- coding: utf8 -*- from numpy import * print("      n     r=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

Результат работы программы (для сокращения вывода результатов приведены первые три и последние четыре значения):

       n   r=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

Анализ дискретной модели показывает, что для r=1,5;2,0;2,5 с ростом количества итераций значение $x_{n}$ стабилизируется и становится практически равным предельному $x_{infty}$, которое определяется соотношением (5). Причём для приведенных значений r величина $x_{infty}$ соответственно равна $x_{infty}=0,3333;0,5000;0,6000$.
Увеличим r=3,1;3,25;3,5 и число итераций n=1008, для этого внесём следующие изменения в программу:

Код программы

 # -*- coding: utf8 -*- from numpy import * print("      n     r=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

Результат работы программы (для сокращения вывода результатов приведены первые четыре и последние восемь значений):

    n   r=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

Как следует из приведенных данных, вместо того, чтобы стабилизироваться возле единственной предельной численности популяции, дробная часть численности популяции колеблется между двумя дробями по мере изменения времени. По сравнению с r=3,1, период цикла для r=3,25 увеличивается вдвое, а для r=3,5 в четыре раза.

Программа для графического отображения циклов роста популяции

 # -*- 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()

Результат выполнения программы:

Благодаря удвоению периода итерация, $x_{n+1}=rx_{n}(1-x_{n})$ стала широко известной. Когда значение скорости роста превосходит r=3,56, удвоение периода ускоряется и уже в точке r=3,57 возникает чрезвычайный хаос. Для отображения наступления хаоса воспользуемся следующей программой:

Код программы

# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * print("     n     r=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()

Результат выполнения программы:

   n      r=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

image

Напишем программу для визуализации зависимости поведения итераций от параметра роста r. Для каждого значения r в интервале $aleqslant rleqslant b$ выполняется 1000 итераций для достижения устойчивости. Затем, каждые 250 значений, полученных в результате итераций, наносятся на график по вертикальной оси, образуя точки (r,x):

Код программы

 # -*- 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()

Результат в виде диаграммы:

Полученный график называется “диаграммой ветвления”, которая позволяет определить, чему соответствует данное значение r — циклу или хаосу. Единственное значение численности популяции определяется до значения $rapprox 3 $, затем цикл с периодом 2 до $rapprox 3,4$, затем цикл с периодом 4, затем цикл с периодом 8 и далее с быстрым приближением к хаосу.

Следует отметить, что вертикальные области незаполненного пространства на графике- это области r=3,6 и r=3,7, между r=3,7 и r=3,8, между r=3,8 и r=3,9, куда возвращается циклический порядок из предыдущего хаоса.
Для рассмотрения появления цикла с периодом кратным 3 в области $3,8leqslant rleqslant 3,9 $ внесём изменения в предыдущую программу:

Код программы

 # -*- 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()

Результат выполнения программы:

Цикл периода 3 появляется в окрестности точки r=3,83, а затем разбивается последовательно на циклы 6,12,24. Существование цикла с периодом 3 подразумевает наличие циклов любого другого конечного периода, а так же хаотических циклов вообще без периода.

Диаграмма ветвления позволяет проследить за развитием системы при плавном изменении параметра. При фиксированном значении параметра за орбитами точек позволяет проследить паутинная диаграмма (диаграмма Ламерея).
Построение паутинной диаграммы позволяет выявить различные эффекты, незаметные на диаграмме ветвления. Напишем программу:

Код программы

 # -*- 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()

Диаграмма Ламерея:
Диаграмма Ламерея

Удвоение периода в механических системах

Рассмотрим дифференциальное уравнение, которое моделирует свободные затухающие колебания материальной точки заданной массы на нелинейной пружине, при которых затухание определяется скоростью.

$mx^{''}+cx^{'}+kx+beta x^{3}=0$ (6)

В уравнении (6) член kx представляет силу линейной пружины, приложенной к материальной точке заданной массы, а член $beta x^{3}$ представляет фактическую нелинейность пружины.

Если на систему свободных колебаний (6) действует сила, то перемещение материальной точки массы, к которой приложена эта сила, описывается дифференциальным уравнением Дуффинга для вынужденных колебаний:

$mx^{''}+cx^{'}+kx+beta x^{3}=F_{0}cosomega t$ (7)

Уравнение (7) для большинства входящих в него параметров решается численным методом. Механическая система для математической модели по уравнению (7) приведена на рисунке:
image

Особенностью приведенной системы является то, что вместо пружины используется гибкая металлическая нить, которая колеблется в вертикальной плоскости, для которой константа Гука k отрицательна. В этой схеме точки устойчивого равновесия (а) и (с), а точка неустойчивого равновесия (b).

При смещении материальной точки из положения (b), действующая на неё сила является отталкивающей. Если периодическая сила, например, созданная осциллирующим магнитным полем частично гасится сопротивлением воздуха. Тогда, уравнение (7) является приемлемой математической моделью для горизонтального перемещения x(t) материальной точки при следующих областях параметров $k<0,c>0,beta >0$.

Для исследования поведения такой нелинейной системы примем $k=-1,m=c=beta =omega =1$, тогда дифференциальное уравнение (7) принимает вид:

$x^{''}+x^{'}-x+x^{3}=F_{0}cos(t)$, (8)

Напишем программу численного интегрирования уравнения (8) при начальных условиях $x(0)=1,x^{'}(0)=0$ в области $100leqslant tleqslant 200$ и для каждого из следующих значений амплитуды $F_{0}=0,6;0,7;0,75;0,8$, причем в каждом случае вывести на график решения для плоскостей $x(t),x^{'}(t)$ и $t,x(t)$:

Код программы

 # -*- 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=80,Хаос')                 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=80,Хаос')                 plt.legend(loc='upper left')                 plt.grid(True)                 plt.show()

Графики как результат работы программы
image
image

Этот переход от удвоения периода к хаосу показывает общий характер поведения нелинейной механической системы в ответ на изменение соответствующего физического параметра, например: $k,m,c,beta ,omega ,F_{0}$. Такие явления не происходят в линейных механических системах.

Аттрактор Лоренца

Подстановка в уравнение Дуффинга для вынужденных колебаний (7) приводят к двумерной нелинейной системе дифференциальных уравнений, что и было приведено в предыдущем листинге. Трёхмерную нелинейную систему дифференциальных уравнений применительно к задачам метеорологии рассматривал Э.Н. Лоренц:

$frac{dx}{dt}=-sx+sy,$
$frac{dy}{dt}=-xz+rx-y,$ (9)
$frac{dz}{dt}=xy-dz$

Решение системы (9) лучше рассматривать в проекции на одну из трёх плоскостей. Напишем программу численного интегрирования при значениях параметров b=frac{8}{3},s=10,r=28 и начальных условиях x(0)=-8, y(0)=8, z(0)=27:

Код программы

 # -*- 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()

Результат работы программы:
image

Рассматривая изображение на графике во времени, можно предположить, что точка P(x(t), y{t), z(t)) совершает случайное число колебаний то справа, то с слева. Для метеорологического приложения системы Лоренца,
после случайного числа ясных дней, следует случайное число дождливых дней.

Рассмотрим программу для отображения аттрактора Лоренца в плоскости xyz для мало различающихся начальных условий:

Код программы

 # -*- 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()

Результаты работы программы показаны на следующих графиках:
image
image

Из приведенных графиков следует, что изменение начального условия для с 1,0 до 1,0001 резко меняет характер изменения аттрактора Лоренца.

Система Росслера

Это очень интенсивно изучаемая нелинейная трехмерная система:
$frac{dx}{dt}=-y-z,$
$frac{dy}{dt}=x-alpha y,$ (10)
$frac{dz}{dt}=b+x(x-c).$

Напишем программу для численного интегрирования системы (10) для следующих параметров a=0,39, b=2, c=4 при начальных условиях x(0)=0, y(0)=0, z(0)=0:

Код программы

 # -*- 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()

Результат работы программы:
График

В плоскости лента Росслера выглядит как петля, но в пространстве она оказывается перекручена подобно ленте Мебиуса.

Выводы

Для демонстрации явлений хаоса приведены простые и интуитивно понятные программы на высокоуровневом языке программирования Python, которые легко модернизировать под новые проекты по данной теме. Статья имеет учебно-методическую направленность и может быть использована в процессе обучения.

Ссылки

  1. Немного о хаосе и о том, как его сотворить
  2. Критический взгляд на аттрактор Лоренца
  3. Генераторы хаоса на ПЛИС
  4. Дифференциальные уравнения и краевые задачи: моделирование и вычисление с помощью Mathematica, Maple и MATLAB. 3-е издание.: Пер. с англ. — М.: ООО «И.Д. Вильяме», 2008. — 1104 с.: ил. — Парал. тит. англ.

 
Источник

python, аттрактор Лоренца, диаграмма ветвления, теория хаоса, удвоение периода

Читайте также

Комментарии для сайта Cackle

Паблик ВКонтакте