読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

ジュリア集合をプロットしてみる

R 力学系

ジュリア集合

$z, c \in \mathbf{C}$に対して、漸化式$z_{n+1} = z_{n} + c$を考える。この漸化式に任意の$z_{0}$を代入して計算していくと発散するものと有界なものに分かれる。有界なものの集合を$K_{c}$、$K_{c}$の境界$\partial K_{c}$を$J_{c}$とする。この境界$J_{c}$がジュリア集合である。

プロット

ジュリア集合を具体的に計算するのは難しいので(多分不可能)、格子点について漸化式を計算して一定の回数以内に発散すると確定しなかったものをプロットすることにする。こうすればジュリア集合とその内部をプロットすることになる。$|z_{0}| > 2 \max(|c|, 1)$で発散するのが明らかなのでこれを条件とする。

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")
}

f:id:kuyata:20140801145750p:plain

f:id:kuyata:20140801145804p:plain

f:id:kuyata:20140803183007p:plain

以前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))
}

f:id:kuyata:20140803182804p:plain

f:id:kuyata:20140803182822p:plain

f:id:kuyata:20140803182826p:plain

追記