マンデルブロ集合の彩色方法（５）

1 発散判定式を変更する

mandelbrot() 内の発散判定式 \(|z_n| > 2\) を変更すると面白い模様になります。

1.1 虚数成分で判定する

発散判定式を \(|Im(z_n)| > 2\) に変更。

mandelbrot1 <- function(
                       center = -0.75+0.0i, ## 計算する範囲の中心
                       radius = 1.5,        ## 計算する範囲の半径
                       mesh.size = 720,     ## 計算する範囲の分割数
                       max.iteration = 50,  ## 反復回数上限
                       escape.radius = 2    ## 脱出半径
                       ) {
    require(tidyverse)

    x <- seq(Re(center)-radius, Re(center)+radius, length.out=mesh.size)
    y <- seq(Im(center)-radius, Im(center)+radius, length.out=mesh.size)
    m <- tibble(c=c(outer(x, y*1.0i, '+')), z=0.0, iter=0)

    for (i in 1:max.iteration) {
        # tmp <- Mod(m$z) < escape.radius
        tmp <- Im(m$z) * Im(m$z) < escape.radius * escape.radius
        m <- m %>%
            mutate(z=if_else(tmp, z^2+c, z), iter=if_else(tmp, iter+1, iter))
    }

    m <- m %>%
        mutate(iter=na_if(iter, max.iteration)) ## 反復回数上限に達した部分を NA に置き換え

    return(m)
}

library(tidyverse)

er <- 2*pi
m1 <- mandelbrot1(escape.radius=er)

### 滑らか coloring ###
m1 <- m1 %>%
    mutate(z1=z^2+c, x=er^2-Re(z)^2+Im(z)^2, y=Re(z1)^2+Im(z1)^2-er^2, f=x/(x+y), iter2=iter+f)

### Histogram coloring ###
h <- m1 %>%
    filter(!is.na(iter2)) %>% ## NA を除外
    count(iter2) %>%          ## iterの値毎に個数をカウント
    mutate(n=cumsum(n)) %>%   ## 累積和を求める
    mutate(n=n/max(n))        ## 値域を[0,1]にする

m1 <- m1 %>%
    left_join(h, by='iter2') %>%
    mutate(n=if_else(is.na(iter2), NA_real_, n))

ggplot(m1) +
    aes(x=Re(c), y=Im(c), fill=n) +
    geom_raster() +
    scale_fill_gradientn(colors=rainbow(13), na.value='black') +
    coord_fixed(ratio=1)

1.2 実数成分で判定する

発散判定式を \(|Re(z_n)| > 2\) に変更。

mandelbrot2 <- function(
                       center = -0.75+0.0i, ## 計算する範囲の中心
                       radius = 1.5,        ## 計算する範囲の半径
                       mesh.size = 720,     ## 計算する範囲の分割数
                       max.iteration = 50,  ## 反復回数上限
                       escape.radius = 2    ## 脱出半径
                       ) {
    require(tidyverse)

    x <- seq(Re(center)-radius, Re(center)+radius, length.out=mesh.size)
    y <- seq(Im(center)-radius, Im(center)+radius, length.out=mesh.size)
    m <- tibble(c=c(outer(x, y*1.0i, '+')), z=0.0, iter=0)

    for (i in 1:max.iteration) {
        # tmp <- Mod(m$z) < escape.radius
        tmp <- Re(m$z) * Re(m$z) < escape.radius * escape.radius
        m <- m %>%
            mutate(z=if_else(tmp, z^2+c, z), iter=if_else(tmp, iter+1, iter))
    }

    m <- m %>%
        mutate(iter=na_if(iter, max.iteration)) ## 反復回数上限に達した部分を NA に置き換え

    return(m)
}

library(tidyverse)

er <- 2*pi
m2 <- mandelbrot2(escape.radius=er)

### 滑らか coloring ###
m2 <- m2 %>%
    mutate(z1=z^2+c, x=er^2-Re(z)^2+Im(z)^2, y=Re(z1)^2+Im(z1)^2-er^2, f=x/(x+y), iter2=iter+f)

### Histogram coloring ###
h <- m2 %>%
    filter(!is.na(iter2)) %>% ## NA を除外
    count(iter2) %>%          ## iterの値毎に個数をカウント
    mutate(n=cumsum(n)) %>%   ## 累積和を求める
    mutate(n=n/max(n))        ## 値域を[0,1]にする

m2 <- m2 %>%
    left_join(h, by='iter2') %>%
    mutate(n=if_else(is.na(iter2), NA_real_, n))

ggplot(m2) +
    aes(x=Re(c), y=Im(c), fill=n) +
    geom_raster() +
    scale_fill_gradientn(colors=rainbow(13), na.value='black') +
    coord_fixed(ratio=1)

参考文献

• 滑らか coloring - マンデルブロ集合の外側を滑らかに彩色する
• Histogram coloring - Plotting algorithms for the Mandelbrot set