反復関数系でフラクタル

シェルピンスキーのガスケット

gasket <- function(n=100000) {
m <- matrix(c(0.5, 0, 0, 0.5) , 2)
lst <- matrix(0, n, 2) # matrixを用意
lst[1,] <- c(100, 0)   # 初期値

prb <- sample(1:3, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i], # switchで振り分け
                    m %*% lst[i-1,],            # 1/3
                    m %*% lst[i-1,] + c(100,0), # 1/3
                    m %*% lst[i-1,] + c(50, 30))# 1/3
}
  plot(lst, pch=".", axes=F, asp=sqrt(3)*1.5)
}

ドラゴンカーブ

dragon <- function(n=100000){
m <- matrix(c(0.5, -0.5, 0.5, 0.5) ,2)
lst <- matrix(0, n, 2)        # matrixを用意
lst[1,] <- c(1, 1) # 初期値

prb <- sample(1:2, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m %*% lst[i-1,] + c( 0.125, 0.625), # 1/2
                    m %*% lst[i-1,] + c(-0.125, 0.375)) # 1/2
}
  plot(lst, pch=".", asp=1,axes=F)
}

コッホ曲線

koho <- function(n=100000){
m1 <- matrix(c(0.5, sqrt(3)/6, sqrt(3)/6, -0.5) ,2)
m2 <- matrix(c(0.5, -sqrt(3)/6, -sqrt(3)/6, -0.5) ,2)
lst <- matrix(0, n, 2)
lst[1,] <- c(1, 1)

prb <- sample(1:2, n,rep=TRUE)
for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,],
                    m2 %*% lst[i-1,] + c(0.5,sqrt(3)/6))
}
  plot(lst, pch=".", axes=F,asp=0.86)
}

コッホ雪片

• http://ecademy.agnesscott.edu/~lriddle/ifs/ksnow/IFSdetailsHexagons.htm

snow <- function(n=100000) {
m1 <- matrix(c(1/2,sqrt(3)/6, -sqrt(3)/6,1/2),2)
m2 <- diag(1/3,2)

lst <- matrix(0, n, 2)      # matrixを用意

prb <- sample(1:7, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,],                     # 1/7
                    m2 %*% lst[i-1,] + c(1/sqrt(3),1/3),  # 1/7
                    m2 %*% lst[i-1,] + c(0,2/3),          # 1/7
                    m2 %*% lst[i-1,] + c(-1/sqrt(3),1/3), # 1/7
                    m2 %*% lst[i-1,] + c(-1/sqrt(3),-1/3),# 1/7
                    m2 %*% lst[i-1,] + c(0,-2/3),         # 1/7
                    m2 %*% lst[i-1,] + c(1/sqrt(3),-1/3)) # 1/7
}
  plot(lst, pch=".", axes=F,asp=1)
}

Cカーブ

cc <- function(n=100000){
m1 <- matrix(c(0.5, 0.5, -0.5, 0.5) ,2)
m2 <- matrix(c(0.5, -0.5, 0.5, 0.5) ,2)
lst <- matrix(0, n, 2)
lst[1,] <- c(1, 1)

prb <- sample(1:2, n,rep=TRUE)
for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,],
                    m2 %*% lst[i-1,] + c(-0.5, 0.5))
}
  plot(lst, pch=".", asp=1,axes=F)
}

楓の葉

kaede <- function(n=100000){
m1 <- matrix(c(0.8, 0, 0, 0.8) ,2)
m2 <- matrix(c(0.5, 0, 0, 0.5) ,2)
m3 <- matrix(c(0.355, 0.355, -0.355, 0.355) ,2)
m4 <- matrix(c(0.355, -0.355, 0.355, 0.355) ,2)

lst <- matrix(0, n, 2)
lst[1,] <- c(1, 1) # 初期値

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.5,0.168, 0.166, 0.166)), labels=1:4))

