反復関数系(IFS)を使ってフラクタルを描いてみた。

switchを使って分岐。
スペシャルサンクス：「カオスとフラクタル◎Excelで体験」臼田昭司、東野勝治、井上祥史、伊藤敏、葭谷安正共著（オーム社出版局）1999

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

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,],
                      m %*% lst[i-1,] + c(100,0),
                      m %*% lst[i-1,] + c(50, 30))
  }
  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),
                      m %*% lst[i-1,] + c(-0.125, 0.375))
  }
  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)        # matrixを用意
  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,],
                      m2 %*% lst[i-1,] + c(1/sqrt(3),1/3),
                      m2 %*% lst[i-1,] + c(0,2/3),
                      m2 %*% lst[i-1,] + c(-1/sqrt(3),1/3),
                      m2 %*% lst[i-1,] + c(-1/sqrt(3),-1/3),
                      m2 %*% lst[i-1,] + c(0,-2/3),
                      m2 %*% lst[i-1,] + c(1/sqrt(3),-1/3))
  }
  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)        # matrixを用意
  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)
}

Cカーブ

楓の葉

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(factor(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),
                      m2 %*% lst[i-1,] + c(0.25,0.4),
                      m3 %*% lst[i-1,] + c(0.266, 0.078),
                      m4 %*% lst[i-1,] + c(0.378, 0.434))
  }
  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(factor(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),
                      m2 %*% lst[i-1,] + c(0.393,-0.102),
                      m3 %*% lst[i-1,] + c(0.527,-0.014),
                      m4 %*% lst[i-1,] + c(0.486, 0.05))
  }
  plot(lst, pch=".", asp=1,axes=F)
}

羊歯

アンモナイト

• http://www.env.kochi-tech.ac.jp/monma/vba_seminar/pages/quizzes/ammonite.html

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)

p1 <- 0.06
p2 <- 0.02
p3 <- 1-p1-p2

prb <- factor(cut(runif(n), cumsum(c(0,p1,p2,p3)), labels=1:3))

for(i in 2:n){
lst[i,]<-switch(as.numeric(prb[i]),
m1 %*% lst[i-1,] + c( 0.59,-0.32),
m3 %*% lst[i-1,] + c( 0.79,-0.06),
m3 %*% lst[i-1,] + c(-0.05, 0.01))
}
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)        # matrixを用意
  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)        # matrixを用意
  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)      # matrixを用意
  
  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)
}

McWorter’s Pentigree IFS

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)      # matrixを用意
  
  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)
}

golden dragon

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)      # matrixを用意
  
  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])
}

penta

クリスタル

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)      # matrixを用意
  
  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)
  p1 <- 0.163
  p2 <- 0.600
  p3 <- 1-p1-p2
  
  prb <- as.numeric(factor(cut(runif(n), cumsum(c(0,p1,p2,p3)), 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)
}

coral

スパイラル

• 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)
  p1<-0.895652
  p2<-0.052174
  p3 <- 1-p1-p2
  
  prb <- as.numeric(factor(cut(runif(n), cumsum(c(0,p1,p2,p3)), 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)
  p1<-0.12
  p2<-0.44
  p3 <- 1-p1-p2
  
  prb <- as.numeric(factor(cut(runif(n), cumsum(c(0,p1,p2,p3)), 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)
  p1<-0.100
  p2<-0.180
  p3<-0.180
  p4<-0.180
  p5<-0.180
  p6 <- 1-p1-p2-p3-p4-p5
  
  prb <- as.numeric(factor(cut(runif(n), cumsum(c(0,p1,p2,p3,p4,p5,p6)), 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(n=10000){
  par(mai=c(0,0,0,0), mfrow=c(3,6))
  koho(n)
  cc(n)
  gasket(n)
  kaede(n)
  ammonite(n)
  dragon(n)
  sida(n)
  pentagon(n)
  carpet(n)
  snow(n)
  pentigree(n)
  gdragon(n)
  penta(n)
  crystal(n)
  coral(n)
  spiral(n)
  sticks(n)
  tree(n)
}
zenbu()

一度に描写