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