for(i in 2:n) {
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,] + c(0.1, 0.04),     # 50.0%
                    m2 %*% lst[i-1,] + c(0.25,0.4),      # 16.8%
                    m3 %*% lst[i-1,] + c(0.266, 0.078),  # 16.6%
                    m4 %*% lst[i-1,] + c(0.378, 0.434))  # 16.6%
}
  plot(lst, pch=".", asp=1,axes=F)
}

羊歯

sida <- function(n=100000){
m1 <- matrix(c( 0.856,-0.0205, 0.0414,0.858) ,2)
m2 <- matrix(c( 0.244, 0.176, -0.385, 0.224) ,2)
m3 <- matrix(c(-0.144, 0.181,  0.39,  0.259) ,2)
m4 <- matrix(c( 0,     0.355,  0,     0.216) ,2)

lst <- matrix(0, n, 2)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.73,0.13, 0.13, 0.01)), labels=1:4))

for(i in 2:n) {
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,] + c(0.07,  0.147), # 73.0%
                    m2 %*% lst[i-1,] + c(0.393,-0.102), # 13.0%
                    m3 %*% lst[i-1,] + c(0.527,-0.014), # 13.0%
                    m4 %*% lst[i-1,] + c(0.486, 0.05))  #  1.0%
}
  plot(lst, pch=".", asp=1,axes=F)
}

アンモナイト

ammonite <- function(n=100000) {

m1 <- matrix(c(-0.29, 0,   0   , 0.2 ), 2)
m2 <- matrix(c(-0.07, 0.01,0.02, 0.29), 2)
m3 <- matrix(c( 0.94,0.21,-0.22, 0.96), 2)

lst <- matrix(0, n, 2)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.06,0.02,0.92)), labels=1:3))

for(i in 2:n){
   lst[i,]<-switch(prb[i],
                   m1 %*% lst[i-1,] + c( 0.59,-0.32), #  6.0%
                   m3 %*% lst[i-1,] + c( 0.79,-0.06), #  2.0%
                   m3 %*% lst[i-1,] + c(-0.05, 0.01)) # 92.0%
}
  plot(lst, pch=".", asp=1, xlim=c(-1,1), ylim=c(-1,1), axes=F)
}

シェルピンスキーのペンタゴン

• http://ecademy.agnesscott.edu/~lriddle/ifs/pentagon/pentagon.htm

pentagon <- function(n=100000) {
m <- diag(0.382,2)
lst <- matrix(0, n, 2)
lst[1,] <- c(0, 0)

prb <- sample(1:5, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m %*% lst[i-1,],
                    m %*% lst[i-1,] + c( 0.618, 0),
                    m %*% lst[i-1,] + c( 0.809, 0.588),
                    m %*% lst[i-1,] + c( 0.309, 0.951),
                    m %*% lst[i-1,] + c(-0.191, 0.588))
}
  col <- c("#0085C7","#F4C300","#000000","#009F3E","#DF0024") # オリンピックカラー
  plot(lst, pch=".", axes=F,asp=1,col=col[prb])
}

シェルピンスキーのカーペット

• http://ecademy.agnesscott.edu/~lriddle/ifs/carpet/carpet.htm

carpet <- function(n=100000) {
m <- diag(1/3,2)
lst <- matrix(0, n, 2)
lst[1,] <- c(0, 0)

prb <- sample(1:8, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m %*% lst[i-1,],
                    m %*% lst[i-1,] + c(0  ,1/3),
                    m %*% lst[i-1,] + c(0  ,2/3),
                    m %*% lst[i-1,] + c(1/3,0),
                    m %*% lst[i-1,] + c(1/3,2/3),
                    m %*% lst[i-1,] + c(2/3,0),
                    m %*% lst[i-1,] + c(2/3,1/3),
                    m %*% lst[i-1,] + c(2/3,2/3))
}

  col = sample(c("#0068B7", "#00693E", "#008DCB", "#009E96", "#00A051", "#00A0E9", "#187FC4", "#1D2088", "#86B81B", "#920783", "#9FD9F6", "#D3DEF1", "#D4ECEA", "#E4007F", "#EA5504", "#EA5532", "#ED6C00", "#F39800", "#FFF100"),8)
  plot(lst, pch=".", axes=F,asp=1,col=col[prb])
}

