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

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

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

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

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

Эту итерационную функцию можно интерпретировать как зависимость численности популяции от ее величины в предыдущий период времени и от параметра r, который обычно называют скоростью размножения. Собственно, сама по себе функция довольно унылая и имеет весьма банальный график. Интересные вещи проявляются, если рассматривать ее бифуркационную диаграмму: изменяя параметр r, можно наблюдать «динамику» неподвижных точек уравнения. Запишем логистическое отображение в 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) 

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

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

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

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

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

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

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, можно получить серию изображений, на которых видно, как происходит эволюция системы от начальных условий и неподвижной точки к, собственно, странному аттрактору.

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», говорит, что, хотя в хаотичных системах трудно предсказать состояние какой-то конкретной точки, статистические же параметры таких систем вполне точно определены. Таким образом, можно утверждать о статистической предсказуемости системы. Например, погода в конкретную минуту, строго говоря, непредсказуема, зато климат имеет вполне себе определенные параметры. И в этом нет никакого противоречия.