Below, you find the R code for computing the Game of life n=20 # Size of matrix mat=matrix(0,ncol=n,nrow=n) # Create a n x n matrix with zeros mat[5:14,10]=1 # Add 10 live cells temp_mat=mat # Create a temporary matrix image(t(apply(mat, 2, rev)),col=c("grey50","seagreen1"),yaxt="n",xaxt="n") # Plot image grid(nx=n,ny=n,col="grey70",lty=1) mat1=mat[c(2:n,1),];mat2=mat[c(n,1:n-1),]; mat3=mat[,c(2:n,1)]; mat4=mat[,c(n,1:n-1)]; mat5=mat[c(2:n,1),c(2:n,1)]; mat6=mat[c(2:n,1),c(n,1:n-1)] mat7=mat[c(n,1:n-1),c(n,1:n-1)];mat8=mat[c(n,1:n-1),c(2:n,1)]; for (k in 1:200){ # Repeat 200 times numb_alive=mat1+mat2+mat3+mat4+mat5+mat6+mat7+mat8 temp_mat[mat==1 & numb_alive<2]=0 temp_mat[mat==1 & numb_alive>3]=0 temp_mat[mat==1 & (numb_alive==2 | numb_alive==3)]=1 temp_mat[mat==0 & numb_alive==3]=1 mat=temp_mat # Update matrix mat1=mat[c(2:n,1),];mat2=mat[c(n,1:n-1),];mat3=mat[,c(2:n,1)];mat4=mat[,c(n,1:n-1)] mat5=mat[c(2:n,1),c(2:n,1)];mat6=mat[c(2:n,1),c(n,1:n-1)];mat7=mat[c(n,1:n-1),c(2:n,1)] mat8=mat[c(n,1:n-1),c(n,1:n-1)] image(t(apply(mat, 2, rev)),col=c("grey50","seagreen1"),add=TRUE) # Plot image grid(nx=n,ny=n,col="grey70",lty=1) Sys.sleep(0.5) # To see changes on the screen we need to pause the loop }
Below, you find the R code for the SIR model rm(list=ls()) q=0.5 # Probability for recovery n=50 # Size of matrix 50 x 50 mat=matrix(ncol=n,nrow=n,0) # Create matrix with zeros (Healthy individuals) mat[25,15]=1 # Place an infected individual at pos (25,15) image(mat,col=c("grey50","deeppink","seagreen1"),yaxt="n",xaxt="n",zlim=c(0,2)) # Plot image grid(nx=n,ny=n,col="grey70",lty=1) temp=mat t=80 # Number of time steps n_helalthy=rep(0,t) # Vector to count healthy at each time step n_infected=rep(0,t) # Vector to count infected at each time step n_resistant=rep(0,t) # Vector to count resistant at each time step for (k in 1:t){ # Repeat t times kH=0 # Initialize counter for healthy kI=0 # Initialize counter for infected kR=0 # Initialize counter for resistant # Step through each element in the matrix for (i in 1:n){ for (j in 1:n){ if (mat[i,j]==0) { kH=kH+1} # Count healthy if (mat[i,j]==1) { kI=kI+1} # Count infected if (mat[i,j]==2) { kR=kR+1} # Count resistant R=0 # Initialize counter for number of infected neighbors if (mat[i,j]==0){ # If healthy individual E=i+1 W=i-1 N=j-1 S=j+1 # Check if outside the matrix if (E==n+1) { E=1} if (W==0) { W=n} if (N==0) { N=n} if (S==n+1) { S=1} # Count number of infected neighbors if(mat[E,j]==1){R=R+1} # East if(mat[W,j]==1){R=R+1} # West if(mat[i,N]==1){R=R+1} # North if(mat[i,S]==1){R=R+1} # South if(mat[E,N]==1){R=R+1} # North East if(mat[E,S]==1){R=R+1} # South East if(mat[W,N]==1){R=R+1} # North West if(mat[W,S]==1){R=R+1} # South West } a=-1.5 b=0.6 Pinfect=(1/(1+exp(-(a+b*R)))) # Calc probability for healthy to become infected g=runif(1) # Draw a random number between 0 and 1 if (g<Pinfect & mat[i,j]==0 & R>0){ temp[i,j]=1 # Healthy becomes infected } if (mat[i,j]==1){ # If infected individual g=runif(1) # Draw a random number between 0 and 1 if (g<q){ temp[i,j]=2 # Infected becomes Resistant } } } } image(mat,col=c("grey50","deeppink","seagreen1"),add=TRUE,,zlim=c(0,2)) # Plot image grid(nx=n,ny=n,col="grey70",lty=1) Sys.sleep(0.1) # To see movement on screen we need to pause the loop mat=temp # Overwrite matrix # Save number of healthy, infected and resistant at each time step n_helalthy[k]=kH n_infected[k]=kI n_resistant[k]=kR } graphics.off() plot(1:k,n_helalthy,type="l",ylab="Number",xlab="Time steps (weeks)",col=1,ylim=c(0,2600)) lines(1:k,n_infected,col=2) lines(1:k,n_resistant,col=3) legend(x=52,y=1599,c("Susceptible","Infected","Recovered"),lty=1,col=1:3)