McWorter’s Pentigree IFS

• http://ecademy.agnesscott.edu/~lriddle/ifs/pentigre/pentigre2ndFormIFSdetails.htm

pentigree <- function(n=100000) {
m <- matrix(c(0.309,0.255,-0.255,0.309),2)

lst <- matrix(0, n, 2)

prb <- sample(1:6, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m %*% lst[i-1,],                   # 1
                    m %*% lst[i-1,]+c( 0.727, 0),      # 2
                    m %*% lst[i-1,]+c( 0.225, 0.691),  # 3
                    m %*% lst[i-1,]+c(-0.588, 0.427),  # 4
                    m %*% lst[i-1,]+c(-0.588,-0.427),  # 5
                    m %*% lst[i-1,]+c( 0.255,-0.691))  # 6
}
  plot(lst, pch=".", axes=F,asp=1,col=prb)
}

golden dragon

• http://ecademy.agnesscott.edu/~lriddle/ifs/heighway/goldenDragon.htm

gdragon <- function(n=100000) {
m  <- matrix(c( 0.62367,0.40337,-0.40337, 0.62367),2)
m1 <- matrix(c(-0.37633,0.40337,-0.40337,-0.37633),2)

lst <- matrix(0, n, 2)

prb <- sample(1:2, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m  %*% lst[i-1,],         # 1
                    m1 %*% lst[i-1,]+c(1, 0)) # 2
}
  plot(lst, pch=".", axes=F,asp=1,col=prb)
}

penta

• http://ecademy.agnesscott.edu/~lriddle/ifs/pentaden/penta.htm

penta <- function(n=100000) {
m136 <- matrix(c( 0.341, 0.071,-0.071, 0.341),2)
m2   <- matrix(c( 0.038, 0.346,-0.346, 0.038),2)
m4   <- matrix(c(-0.234,-0.258, 0.258,-0.234),2)
m5   <- matrix(c( 0.173,-0.302, 0.302, 0.173),2)

lst <- matrix(0, n, 2)

prb <- sample(1:6, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m136  %*% lst[i-1,],                 # 1
                    m2    %*% lst[i-1,]+c(0.341, 0.071), # 2
                    m136  %*% lst[i-1,]+c(0.379, 0.418), # 3
                    m4    %*% lst[i-1,]+c(0.720, 0.489), # 4
                    m5    %*% lst[i-1,]+c(0.486, 0.231), # 5
                    m136  %*% lst[i-1,]+c(0.659,-0.071)) # 6
}
  col <- sample(c("#E60012", "#EB6100", "#F39800", "#FCC800", "#FFF100", "#CFDB00", "#8FC31F", "#22AC38", "#009944", "#009B6B", "#009E96", "#00A0C1", "#00A0E9", "#0086D1", "#0068B7", "#00479D", "#1D2088", "#601986", "#920783", "#BE0081", "#E4007F", "#E5006A", "#E5004F", "#E60033"),6)
  plot(lst, pch=".", axes=F,asp=1,col=col[prb])
}

クリスタル

crystal <- function(n=100000) {
m1 <- matrix(c(0,.5,-.5,0),2)
m2 <- matrix(c(0,-.5,.5,0),2)
m3 <- matrix(c(.5,0,0,.5),2)

lst <- matrix(0, n, 2)

prb <- sample(1:3, n,rep=TRUE)

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,]+c(0.5,0),    # 1
                    m2 %*% lst[i-1,]+c(0.5,0.5),  # 2
                    m3 %*% lst[i-1,]+c(0.25,0.5)) # 3
}
col <- c("#00A0E9", "#E60012", "#1E2C5C", "#910000") # ４色
  plot(lst, pch=".", axes=F,asp=1,col=col[prb])
}

coral

