用R来玩几个经典的分形图,很好玩。感兴趣的运行一下,跟帖发个图片上来,多谢支持,谢谢:)
第一个:
<br />
plot.tri <- function(n = 1000, col ="blue", ani=FALSE, cex=1.2){<br />
p <- runif(n); <br />
X <- rbind(rep(0, n), rep(0, n))<br />
B <- cbind(c(0,0),c(0.25,0.433),c(0.5,0))<br />
if(ani) plot(0,0,xlim=c(0,1),ylim=c(0,0.85),type="n",xlab="",ylab="")<br />
for(i in 2:n){ <br />
pp <- p[i]; <br />
ind <- rank(c(c(1/3,2/3,1), pp), ties.method="min")[4]<br />
X[,i] <- 0.5*X[,i-1] + B[,ind]<br />
if(ani) points(X[1,i], X[2,i],pch = ".", cex = 1, col = col)<br />
}<br />
if(!ani) plot(X[1,], X[2,],pch = ".", cex = cex, col = col, xlab="", ylab="") <br />
}<br />
<br />
### example<br />
plot.tri(100000, ani=TRUE)<br />
第二个:
<br />
plot.koch <- function(k,col="blue"){ <br />
plot(0,0,xlim=c(0,1), ylim=c(-sqrt(3)/6,sqrt(3)/2), asp = 1,type="n",xlab="", ylab="")<br />
plotkoch <- function(x1,y1,x2,y2,n){<br />
if (n > 1){<br />
plotkoch(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1); <br />
plotkoch((2*x1+x2)/3,(2*y1+y2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1) *sqrt(3)/6,n-1);<br />
plotkoch((x1+x2)/2-(y1-y2)*sqrt(3)/6,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*x2+x1)/3,(2 *y2+y1)/3,n-1); <br />
plotkoch((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1)<br />
} else { <br />
x=c(x1,(2*x1+x2)/3,(x1+x2)/2-(y1-y2)*sqrt(3)/6,(2*x2+x1)/3,x2); <br />
y=c(y1,(2*y1+y2)/3,(y1+y2)/2-(x2-x1)*sqrt(3)/6,(2*y2+y1)/3,y2); <br />
lines(x,y,type="l",col=col) <br />
}<br />
}<br />
plotkoch(0,0,1,0,k) <br />
plotkoch(0.5,sqrt(3)/2,0,0,k) <br />
plotkoch(1,0,0.5,sqrt(3)/2,k)<br />
}<br />
<br />
## example<br />
for(i in 1:5){<br />
plot.koch(i,col=i) <br />
Sys.sleep(1)<br />
}<br />
第三个:
<br />
plot.leaf <- function(n=100000, col="green",cex=2){<br />
x <- c(.5, .5);<br />
plot(x[1], x[2], xlim=c(-3, 3), ylim =c(0, 10),type="n")<br />
p <- c( .85, .92, .99, 1.00);<br />
A <- rbind(c(.85, .04), c(-.04,.85), c(.20,-.26), c(.23,.22),<br />
c(-.15,.28), c(.26,.24), c(0, 0), c(0, .16))<br />
B <- cbind(c(0, 1.6), c(0, 1.6), c(0,.44), c(0,0))<br />
<br />
for (i in 1:n){<br />
ran <- runif(1);<br />
ind <- rank(c(p, ran), ties.method="min")[5]<br />
x <- A[(2*ind-1):(2*ind),]%*%x + B[,ind]<br />
points(x[1],x[2], pch=".", cex=cex, col=col)<br />
}<br />
} <br />
plot.leaf(cex=1.6)<br />
第四个:
<br />
plot.tree <- function(x1,y1,x2,y2,n,xlim=c(-1,1), ylim=c(0,2), col="blue", add=FALSE){<br />
if(!add)<br />
plot(0,0,xlim=xlim, ylim=ylim, type="n",xlab="", ylab="",asp=1)<br />
tree <- function(x1,y1,x2,y2,n){<br />
flag <- 0;<br />
theta <- pi/6;<br />
if (x2<x1) flag <- 1;<br />
if (n>1){<br />
tree(x1,y1,(2*x1+x2)/3,(2*y1+y2)/3,n-1);<br />
tree((2*x1+x2)/3,(2*y1+y2)/3,(2*x2+x1)/3,(2*y2+y1)/3,n-1);<br />
tree((2*x2+x1)/3,(2*y2+y1)/3,x2,y2,n-1);<br />
tree((2*x1+x2)/3,(2*y1+y2)/3,(2*x1+x2)/3+sin(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),(2*y1+y2)/3+cos(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),n-1);<br />
tree((2*x2+x1)/3,(2*y2+y1)/3,(2*x2+x1)/3+sin(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),(2*y2+y1)/3+cos(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3),n-1);<br />
} else {<br />
x <- c(x1,x2);<br />
y <- c(y1,y2);<br />
xx <- c((2*x1+x2)/3,(2*x1+x2)/3+sin(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));<br />
yy <- c((2*y1+y2)/3,(2*y1+y2)/3+cos(pi/2-atan((y2-y1)/(x2-x1))-theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));<br />
xxx <- c((2*x2+x1)/3,(2*x2+x1)/3+sin(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));<br />
yyy <- c((2*y2+y1)/3,(2*y2+y1)/3+cos(pi/2-atan((y2-y1)/(x2-x1))+theta+flag*pi)*sqrt(((y2-y1)^2+(x2-x1)^2)/3));<br />
lines(x, y, type="l",col=col)<br />
lines(xx, yy, type="l",col=col)<br />
lines(xxx, yyy, type="l",col=col)<br />
}<br />
}<br />
tree(x1,y1,x2,y2,n)<br />
}<br />
<br />
## example<br />
for(i in 1:5){<br />
plot.tree(0,0,0,1.5,i,col=i)<br />
Sys.sleep(1)<br />
}<br />