library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(gganimate)
## Warning: package 'gganimate' was built under R version 4.3.3
As K sites are added and K-1 deleted at each step, \(|A_n|-|A_{n-1}|=1\), i.e., the number of sites increases by 1 per step. Hence \(|A_n|=n+1\) (\(A_0=1\)).
A=c(0);K=3;B=list();B[[1]]=A;
anim=data.frame("trial"=1,"proc"=B[[1]])
for(n in 1:100){
  x=length(A);A_ng=c();
  while(length(A)!=x+K){
    for(y in A){
   A_ng[length(A_ng)+1]=y+1
   A_ng[length(A_ng)+1]=y-1
    }
    A_bd=setdiff(A_ng,A)
    t=sample(A_bd,size=1)
    A=append(A,t)
    }
    d=sample(1:length(A),size=K-1,replace=FALSE)
    A=A[-d]
    B[[n+1]]=A
    anim=rbind(anim,data.frame("trial"=n+1,"proc"=B[[n+1]]))
}
x=1:length(A)
plot(x,A,type="b",col="#0D7072",cex=0.3)
ggplot(anim,fill="darkblue")+
  geom_point(aes(trial,proc))+transition_manual(trial,cumulative=TRUE)+ease_aes("sine")
L<- function(A) {
  A1=sort(A)
  A_g= diff(A1)-1
  sum(A_g)
}

G2<- function(A) {
  A1= sort(A)
  A_g= diff(A1)-1
  sum(A_g >=2)
}
lm=c();g2m=c();
for(m in 1:100){
LA=c();G2A=c();
A=c(0);K=3;
for(n in 1:100){
  x=length(A);A_ng=c();
  while(length(A)!=x+K){
    for(y in A){
   A_ng[length(A_ng)+1]=y+1
   A_ng[length(A_ng)+1]=y-1
    }
    A_bd=setdiff(A_ng,A)
    t=sample(A_bd,size=1)
    A=append(A,t)
    }
    d=sample(1:length(A),size=K-1,replace=FALSE)
    A=A[-d]
   LA[n]=L(A)
   G2A[n]=G2(A)
  }
lm[m]=mean(LA)
g2m[m]=mean(G2A)
}
x_axis=1:100
y_axis=lm[x_axis]/(log(x_axis)/log(log(x_axis)))
plot(x_axis,y_axis,type="l",xlab="n",ylab="L_n/h_n")

As is clear from the plot, \(|\frac{L_n}{h_n}|\ \leq\ 2\ \forall\ n.\) Hence, experimentally, \(L_n\ =\ O(h_n).\) Also, referring to the thesis on “Discrete Growth Models” by Dorothea M. Eberz-Wagner, \(\mathbb{P}[lim\ sup_{n\to\ \infty}\frac{L(A_n)}{h_n}=\frac{K-1}{K}]=1\), which further confirms the claim.

plot(1:100,g2m,type="l")