I am listening to the audiobook Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Lifeby Steven Strogatz which I got from Audible. Obviously a mathematical book is not ideal to listen to, but lacking illustrations I can make them myself. On top of that, you learn more from programming than listening.

He mentioned a simple system which synchronizes itself. As far as I understand a large number of items are considered. Each  item has more or less the same properties. Small steps are taken. Every step has a small additive effect, a multiplicative decrease. When an item gets over a limit, then it resets to zero and causes a small increase in the other items.

What is supposed to happen is that at low values you get steep slopes and at high values small slopes. As time progresses the small increases due to resets cause a synchronization.

To program this was not difficult, except for the initialization. If you start the items at random levels you have too many at low levels. In the end I decided to add an item at each time point for an init period. What did surprise me, was that calculating what happens actually takes less computer time than plotting it. But maybe that is because I wanted to use ggplot and a low alpha to show the bunching of items. I also discovered it is also not wise to display too complex a figure in ggplot running in Eclipse. It can freeze the whole program. I took to saving all plots rather than displaying within Eclipse.

In the init phase every item has different start and ending times. Note the black line at level 0 which consists of not yet running items

After 10000 steps it is still chaos, with a few somewhat fatter lines.

Not that they are perfect after 300000 steps, but it gets more and more closer.

R code

library(ggplot2)

plotpart <- function(xmat,fname='een.png',xstart=1) {
  fin <- ncol(xmat)
  nitems <- nrow(xmat)
  df <- data.frame(score=as.numeric(xmat),
      time=rep((1:fin)+xstart,each=nitems),
      item =rep(1:nitems,fin))
  g<- ggplot(df,aes(x=time,y=score,group=item,alpha=score)) + 
      geom_line() + 
      scale_alpha_continuous(range=c(.05,.0501) )+
      theme(legend.position=’none’)
  png(fname)
  print(g)
  dev.off()
}

nextgen <- function(x,add,shrink=.99,limit=1,spil=.0105) {
  x <- x*shrink+add
  maxed <- x>limit
  x[maxed]<- .01#add[maxed]
  x[!maxed] <- x[!maxed] + sum(maxed)*spil
  x[x>limit] <- limit
  x
}

ngenmat <- function(xstart,niter,add,spil) {
  nitems <- length(xstart)
  xmat <- matrix(0,nrow=nitems,ncol=niter)
  xmat[,1 ] <- xstart
  for ( i in 2:niter) {
    xmat[,i] <- nextgen(xmat[,i-1],add=add,spil=spil)
  }
  xmat
}

ngenmove <- function(xstart,niter,add,spil) {
  for ( i in 1:niter) {
    xstart <- nextgen(xstart,add=add,spil=spil)
  }
  xstart
}

# init a run
nitems <- 1000
niter <- 1000
xmat <- matrix(0,nrow=nitems,ncol=niter)
for (i in 1:nitems) xmat[i,i] <- runif(1,min=.01,max=.02)
add <- rnorm(nitems,.0103,.00005)
spil=.000025
for ( i in 2:niter) {
  todo <- xmat[,i-1]>0
  xmat[todo,i] <- nextgen(xmat[todo,i-1],add=add[todo],spil=spil)
}
#and show it
show <- sample(1:nrow(xmat),50)
plotpart(xmat[show,],fname=’een.png’,xstart=1)

# move 9000 steps
x1 <- ngenmove(xmat[,ncol(xmat)],niter=9000,add=add,spil=spil)
xm1 <- ngenmat(x1,niter=1000,add=add,spil=spil)
plotpart(xm1[show,],fname=’twee.png’,xstart=10000)

x2 <- ngenmove(xm1[,ncol(xm1)],niter=39000,add=add,spil=spil)
xm2 <- ngenmat(x2,niter=1000,add=add,spil=spil)
plotpart(xm2[show,],fname=’drie.png’,xstart=50000)

x3 <- ngenmove(xm2[,ncol(xm2)],niter=49000,add=add,spil=spil)
xm3 <- ngenmat(x3,niter=1000,add=add,spil=spil)
plotpart(xm3[show,],99000,fname=’vier.png’)

x4 <- ngenmove(xm3[,ncol(xm3)],niter=49000,add=add,spil=spil)
xm4 <- ngenmat(x4,niter=1000,add=add,spil=spil)
plotpart(xm4[show,],150000,fname=’vijf.png’)

x5 <- ngenmove(xm4[,ncol(xm4)],niter=49000,add=add,spil=spil)
xm5 <- ngenmat(x5,niter=1000,add=add,spil=spil)
plotpart(xm5[show,],200000,fname=’zes.png’)

x6 <- ngenmove(xm5[,ncol(xm5)],niter=99000,add=add,spil=spil)
xm6 <- ngenmat(x6,niter=1000,add=add,spil=spil)
plotpart(xm6[show,],300000,fname=’zeven.png’)