Анимированные графики в R (и немного про бифуркацию, хаос и аттракторы)

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

install.packages("animation")
library(animation)

Стоит отметить, что animation в процессе обработки изображений использует пакет программ ImageMagick ^[1], так что его желательно установить заранее. Под Windows работоспособность этого решения я не проверял.
Для создания анимированного графика нам, в общем-то, понадобится всего одна функция такого вида:

saveGIF({
  # Какой-то код, генерирующий последовательность графиков
}, movie.name=..., interval=..., ani.width=...,  ani.height=...)

Так получилось, что в то время я проходил весьма познавательный курс Introduction to Dynamical Systems and Chaos ^[2], и мне было интересно, как от относительно простых математических объектов переходят к весьма причудливым изображениям. Взять хотя бы логистическое отображение ^[3] такого вида:

Эту итерационную функцию можно интерпретировать как зависимость численности популяции от ее величины в предыдущий период времени и от параметра r, который обычно называют скоростью размножения. Собственно, сама по себе функция довольно унылая и имеет весьма банальный график. Интересные вещи проявляются, если рассматривать ее бифуркационную диаграмму ^[4]: изменяя параметр r, можно наблюдать «динамику» неподвижных точек уравнения ^[5]. Запишем логистическое отображение в R в виде такой функции:

logistic.map <- function(r, x0, n, m){
  x <- rep(x0, n)
  for(i in 1:(n-1)) {
    x[i+1] <- r * x[i]  * (1 - x[i])
  }
  return(x[(n-m):n])
}

Зададим некоторые интересные точки и начальные параметры для отображения и построим бифуркационную диаграмму, изменяя масштаб изображения по ходу дела:

nrows <- 6
r.len <- 1500;
# Points of interest on the plot
R <- matrix(c(
  seq(2.4, 4, length.out=r.len),
  seq(3.442420, 3.639398, length.out=r.len),
  seq(3.562297, 3.572910, length.out=r.len),
  seq(3.569792, 3.570244, length.out=r.len),
  seq(3.570005, 3.571369, length.out=r.len),
  seq(3.631992, 3.633301, length.out=r.len)
), nrow=nrows, byrow=T)
X <- matrix(c(
  0, 1,
  0.8567335, 0.9140401,
  0.8887529, 0.8936790,
  0.8920580, 0.8925577,
  0.8911242, 0.8927333,
  0.9066966, 0.9083943
), nrow=nrows, byrow=T)

x0 <- 0.5
n <- 200
m <- 170
saveGIF({
  for (i in 1:nrows){
    r <- R[i,]
    x <- as.vector(sapply(r, logistic.map, x0, n, m))
    r <- sort(rep(r, (m+1)))
    del_idx <- unlist(sapply(1:length(x), function(j) if (x[j] < X[i, 1] | x[j] > X[i, 2]) j))
    if (length(del_idx > 0)){
      x <- x[-del_idx]
      r <- r[-del_idx]
    }
    plot(x ~ r, col="gray66", pch=".", main="Bifurcation Diagram for the Logistic Map")
  }
}, movie.name = "bifur.gif", interval=2.4, ani.width=600,  ani.height=500)

В результате получим примерно такую картинку:

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

Первое уравнение определяет расстояние от одной точки бифуркации до следующей в единицах r. Второе же показывает, насколько ответвление n длиннее ответвления n+1. Так вот, оказывается, что

Число 4.669201… называется универсальной постоянной Фейгенбаума ^[7], которая как раз и характеризует скорость перехода динамических систем от порядка к детерминированному хаосу.

Другой интересный и не менее известный объект — аттрактор Лоренца ^[8]. На Хабре ему даже посвящена отдельная статья ^[9]. Для описания движения воздушных потоков Эдвард Лоренц использовал систему из трех обыкновенных дифференциальных уравнений, известных теперь как уравнения Лоренца:

Не будем изголяться и решим эту систему численно — методом Эйлера ^[10]:

lorenz.solution <- function(sigma=10, r=28, beta=8/3, x=0.01, y=0.01, z=0.01, dt=0.001, n=30000){
  sol <- array(0, dim=c(n,3))
  t <- 0
  for(i in 1:n){
    x0 <- x; y0 <- y; z0 <- z
    x <- x0 + (y0 - x0) * sigma * dt
    y <- y0 + ((r - z0) * x0 - y0) * dt
    z <- z0 + (x0 * y0 - beta * z0) * dt
    t <- t + dt
    sol[i,] <- c(x, y, z)
  }
  return(sol)
}

Решение системы для параметров, заданных по умолчанию, выглядит так:

Изменяя параметр r, можно получить серию изображений, на которых видно, как происходит эволюция системы от начальных условий и неподвижной точки к, собственно, странному аттрактору ^[11].

library(scatterplot3d)
saveGIF({
  for (r in 2:34){ 
    sol <- lorenz.solution(r=r)
    s3d<-scatterplot3d(sol[,1], sol[,2], sol[,3], color="gray66", angle=15, box=F, grid=F, axis=F, pch=".", main=paste0("Lorenz Attractor with rho=", r))
  }
}, movie.name = "lorenz.gif", interval=.3, ani.width=500,  ani.height=500)

Аттрактора Лоренца не единственный в своем роде. Например, аттрактор Чена описывается так:

Скрытый текст

chen.solution <- function(a=40, c=28, b=3, x=-0.1, y=0.5, z=-0.6, dt=0.001, n=30000){
  sol <- array(0, dim=c(n,3))
  t <- 0
  for(i in 1:n){
    x0 <- x; y0 <- y; z0 <- z
    x <- x0 + (y0 - x0) * a * dt
    y <- y0 + ((c - a) * x0 - x0 * z0 + c * y0) * dt
    z <- z0 + (x0 * y0 - b * z0) * dt
    t <- t + dt
    sol[i,] <- c(x, y, z)
  }
  return(sol)
}

saveGIF({
  for (a in 32:45){ 
    sol <- chen.solution(a=a)
    s3d<-scatterplot3d(sol[,1], sol[,2], sol[,3], color="gray66", angle=15, box=F, grid=F, axis=F, pch=".", main=paste0("Chen Attractor with a=", a))
  }
}, movie.name = "chen.gif", interval=.25, ani.width=500,  ani.height=500)

Орбиты, по которым происходит движение, «стягиваются» к аттрактору; это движение хаотично и чувствительно к начальным условиям, в то же время оно стабильно в глобальном плане. Дэвид Фельдман, который ведет курс «Introduction to Dynamical Systems and Chaos», говорит, что, хотя в хаотичных системах трудно предсказать состояние какой-то конкретной точки, статистические же параметры таких систем вполне точно определены. Таким образом, можно утверждать о статистической предсказуемости системы. Например, погода в конкретную минуту, строго говоря, непредсказуема, зато климат имеет вполне себе определенные параметры. И в этом нет никакого противоречия.

Автор: kxx

Источник ^[12]