これがカオス!?
単純な方程式の挙動をRで可視化し、カオスを体感してみましょう。
実数\(r\)を変化させた際の次の方程式の挙動を調べてみましょう。
\[ x_{t+1}=r x_t (1- x_t) \tag{1} \]
準備
使うパッケージを読み込みましょう。
library(tidyverse)
library(patchwork)必要に応じて、事前にインストールしておいてください。
install.packages("tidyverse")
install.packages("patchwork")特定のrでの挙動を可視化する
この式の挙動を視覚的に理解できるように可視化してみましょう。
\(x_t\)を代入して\(x_{t+1}\)を得る関数fを用意します。
f <- function(x, r){ r * x * (1 - x) }適当な初期値を用意します。
x0 <- 0.9
r <- 3.5100 stepまで計算する事にします。
t_max <- 100逐次計算はforで回してしまいましょう。そのための関数を用意します。
f_simulation <- function(r, x0, t_max){
x_all <- NULL
x <- x0
for(t in 1:t_max){
x <- f(x, r) # 計算して
x_all <- c(x_all, x) # 繋げる
}
return(x_all)
}結果をテーブルデータの形式に整えて取り扱いやすいtibble(data.frameの拡張版)に変換する関数を用意します。
f_arrange <- function(x_all){
t_max <- length(x_all)
dat <-
x_all %>%
tibble(t = 2:(t_max+1), x_pre = .) %>%
bind_rows(tibble(t = 1, x_pre = x0)) %>%
arrange(t) %>%
mutate(x_post = lead(x_pre)) %>% # 1つずらしてx_{n+1}を作る
filter(!is.na(x_post)) # NA行の除去
return(dat)
}さてお試し。
dat <-
f_simulation(r, x0, t_max) %>%
f_arrange()
dat %>% head(5)## # A tibble: 5 × 3
## t x_pre x_post
## <dbl> <dbl> <dbl>
## 1 1 0.9 0.315
## 2 2 0.315 0.755
## 3 3 0.755 0.647
## 4 4 0.647 0.799
## 5 5 0.799 0.561x_pre列が\(x_t\)、x_post列が\(x_{t+1}\)に相当します。
これを可視化します。まず、横軸に\(t\)、縦軸に\(x_t\)を取ったグラフを作成します。この作図関数を後でも使いたいので独自の関数gg_line()として定義しておきます。
# 作図関数
gg_line <- function(dat, t_upper = 100){
ggplot(data = dat %>% filter(t <= t_upper)) +
aes(x = t, y = x_pre) +
geom_path() +
geom_point() +
scale_y_continuous(limits = c(0, 1)) +
ylab("x_t")
}
# 作図
dat %>% gg_line()
明らかに4つの値を周期的に取るパタンに収束することが分かります。
別の作図として、横軸に\(x_t\)、縦軸に\(x_{t+1}\)を取ってみましょう。
ggplot(data = dat) +
aes(x = x_pre, y = x_post) +
geom_point() +
scale_x_continuous(limits = c(0, 1)) +
scale_y_continuous(limits = c(0, 1)) +
xlab("x_t") + ylab("x_t+1")
二次関数上に点が分布しています。この関数は式(1)で表されるはずです。図に加えてみましょう。
曲線の座標を計算しておいて
dat_line <-
seq(0, 1, by = 0.01) %>%
tibble(x_pre = .) %>%
mutate(x_post = r * x_pre * (1 - x_pre))geom_line()関数で図に青色の線で重ねます。プロットが前面に描画されるようにラインの作図はgeom_point()関数の前に挿入します。
ggplot(data = dat) +
aes(x = x_pre, y = x_post) +
geom_path(data = dat_line, color = "blue") + # ←ここ
geom_point() +
scale_x_continuous(limits = c(0, 1)) +
scale_y_continuous(limits = c(0, 1)) +
xlab("x_t") + ylab("x_t+1")
良いですね。更にポイント同士を結び、step-by-stepで値が変化する様子を可視化してみましょう。ただし単に順番に点を結ぶのではなく、ある点が次のステップの点のx座標に利用されている様子が分かるように工夫しましょう。
これも後で使いまわしたいので、作図関数gg_chaos()を定義しておきます。
# 作図関数
gg_chaos <- function(dat, dat_line){
ggplot(data = dat) +
aes(x = x_pre, y = x_post) +
geom_path(data = dat_line, color = "blue") + # 二次曲線
# 最初(step 0)の縦線
geom_segment(data = dat %>% filter(t == 1),
aes(xend = x_pre, yend = 0), alpha = 0.4, linewidth = 0.2) +
# step 1 以降の縦線
geom_segment(data = dat %>% filter(t > 1),
aes(xend = x_pre, yend = x_pre), alpha = 0.4, linewidth = 0.2) +
# 横線
geom_segment(aes(xend = x_post, yend = x_post), alpha = 0.4, linewidth = 0.2, linetype = "dashed") +
# データポイント
geom_point()+
# スケール調整
scale_x_continuous(limits = c(0, 1)) +
scale_y_continuous(limits = c(0, 1)) +
# 軸ラベル
xlab("x_t") + ylab("x_t+1")
}
# 作図
dat %>%
gg_chaos(dat_line)
良い感じですね。
もう少し工夫すると次のような動画を作ることができます。
rを変化させる
例えば\(r \in (1.5, 2, 2.5, 3, 3.5, 4)\)とします。その時の応答を効率的に計算するには、畳み込みデータ(nested data)を上手く利用しましょう。
cf. map脳になろう、もしくはnested dataのハンドリング
dat <-
seq(1.5, 4, by = 0.5) %>%
tibble(r = .) %>%
mutate(data = map(
r, ~ f_simulation(., x0, t_max) %>% f_arrange())
) %>%
mutate(data_line = map(
r, ~ tibble(x_pre = seq(0, 1, by = 0.01)) %>% mutate(x_post = .x * x_pre * (1 - x_pre)))
)
dat## # A tibble: 6 × 3
## r data data_line
## <dbl> <list> <list>
## 1 1.5 <tibble [100 × 3]> <tibble [101 × 2]>
## 2 2 <tibble [100 × 3]> <tibble [101 × 2]>
## 3 2.5 <tibble [100 × 3]> <tibble [101 × 2]>
## 4 3 <tibble [100 × 3]> <tibble [101 × 2]>
## 5 3.5 <tibble [100 × 3]> <tibble [101 × 2]>
## 6 4 <tibble [100 × 3]> <tibble [101 × 2]>このdata列に格納されているのが、それぞれのrに対応するシミュレーションの結果です。data_lineが青色の曲線のための数値ですね。
このまま先ほど作ったgg_chaos()関数で作図してしまいます。
dat_g <-
dat %>%
mutate(g = map2(data, data_line, ~ gg_chaos(.x, .y))) %>%
mutate(g = map2(g, r, ~ .x + ggtitle(str_c("r=", .y))))
dat_g## # A tibble: 6 × 4
## r data data_line g
## <dbl> <list> <list> <list>
## 1 1.5 <tibble [100 × 3]> <tibble [101 × 2]> <gg>
## 2 2 <tibble [100 × 3]> <tibble [101 × 2]> <gg>
## 3 2.5 <tibble [100 × 3]> <tibble [101 × 2]> <gg>
## 4 3 <tibble [100 × 3]> <tibble [101 × 2]> <gg>
## 5 3.5 <tibble [100 × 3]> <tibble [101 × 2]> <gg>
## 6 4 <tibble [100 × 3]> <tibble [101 × 2]> <gg>このg列にグラフが入っているのでwrap_plots()関数を使って1つのグラフにまとめます。
dat_g$g %>% wrap_plots()
\(x_t\)の時系列も同じ方法で描画してみましょう。
dat_g %>%
mutate(g_line = map(data, gg_line)) %>%
mutate(g_line = map2(g_line, r, ~ .x + ggtitle(str_c("r=", .y)))) %>%
.$g_line %>%
wrap_plots(ncol = 1)
\(r \le 3\)の範囲では\(x_t\)は1つの値に収束するのに対し、\(r\)が大きくなると取り得る値の範囲が広がることが分かります（ただし式(1)は常に満たされます）。
rに応じたxの収束性を可視化する
\(r\)をもっと細かく変化させた時に、step 300から500の\(x_t\)の取り得る値を可視化してみましょう。
dat <-
tibble(r = seq(2.5, 4, by = 0.005)) %>%
mutate(data = map(r, ~ f_simulation(., x0, t_max = 500) %>% f_arrange())) %>%
mutate(r = round(r, 3)) %>%
unnest(data)dat %>%
filter(t >= 300) %>%
ggplot() +
aes(r, x_pre) +
geom_vline(xintercept = seq(2.5, 4, by = 0.5), color = "lightblue") +
geom_point(size = 0.1) +
ylab("x_t")
美しいグラフが描けました。収束性は単純に拡大するのではなく、ところどころ極端な収束を見せることがあるようです。例えば\(r=3.8\)ぐらいのところで3つの値に収束する傾向があるようです。細かく値を変化させながら、この範囲の\(r\)の値がいくつかを調べてみると、\(r \simeq 3.83\)が該当しました。
dat %>%
filter(t >= 300) %>%
ggplot() +
aes(r, x_pre) +
geom_vline(xintercept = seq(2.5, 4, by = 0.5), color = "lightblue") +
geom_point(size = 0.1)+
geom_vline(xintercept = 3.83, color = "red", linewidth = 1)+
ylab("x_t")
そこで、その近傍で\(x_t\)の時系列プロットを描いてみると…
dat %>%
filter(t <= 150) %>%
filter(r %in% c(3.82, 3.84, 3.86)) %>%
ggplot() +
aes(t, x_pre) +
geom_point() +
geom_path() +
facet_wrap(~r, ncol = 1)+
ylab("x_t")
\(r\)の値を僅かに変更しているだけですが、\(r\)が特定の値になった時だけ突然に系全体が3つの値に収束する様子が見て取れます。
まとめ
式(1)のように非常に単純な系でも「僅かな初期値の違いで系全体の振る舞いが、予測できない変化を示す」というカオスを体験することができました。
補足
微分方程式をRで取り扱う方法は以下でも紹介しています。
謝辞
本記事は統計数理研究所の公開講座「カオス農学実践編―トラクタのカオス動力学と農作業安全対策―」にて講師の酒井憲司先生が取り上げていたトピックの一部をRで再現したものです。大変興味深いご講義をありがとうございました。感謝申し上げます。