coral <- function(n=100000) {

m1 <- matrix(c(-0.16666667, 0.16666667,-0.1666667,-0.1666667),2)
m2 <- matrix(c( 0.83333333,-0.25000000, 0.2500000, 0.8333333),2)
m3 <- matrix(c( 0.33333333, 0.08333333,-0.0833333, 0.3333333),2)

lst <- matrix(0, n, 2)
lst[1,] <- c(1,1)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.163,0.6,0.237)), labels=1:3))

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,],
                    m2 %*% lst[i-1,] + c(-0.1666667,-0.166667),
                    m3 %*% lst[i-1,] + c( 0.0833333, 0.666667))
}
  plot(lst, pch=".", asp=1, axes=F)
}

スパイラル

• http://cs.lmu.edu/~ray/notes/ifs/

spiral <- function(n=100000) {

m1 <- matrix(c(0.787879,0.242424,-0.424242,0.859848),2)
m2 <- matrix(c(-0.121212,0.090909,0.257576,0.053030),2)
m3 <- matrix(c(0.252525,0.252525,-0.136364,0.181818),2)

lst <- matrix(0, n, 2)
lst[1,] <- c(1,1)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.895652,0.052174,0.052174)), labels=1:3))

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,] + c(1.758647,1.408065),
                    m2 %*% lst[i-1,] + c(-3.721654,1.377236),
                    m3 %*% lst[i-1,] + c(3.086107,1.568035))
}
  plot(lst, pch=".", asp=1, axes=F)
}

スティック

• http://cs.lmu.edu/~ray/notes/ifs/

sticks <- function(n=100000) {
m1 <- matrix(c(0.005,0.000,0.000,0.500),2)
m2 <- matrix(c(0.414,0.414,-0.414,0.414),2)
m3 <- matrix(c(0.414,-0.414,0.414,0.414),2)

lst <- matrix(0, n, 2)
lst[1,] <- c(1,1)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.12,0.44,0.44)), labels=1:3))

for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,] + c(0.0,0.0),
                    m2 %*% lst[i-1,] + c(0.0,0.5),
                    m3 %*% lst[i-1,] + c(0.0,0.5))
}
  plot(lst, pch=".", asp=1, axes=F)
}

ツリー

• http://cs.lmu.edu/~ray/notes/ifs/

tree <- function(n=100000) {
m1 <- matrix(c(0.000,0.000,0.000,0.600),2)
m2 <- matrix(c(0.440,0.000,0.000,0.550),2)
m3 <- matrix(c(0.343,0.199,-0.248,0.429),2)
m4 <- matrix(c(0.343,-0.199,0.248,0.429),2)
m5 <- matrix(c(0.280,0.280,-0.350,0.350),2)
m6 <- matrix(c(0.280,-0.280,0.350,0.350),2)

lst <- matrix(0, n, 2)
lst[1,] <- c(1,1)

prb <- as.numeric(cut(runif(n), cumsum(c(0,0.1,0.18,0.18,0.18,0.18,0.18)), labels=1:6))


for(i in 2:n){
  lst[i,] <- switch(prb[i],
                    m1 %*% lst[i-1,] + c(0.00,-0.065),
                    m2 %*% lst[i-1,] + c(0.00,0.200),
                    m3 %*% lst[i-1,] + c(-0.03,0.100),
                    m4 %*% lst[i-1,] + c(0.03,0.100),
                    m5 %*% lst[i-1,] + c(-0.05,0.000),
                    m6 %*% lst[i-1,] + c(0.05,0.000))
}
  col <- c("#89C997", "#69BD83", "#3EB370", "#00A95F", "#00A051", "#009944")
  plot(lst, pch=".", asp=1, axes=F,xlim=c(-0.2,0.2),ylim=c(-0.2,0.5),col=col[prb])
}

一度に描写

zenbu<-function(){
par(mai=c(0,0,0,0), mfrow=c(3,6))
  koho(10000)
  cc(10000)
  gasket(10000)
  kaede(10000)
  ammonite(10000)
  dragon(10000)
  sida(10000)
  pentagon(10000)
  carpet(10000)
  snow(10000)
  pentigree(10000)
  gdragon(10000)
  penta(10000)
  crystal(10000)
  coral(10000)
  spiral(10000)
  sticks(10000)
  tree(10000)
}