# source("http://stat-help.net/n1b.txt") # f'n is f(z)=z^2-2^z get.n1b=function(nx, ny=nx, xlimss=c(-3, 10), ylimss=c(-10, 10)){ x.coords=seq(xlimss[1], xlimss[2], length=nx) y.coords=seq(ylimss[1], ylimss[2], length=ny) counts=rep(0, n<-nx*ny) log2=log(2) z=outer(x.coords,y.coords, FUN=function(x,y) complex(real=x, imag=y)) dist=3 count=0 tmat=is.finite(z) cc=z countmat=matrix(1, nrow=nrow(z), ncol=ncol(z)) for(count in 1:22){ f.z=z*z-2^z der.f.z=2*z-log2*2^z new.z=z-f.z/der.f.z tmat[tmat]=is.finite(new.z[tmat]) & Mod(z[tmat]-new.z[tmat])>0.1 countmat[tmat]=countmat[tmat]+1 z=new.z } x.coords=c(t(Re(cc))); y.coords=c(t(Im(cc))) counts=c(t(countmat)) return(list(map=cbind(x.coords, y.coords, counts), xlimss=xlimss, ylimss=ylimss)) } newton1b.res=get.n1b(480) # colorss=topo.colors(25); # colorss=heat.colors(25); new.newton1b=newton1b.res$map; colorss=rainbow(length(unique(new.newton1b[,3])), start=1/3, end=3/4) colorss=colorss[sample(length(colorss))] # par(bg="grey") # jpeg("Newton1b.jpeg") par(mar=c(0,0,0,0)) xlimss=newton1b.res$xlimss ylimss=newton1b.res$ylimss plot(new.newton1b[, 1:2], pch=".", col=colorss[new.newton1b[,3]], xlab=" ",ylab=" ", cex.lab=1.5, cex.main=1.3, main=" ", xlim=c(xlimss[1]+0.2, xlimss[2]-0.2), ylim=c(ylimss[1]+0.2, ylimss[2]-0.2)) # dev.off()