ジュリア集合
$z, c \in \mathbf{C}$に対して、漸化式$z_{n+1} = z_{n}^{2} + c$を考える。この漸化式に任意の$z_{0}$を代入して計算していくと発散するものと有界なものに分かれる。有界なものの集合を$K_{c}$、$K_{c}$の境界$\partial K_{c}$を$J_{c}$とする。この境界$J_{c}$がジュリア集合である。(2019/2/18 追記: $z_{0} = 0$のものをマンデルブロ集合という)
プロット
ジュリア集合を具体的に計算するのは難しいので(多分不可能)、格子点について漸化式を計算して一定の回数以内に発散すると確定しなかったものをプロットすることにする。こうすればジュリア集合とその内部をプロットすることになる。$|z_{0}| > 2 \max(|c|, 1)$で発散するのが明らか(2019/2/18 なぜそう思ったのか思い出せないので参考文献を追記)なのでこれを条件とする。
size <- 1000 x <- y <- seq(-2, 2, length=size) z <- outer(x, y, function(a, b){a + b * 1i}) for (c in c(0.285 + 0.01i, -0.618 + 0.0i, -0.8 + 0.156i)) { julia <- matrix(rep(1, size ^ 2), nc=size) tmp <- z for (i in 1:100) { tmp <- tmp ^ 2 + c julia[which(Mod(tmp) > 2)] <- 0 tmp[which(Mod(tmp) > 2)] <- NaN } plot(z[which(julia == 1)], pch=20, xlab='', ylab='', main=as.character(c), col="blue") }
以前contourプロットもやったので何回目に発散したかもプロットしてみる。
size <- 1500 x <- y <- seq(-2, 2, length=size) z <- outer(x, y, function(a, b){a + b * 1i}) for (c in c(0.285 + 0.01i, -0.618 + 0.0i, -0.8 + 0.156i)) { julia <- matrix(rep(1, size ^ 2), nc=size) tmp <- z for (i in 1:500) { julia[which(Mod(tmp) > 2)] <- log(i - 1) tmp[which(Mod(tmp) > 2)] <- NaN tmp <- tmp ^ 2 + c } image(x, y, julia, xlab='', ylab='', col=colorRampPalette(c("white", "yellow", "red"), space = "rgb")(10), main=as.character(c)) }
追記
- How To Quickly Compute The Mandelbrot Set In Python (IT Best Kept Secret Is Optimization): マンデルブロ集合を描く例。numpyを使った普通のコードからcythonやopencl、tensorflowで最適化した例まである。
- マンデルブロ集合とは - マンデルブロ集合の不思議な世界: $|z_{n}| > 2$で発散することの証明あり