今回は『バーニングシップフラクタル』についての備忘録です。
バーニングシップフラクタルとは？
Wikipedia - Bruning Ship Fractal より抜粋
バーニングシップフラクタルは、1992年にマイケル・ミケリッチとオットー・E・レスラーによって初めて記述・生成されたもので、以下の関数を反復することで生成される：
\[ z_{n+1} = (|Re(z_n)| + i|Im(z_n)|)^2 + c,\ z_0=0 \]
マンデルブロ集合との違い
下記はマンデルブロ集合の定義式です。前項の定義式と見比べると二乗の中の式が異なる以外は一緒です。描画方法も同じ。
\[ z_{n+1} = {z_n}^2 + c,\ z_0=0 \]
プログラム
下記はバーニングシップフラクタルを生成するR言語用プログラムです。 以前に紹介した マンデルブロ集合用のプログラム を少し書き換えただけです。
まずは、関数の定義。
bsf <- function(
center =-0.25-0.5i, ## 計算する範囲の中心
radius = 2.0, ## 計算する範囲の半径
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*1i, '+')), z=0, iter=0)
for (i in 1:max.iteration) {
tmp <- Mod(m$z) < escape.radius
m <- m %>%
mutate(
z=if_else(tmp, (abs(Re(z)) + 1i*abs(Im(z)))^2 + c, z), ### 式変更 ###
iter=if_else(tmp, iter+1, iter)
)
}
m <- m %>%
mutate(iter=na_if(iter, max.iteration)) ## 反復回数上限に達した部分を NA に置き換え
return(m)
}次に計算を実行します。
描画方法はマンデルブロ集合と同じという事で、Continuous (smooth) coloring や Rank-Order Mapping も適用します。
library(tidyverse)
library(viridis)
b1 <- bsf() %>%
### Continuous (smooth) coloring ###
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(iter=iter-log10(log10(Mod(z)))/log10(2)) %>%
### Rank-Order Mapping ###
mutate(ci=percent_rank(iter))
b1 %>%
ggplot() +
aes(x=Re(c), y=Im(c), fill=ci) +
geom_raster() +
scale_fill_viridis(option='magma', na.value='black') +
guides(fill='none') +
coord_fixed(ratio=1)
想像していた図と違う…上下逆さま。
調べてみると 縦軸反転が慣例 のようです。 以降の図は縦軸を反転します。
b1 %>%
ggplot() +
aes(x=Re(c), y=Im(c), fill=ci) +
geom_raster() +
scale_fill_viridis(option='magma', na.value='black') +
guides(fill='none') +
scale_y_reverse() + ### 縦軸反転 ##
coord_fixed(ratio=1)
左側の尻尾部分を拡大
bsf(center =-1.755-0.04i,
radius = 0.055,
max.iteration = 500,
escape.radius = 5) %>%
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(iter=iter-log10(log10(Mod(z)))/log10(2)) %>%
mutate(ci=percent_rank(iter)) %>%
ggplot() +
aes(x=Re(c), y=Im(c), fill=ci) +
geom_raster() +
scale_fill_viridis(option='magma', na.value='black') +
guides(fill='none') +
scale_y_reverse() +
coord_fixed(ratio=1)
さらに左側
bsf(
center =-1.861-0.0025i,
radius = 0.004,
max.iteration = 500,
escape.radius = 5 ) %>%
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(z=(abs(Re(z)) + 1i*abs(Im(z)))^2+c, iter=iter+1) %>%
mutate(iter=iter-log10(log10(Mod(z)))/log10(2)) %>%
mutate(ci=percent_rank(iter)) %>%
ggplot() +
aes(x=Re(c), y=Im(c), fill=ci) +
geom_raster() +
scale_fill_viridis(option='magma', na.value='black') +
guides(fill='none') +
scale_y_reverse() +
coord_fixed(ratio=1)
すてき。
