Облик странных аттракторов необычен и притягателен даже в двумерном измерении. Plotly позволяет строить их в трех измерениях, причем он дает возможность очень легко получить именно 3D-модель, которую можно «вертеть» и сквозь которую можно «пролетать» — ощущение «прикосновения».
С чего все началось
Все началось очень давно, где-то в 2007 году в университете я познакомился с теорией самоорганизации и впервые увидел аттрактор Лоренца, его черно-белую иллюстрацию в какой-то книге. Тогда мне показалось слишком странным, что нечто может двигаться по столь необычной траектории. Еще более странным мне казалась сама идея того, что практически все на свете может быть описано одной-единственной теорией.
В общем, все как обычно — мое мировоззрение поменялось, жизнь продолжалась, время шло. И вот, совсем недавно я натыкаюсь на ссылку и вижу вот это:
Изображение взято с сайта chaoticatmospheres.com
«Красиво.» — подумал я. Мысль о том, что все это можно построить в Matplotlib тоже промелькнула, но я уже заранее знал, что ничего впечатляющего не получится. И вот совсем недавно, буквально две недели назад, я познакомился с Plotly и сразу понял, что из этого может что-то получиться.
Первая попытка построения сразу же провалилась. Оказалось, что формулы на некоторых изображениях «Галереи странных аттракторов» содержат ошибки. Впрочем, автор галереи, честно предупреждает о том, что не является математиком, как и автор этой статьи.
Недолгое «гугление» позволило найти вот этот код, который оказался чрезвычайно полезен, а создал его Michael Tyka. Этот замечательный человек сделал целый плагин для Blender-а, позволяющий строить модели (!) 60-ти аттракторов. Фактически, их можно распечатать на 3D-принтере, а учитывая, что есть технологии печати воском, то можно довольно легко получить форму для отливки в бронзе.
Код для визуализации
Что ж, помимо, того что я являюсь математиком-любителем, я так же являюсь любителем-программистом. Так что не судите строго за качество кода.
################################
### ИМПОРТИРУЕМ БИБЛИОТЕКИ ###
################################
import numpy as np
from scipy.integrate import odeint
import plotly.graph_objects as go
##################################
### РЕШАЕМ СИСТЕМУ УРАВНЕНИЙ ###
##################################
# Система уравнений:
def LorenzMod1(XYZ, t, alpha, beta, xi, delta):
x, y, z = XYZ
x_dt = -alpha*x + y*y - z*z + alpha*xi
y_dt = x*(y - beta*z) + delta
z_dt = -z + x*(beta*y + z)
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.1
beta = 4
xi = 14
delta = 0.08
x_0, y_0, z_0 = 0, 1, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 100, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(LorenzMod1, (x_0, y_0, z_0), t, args=(alpha, beta, xi, delta))
X, Y, Z = f.T
#######################
### ВИЗУАЛИЗИРУЕМ ###
#######################
# Массив, отвечающий за изменение цвета:
c = np.linspace(0, 1, n)
# Готовим данные для отрисовки:
DATA = go.Scatter3d(x=X, y=Y, z=Z,
line=dict(color= c,
width=3,
# Выбираем цветовую палитру:
# Greys,YlGnBu,Greens,YlOrRd,Bluered,RdBu,
# Reds,Blues,Picnic,Rainbow,Portland,Jet,
# Hot,Blackbody,Earth,Electric,Viridis,Cividis.
colorscale="Cividis"),
# Рисуем только линии:
mode='lines')
fig = go.Figure(data=DATA)
# Задаем параметры отрисовки:
fig.update_layout(width=1000, height=1000,
margin=dict(r=10, l=10, b=10, t=10),
# Устанавливаем цвет фона:
paper_bgcolor='rgb(0,0,0)',
scene=dict(camera=dict(up=dict(x=0, y=0, z=1),
eye=dict(x=0, y=1, z=1)),
# Устанавливаем пропорциональное
# соотношение осей друг к другу:
aspectratio = dict(x=1, y=1, z=1),
# Отображаем, как указано в "aspectratio"
aspectmode = 'manual',
# Скрываем оси:
xaxis=dict(visible=False),
yaxis=dict(visible=False),
zaxis=dict(visible=False)
)
)
######################
#!! ВОСТОРГАЕМСЯ !!#
######################
fig.show()
В результате должна появиться 3D-модель странного аттрактора, называемого Lorenz Mod 1:
Нужно сразу отметить, что для решения систем дифференциальных уравнений была выбрана функция odeint из модуля SciPy, которая показалась мне самым простым и быстрым выходом для создания работающего кода. Однако, все уравнения могут решаться обычным методом Эйлера.
Для обозначения коэффициентов в коде я по привычке использовал названия греческих букв принятых в LaTeX-е. В работе с блокнотами Jupyter это иногда бывает очень полезно, так-как формулы могут быстро стать кодом, а код может быстро превратиться в формулы.
Если вы новичок в экосистеме Python, но хотите, что бы код гарантированно выполнился, то лучше всего установите последнюю версию дистрибутива Python Anaconda, а пакет Plotly через conda — встроенный менеджер пакетов дистрибутива.
Учитывая огромное количество странных аттракторов, построить их все, мне не представляется возможным. Поэтому, в этой статье я приведу лишь самые интересные из тех, что мне удалось построить.
The Chen-Lee Attractor
# Система уравнений:
def ChenLee(XYZ, t, alpha, beta, delta):
x, y, z = XYZ
x_dt = alpha*x - y*z
y_dt = beta*y + x*z
z_dt = delta*z + x*y/3
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 5
beta = -10
delta = -0.38
x_0, y_0, z_0 = 1, 1, 1
# Максимальное время и общее количество
# временных точек:
tmax, n = 200, 30000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(ChenLee, (x_0, y_0, z_0), t,
args=(alpha, beta, delta))
The Chua Attractor
# Система уравнений:
def ChuaAttractor(XYZ, t, alpha, beta, zeta, delta):
x, y, z = XYZ
h = zeta*x + (0.5*(delta - zeta))*(np.abs(x + 1) - np.abs(x - 1))
x_dt = alpha*(-x + y - h)
y_dt = x - y + z
z_dt = -beta*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 15.6
beta = 25.58
zeta = -5/7
delta = -8/7
x_0, y_0, z_0 = 1.8, -0.7, -2.85
# Максимальное время и общее количество
# временных точек:
tmax, n = 200, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(ChuaAttractor, (x_0, y_0, z_0), t,
args=(alpha, beta, zeta, delta))
The Coullet Attractor
# Система уравнений:
def Coullet(XYZ, t, alpha, beta, zeta, delta):
x, y, z = XYZ
x_dt = y
y_dt = z
z_dt = alpha*x + beta*y + zeta*z + delta*x**3
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.8
beta = -1.1
zeta = -0.4
delta = -1
x_0, y_0, z_0 = 0.1, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 200, 20000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(Coullet, (x_0, y_0, z_0), t,
args=(alpha, beta, zeta, delta))
The Dadras Attractor
# Система уравнений:
def DadrasAttractor(XYZ, t, rho, sigma, tau, zeta, epsilon):
x, y, z = XYZ
x_dt = y - rho*x + sigma*y*z
y_dt = tau*y - x*z + z
z_dt = zeta*x*y - epsilon*z
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
rho = 3
sigma = 2.7
tau = 1.7
zeta = 2
epsilon = 9
x_0, y_0, z_0 = 0.1, 0.03, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 220, 40000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(DadrasAttractor, (x_0, y_0, z_0), t,
args=(rho, sigma, tau, zeta, epsilon))
The Dequan Li Attractor
# Система уравнений:
def DequanLi(XYZ, t, alpha, beta, delta, epsilon, rho, xi):
x, y, z = XYZ
x_dt = alpha*(y - x) + delta*x*z
y_dt = rho*x + xi*y -x*z
z_dt = beta*z + x*y - epsilon*x*x
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 40
beta = 1.833
delta = 0.16
epsilon = 0.65
rho = 55
xi = 20
x_0, y_0, z_0 = 0.01, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 50, 40000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(DequanLi, (x_0, y_0, z_0), t,
args=(alpha, beta, delta, epsilon, rho, xi))
The Finance Attractor
# Система уравнений:
def FinanceAttractor(XYZ, t, alpha, beta, zeta):
x, y, z = XYZ
x_dt = (1/beta - alpha)*x + x*y + z
y_dt = -beta*y - x**2
z_dt = -x - zeta*z
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.001
beta = 0.2
zeta = 1.1
x_0, y_0, z_0 = 0.1, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 300, 40000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(FinanceAttractor, (x_0, y_0, z_0), t,
args=(alpha, beta, zeta))
The Four-Wing Attractor
# Система уравнений:
def FourWing(XYZ, t, alpha, beta, zeta):
x, y, z = XYZ
x_dt = alpha*x + y + y*z
y_dt = -x*z + y*z
z_dt = -z - zeta*x*y + beta
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 5
beta = 16
zeta = 2
x_0, y_0, z_0 = 1, -1, 1
# Максимальное время и общее количество
# временных точек:
tmax, n = 100, 60000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(FourWing, (x_0, y_0, z_0), t,
args=(alpha, beta, zeta))
The Hadley Attractor
# Система уравнений:
def HadleyAttractor(XYZ, t, alpha, beta, xi, delta):
x, y, z = XYZ
x_dt = -y*y - z*z - alpha*x + alpha*xi
y_dt = x*y - beta*x*z - y + delta
z_dt = beta*x*y + x*z-z
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.2
beta = 4
xi = 8
delta = 1
x_0, y_0, z_0 = 0.39, -1, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 100, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(HadleyAttractor, (x_0, y_0, z_0), t,
args=(alpha, beta, xi, delta))
The Halvorsen Attractor
# Система уравнений:
def HalvorsenAttractor(XYZ, t, alpha):
x, y, z = XYZ
x_dt = -alpha*x - 4*y - 4*z - y*y
y_dt = -alpha*y - 4*z - 4*x - z*z
z_dt = -alpha*z - 4*x - 4*y - x*x
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 1.4
x_0, y_0, z_0 = -5, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 100, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(HalvorsenAttractor, (x_0, y_0, z_0), t,
args=(alpha,))
The Liu-Chen Attractor
# Система уравнений:
def LiuChen(XYZ, t, alpha, beta, sigma, delta, epsilon, xi):
x, y, z = XYZ
x_dt = alpha*y + beta*x + sigma*y*z
y_dt = delta*y - z + epsilon*x*z
z_dt = xi*z - x*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 2.4
beta = -3.75
sigma = 14
delta = -11
epsilon = 4
xi = 5.58
x_0, y_0, z_0 = 1, 3, 5
# Максимальное время и общее количество
# временных точек:
tmax, n = 55, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(LiuChen, (x_0, y_0, z_0), t,
args=(alpha, beta, sigma, delta, epsilon, xi))
The Lorenz Mod 2 Attractor
# Система уравнений:
def LorenzMod2(XYZ, t, alpha, beta, xi, delta):
x, y, z = XYZ
x_dt = -alpha*x + y**2 -z**2 + alpha*xi
y_dt = x*(y - beta*z) + delta
z_dt = -z + x*(beta*y + z)
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.9
beta = 5
xi = 9.9
delta = 1
x_0, y_0, z_0 = 5, 5, 5
# Максимальное время и общее количество
# временных точек:
tmax, n = 50, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(LorenzMod2, (x_0, y_0, z_0), t,
args=(alpha, beta, xi, delta))
The Modified Chua Chaotic Attractor
# Система уравнений:
def ChuaModified(XYZ, t, alpha, beta, gamma, delta, zeta):
x, y, z = XYZ
h = -delta*np.sin((np.pi*x)/(2*gamma))
x_dt = alpha*(y - h)
y_dt = x - y + z
z_dt = -beta*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 10.82
beta = 14.286
gamma = 1.3
delta = 0.11
zeta = 7
x_0, y_0, z_0 = 1, 1, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 200, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(ChuaModified, (x_0, y_0, z_0), t,
args=(alpha, beta, gamma, delta, zeta))
The Newton Leipnik Attractor
# Система уравнений:
def NewtonLeipnik(XYZ, t, alpha, beta):
x, y, z = XYZ
x_dt = -alpha*x + y + 10*y*z
y_dt = -x - 0.4*y + 5*x*z
z_dt = beta*z - 5*x*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.4
beta = 0.175
x_0, y_0, z_0 = 0.349, 0, -0.16
# Максимальное время и общее количество
# временных точек:
tmax, n = 300, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(NewtonLeipnik, (x_0, y_0, z_0), t,
args=(alpha, beta))
The Nose-Hoover Attractor
# Система уравнений:
def NoseHoover(XYZ, t, alpha):
x, y, z = XYZ
x_dt = y
y_dt = -x + y*z
z_dt = alpha - y*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 1.5
x_0, y_0, z_0 = 1, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 150, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(NoseHoover, (x_0, y_0, z_0), t,
args=(alpha,))
The Roessler Attractor
# Система уравнений:
def Roessler(XYZ, t, alpha, beta, sigma):
x, y, z = XYZ
x_dt = -(y + z)
y_dt = x + alpha*y
z_dt = beta + z*(x - sigma)
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.2
beta = 0.2
sigma = 5.7
x_0, y_0, z_0 = 1, 1, 1
# Максимальное время и общее количество
# временных точек:
tmax, n = 300, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(Roessler, (x_0, y_0, z_0), t,
args=(alpha, beta, sigma))
The Sakarya Attractor
# Система уравнений:
def SakaryaAttractor(XYZ, t, alpha, beta):
x, y, z = XYZ
x_dt = -x + y + y*z
y_dt = -x - y + alpha*x*z
z_dt = z - beta*x*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.4
beta = 0.3
x_0, y_0, z_0 = 1, -1, 1
# Максимальное время и общее количество
# временных точек:
tmax, n = 100, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(SakaryaAttractor, (x_0, y_0, z_0), t,
args=(alpha, beta))
The Thomas Attractor
# Система уравнений:
def Thomas(XYZ, t, beta):
x, y, z = XYZ
x_dt = -beta*x + np.sin(y)
y_dt = -beta*y + np.sin(z)
z_dt = -beta*z + np.sin(x)
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
beta = 0.19
x_0, y_0, z_0 = 0.1, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 185, 10000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(Thomas, (x_0, y_0, z_0), t,
args=(beta,))
The Three-Scroll Unified Chaotic System Attractor (TSUCS1)
# Система уравнений:
def TSUCS1(XYZ, t, alpha, beta, delta, epsilon, xi):
x, y, z = XYZ
x_dt = alpha*(y - x) + delta*x*z
y_dt = xi*y - x*z
z_dt = beta*z + x*y - epsilon*x*x
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 40
beta = 0.833
delta = 0.5
epsilon = 0.65
xi = 20
x_0, y_0, z_0 = 0.01, 0, 0
# Максимальное время и общее количество
# временных точек:
tmax, n = 70, 50000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(TSUCS1, (x_0, y_0, z_0), t,
args=(alpha, beta, delta, epsilon, xi))
The Wang-Sun Attractor
# Система уравнений:
def WangSunAttractor(XYZ, t, alpha, beta, zeta, delta, epsilon, xi):
x, y, z = XYZ
x_dt = alpha*x + zeta*y*z
y_dt = beta*x + delta*y - x*z
z_dt = epsilon*z + xi*x*y
return x_dt, y_dt, z_dt
# Параметры системы и начальные условия:
alpha = 0.2
beta = -0.01
zeta = 1
delta = -0.4
epsilon = -1
xi = -1
x_0, y_0, z_0 = 0.5, 0.1, 0.1
# Максимальное время и общее количество
# временных точек:
tmax, n = 500, 30000
# Интегрируем систему уравнений в каждой точке
# временного интервала t:
t = np.linspace(0, tmax, n)
f = odeint(WangSunAttractor, (x_0, y_0, z_0), t,
args=(alpha, beta, zeta, delta, epsilon, xi))
В заключение
Огонь, вода, земля, небо, солнце, луна, звезды — все это самые древние поэтические сущности. Очень часто мне удается найти в математике что-то такое же прекрасное. Но гораздо чаще я даже не понимаю, как обо всем этом говорить на математическом языке и языке обычном. Не понимаю, но хочется научиться.
Но вот что я осознал на все 100%, так это то, что современные инструменты визуализации дают фантастическую возможность выразить свое отношение к тому чем ты сейчас занят, возможность показать, как это важно для тебя и как тебе это интересно. Сделать все это без слов.
Комментарии (24)
LibrarianOok
18.12.2019 11:20+1Шикарно! Попробую перетолмачить это на Lua для движка Instead. Спасибо!
veledzimovich
18.12.2019 03:15Добрый день. А вы можете пояснить почему собираетесь использовать Instead, а не Love2d?
LibrarianOok
19.12.2019 10:01Потому, что мне нравится API этого движка, а Love2d — скорее фреймворк для написания поверх него своего движка. В Instead можно без мороки раскидать визуализации по разным комнатам или объектам и получить удобное меню для навигации. Дополнить текстовыми пояснениями, если надо. Кроме того, у Instead есть репозиторий, куда можно будет эту программу выложить и она точно не затеряется среди всех этих сотен энтертейнмента на итче, стиме, гоге и т.д.
Smirnov_Nick
18.12.2019 16:47+6Класс! Немного покопался и смог нарисовать аттрактор Aizawa. Спасибо автору за статью!
WinterSt
18.12.2019 04:19Не подскажите как вы сделали анимацию?
Smirnov_Nick
18.12.2019 12:58Вращаю камеру вокруг оси Z и с помощью imageio генерирую серию изображений, в конце объединяю изображения в видео.
idimus
18.12.2019 18:35Есть еще такая программулина как «Chaoscope», можно подглядеть режимы визуализации, для еще более красивых картинок.
uchitel Автор
18.12.2019 05:29Да, занимательная программа. У моей визуализации есть большой минус, который портит все впечатление — приходится каждый раз возвращаться в код, что бы изменить какой-нибудь параметр. Можно добавить виджетов и сделать интерфейс подобный «Chaoscope». Даже не знаю почему мне эта идея не пришла раньше. Но, насколько я помню, в «Chaoscope» изображение является статичным, а в Plotly его можно «вертеть» — мне это очень сильно понравилось.
Arneodo Attractor:
axi1
18.12.2019 12:36+1Там в chaoscope, кроме слайдеров, ещё можно рендерить результат с повышенным качеством:
Медленно, правда, он рендерит.
idimus
19.12.2019 00:30Да, нельзя + он уже 10-ть лет, как не разрабатывается, но ведь ничто не мешает взять лучшее)
SMC
18.12.2019 21:07Боже мой, это великолепно! Спасибо, за то, что напомнили об этом мире. (удивительном мире и самой поэзии математики)
Как бы хотелось иметь как минимум ещё одну жизнь, что бы посветить её познанию именно её (математики) безграничных горизонтов!
Она потрясающа!
Спасибо!uchitel Автор
18.12.2019 05:56Да, математикой можно наслаждаться. Даже не помню, когда в последний раз смотрел какой-нибудь фильм или сериал. Все свободное время — это музыка, крепкий чай и математика.
QDeathNick
18.12.2019 00:06+1В школе на спеки рисовал подобные вещи. Вот бы тому мне из прошлого тысячелетия рассказать, что оказывается это были аттракторы. Был бы отличный эффект бабочки.
vitalijlysanov
18.12.2019 02:10+1
Если запускать код по методике:
https://habr.com/ru/post/448316/
то все уже установленоuchitel Автор
18.12.2019 04:34Да, в Google Colaboratory все уже установлено — можно с ходу начинать работать. Я как-то забыл об этом. Спасибо, что напомнили.
axi1
18.12.2019 08:57+1Я как раз такое задание давал на экзамене, только мы делали сетку начальных условий и со слайдерами.
Видео: www.youtube.com/watch?v=B9-sxUs-3wo.
Тоже питон, но не в браузере: pyqt, pyqtgraph, numba.
ivantgam
Это великолепно, спасибо!