反復関数系でフラクタル

@manozo

2015年9月27日

反復関数系(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)
}

コッホ曲線

コッホ雪片

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

羊歯

アンモナイト

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

アンモナイト

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

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

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

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

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

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