Source code of the figures in the paper "And yet it rocks! Fluctuations and growth in Ragnar Frisch's rocking horse model"

Vincent Carret - Frisch (1933): And yet it rocks.r

# Author: Vincent Carret
# Contact: vincent.carret@univ-lyon2.fr
# This R script reproduces the figures in the paper "And yet it rocks! Fluctuations and growth in Ragnar Frisch's rocking horse model" (version 3 - online 5 Dec. 2020)
# There is an app version accessible here: https://cbheem.shinyapps.io/Frisch/
# The libraries at the beginning are necessary and must be installed

library(plotly)
library(dplyr)

# Functions ----
# Functions to compute the Lambert W Function used in the Laplace transform solution ----
# This algorithm is based on Corless et al. (1996) and the c++ implementation by Istvan Mezo (https://github.com/IstvanMezo/LambertW-function)
initPoint <- function(z, k){
  I <- complex(real = 0, imaginary = 1)
  tPiKI <- complex(real = 0, imaginary = 2*pi*k)
  ip <- log(z) + tPiKI - log(log(z) + tPiKI)
  p <- sqrt(2*(exp(1)*z+1))
  
  if(abs(z - (-exp(-1))) <= 1){
    if(k == 0) ip <- -1 + p - 1/3 * p^2 + 11/72 * p^3
    if(k == 1 && Im(z) < 0) ip <- -1 - p - 1/3 * p^2 - 11/72 * p^3
    if(k == -1 && Im(z) > 0) ip <- -1 - p - 1/3 * p^2 - 11/72 * p^3
  }
  
  if(k == 0 && abs(z - 0.5) <= 0.5) ip <- (0.35173371 * (0.1237166 + 7.061302897 * z)) / (2 + 0.827184 * (1 + 2*z))
  
  if (k == -1 && abs(z - 0.5) <= 0.5) ip <- -(((2.2591588985 + 4.22096*I) * ((-14.073271 - 33.767687754*I) * z - (12.7127 - 19.071643*I) * (1 + 2*z))) / (2 - (17.23103 - 10.629721*I) * (1 + 2*z)))
  
  return(ip)
}

lambertW <- function(z, k = 0){
  if(class(z) != "complex") z <- complex(real = z, imaginary = 0)
  # Some values of z and k yield known results
  if(z == 0) return(if(k == 0) 0 else -Inf)
  if(z == -exp(-1) && (k == 0 || k == -1)) return(-1)
  if(z == exp(1) && k == 0) return(1)
  
  # Halley's method
  w <- initPoint(z, k)
  maxiter = 30
  prec <- 10^-30
  
  for(i in 1:maxiter){
    wprev <- w
    w <- w-(w*exp(w)-z)/(exp(w)*(w+1) - ((w+2)*(w*exp(w)-z))/(2*w+2))
    if(abs(w-wprev) < prec) break
  }
  return(w)
}

# Functions to compute the coefficients and roots of Frisch's model (in the case of a disturbed equilibrium) ----
coefs <- function(s,a,b,c,d,eps,disc,ic){
  if(s != 0) num <- c+s*disc+(d/s-b)*ic*(1-exp(-eps*s))-d*ic*eps
  else num <- c
  denom <- 2*s+a+b*exp(-eps*s)-b*s*eps*exp(-eps*s)+d*eps*exp(-eps*s)
  return(num/denom)
}

polyFrisch <- function(x,a,b,d,eps=6, deriv = 0){
  if(deriv == 0) return(x^2+a*x+b*x*exp(-eps*x)+d-d*exp(-eps*x))
  else if(deriv == 1) return(2*x+a+b*exp(-eps*x)-eps*b*x*exp(-eps*x)+eps*d*exp(-eps*x))
}

newton <- function(x,a,b,d,eps=6){
  return(x-polyFrisch(x,a,b,d,eps)/polyFrisch(x,a,b,d,eps, deriv = 1))
}

compute_roots <- function(rho){
  roots <- c(0)
  for(x in rho){
    for(i in 1:100){
      xp <- newton(x,a,b,d,eps)
      if(abs(x - xp) < 1e-10) break
      else x = xp
    }
    if(round(x,5) == 0){
      x = -1
      for(i in 1:100){
        xp <- newton(x,a,b,d,eps)
        if(abs(x - xp) < 1e-10) break
        else x = xp
      }
    }
    roots <- c(roots,x)
  }
  return(roots)
}

# Functions and parameters reproducing the step by step solution ----
iFT <- function(t,h) return(round(t/h+1))

yt = function(m,mu,x,t,h){
  return(m*x[iFT(t,h)]+mu*(x[iFT(t+h,h)]-x[iFT(t,h)])/h)
}

h = 1/6

# Figures ----
# Figure 1 & 2 ----
eps = 6
mu = 10
m = 0.5
lam = 0.05
r = 2
s = 1

N = 1000
ft = 75
t <- seq(from = 0, to = ft, by = 1/eps)

a = lam*r+lam*s*mu/eps
b = -lam*s*mu/eps
c = 0.165
d = lam*s*m/eps

rho = sapply(0:N, function(x) lambertW(-b*exp(a*eps)*eps, x)/eps-a)
roots = compute_roots(rho)

eq = c/(lam*(r+s*m))
disc = 1.1*eq
init_vals = eq
ks = c()
for(i in 1:N){
  ks <- c(ks, coefs(roots[i],a,b,c,d,eps,disc,init_vals))
}

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

if(N > 0){
  superpos <- rep(Re(ks[1]), length(t))
  if(N > 1){
    for(i in 2:N){
      if(i == 2){
        superpos <- superpos + Re(ks[i])*exp(t*Re(roots[i]))
      } else{
        superpos <- superpos + Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))
      }
    }
  }
  superpos = c(rep(init_vals,eps^2+1), superpos)
  fig = fig %>% add_trace(y = superpos, type = "scatter", mode = "lines", name = "components sum", color = I("black"))
}

fig

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

comps = c(2,3,4,5)
if(length(comps) != 0){
  for(i in as.numeric(comps)){
    if(i == 2){
      comp = c(rep(init_vals,eps^2), NA, init_vals+Re(ks[i])*exp(t*Re(roots[i])))
      fig <- fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = "Trend")
    } else{
      comp = c(rep(init_vals,eps^2), NA, init_vals+Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))) 
      fig = fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = paste0("Cycle ", i-2))
    }
  }
}

y1 = init_vals+Re(ks[2])*exp(t[1]*Re(roots[2]))
y2 = init_vals+Mod(2*ks[3])*exp(t[1]*Re(roots[3]))*cos(t[1]*Im(roots[3]) + Arg(2*ks[3]))

fig <- fig %>% add_segments(x = eps, xend = eps, y = y1, yend = y2, line = list(color = 'rgb(173, 173, 173)', width = 1, dash = 'dash'), showlegend=F) 
fig %>% layout(legend = list(x = 0.8,y=0.98, font = list(size=24)))

# Figure 3 & 4 ----
eps = 6
mu = 10
m = 0.5
lam = 0.3
r = 2
s = 1

N = 1000
ft = 50
t <- seq(from = 0, to = ft, by = 1/eps)

a = lam*r+lam*s*mu/eps
b = -lam*s*mu/eps
c = 0.165
d = lam*s*m/eps

rho = sapply(0:N, function(x) lambertW(-b*exp(a*eps)*eps, x)/eps-a)
roots = compute_roots(rho)

eq = c/(lam*(r+s*m))
disc = 1.1*eq
init_vals = eq
ks = c()
for(i in 1:N){
  ks <- c(ks, coefs(roots[i],a,b,c,d,eps,disc,init_vals))
}

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

if(N > 0){
  superpos <- rep(Re(ks[1]), length(t))
  if(N > 1){
    for(i in 2:N){
      if(i == 2){
        superpos <- superpos + Re(ks[i])*exp(t*Re(roots[i]))
      } else{
        superpos <- superpos + Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))
      }
    }
  }
  superpos = c(rep(init_vals,eps^2+1), superpos)
  fig = fig %>% add_trace(y = superpos, type = "scatter", mode = "lines", name = "components sum", color = I("black"))
}

fig

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

comps = c(2,3)
if(length(comps) != 0){
  for(i in as.numeric(comps)){
    if(i == 2){
      comp = c(rep(init_vals,eps^2), NA, init_vals+Re(ks[i])*exp(t*Re(roots[i])))
      fig <- fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = "Trend")
    } else{
      comp = c(rep(init_vals,eps^2), NA, init_vals+Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))) 
      fig = fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = paste0("Cycle ", i-2))
    }
  }
}

y1 = init_vals+Re(ks[2])*exp(t[1]*Re(roots[2]))
y2 = init_vals+Mod(2*ks[3])*exp(t*Re(roots[3]))*cos(t*Im(roots[3]) + Arg(2*ks[3]))

fig <- fig %>% add_segments(x = eps, xend = eps, y = y1, yend = y2, line = list(color = 'rgb(173, 173, 173)', width = 1, dash = 'dash'), showlegend=F) 
fig %>% layout(legend = list(x = 0.8,y=0.98, font = list(size=24)))

# Figure 5 & 6 ----
eps = 6
mu = 15
m = 1
lam = 0.3
r = 1
s = 2

N = 1000
ft = 100
t <- seq(from = 0, to = ft, by = 1/eps)

a = lam*r+lam*s*mu/eps
b = -lam*s*mu/eps
c = 0.165
d = lam*s*m/eps

rho = sapply(0:N, function(x) lambertW(-b*exp(a*eps)*eps, x)/eps-a)
roots = compute_roots(rho)

eq = c/(lam*(r+s*m))
disc = 1.1*eq
init_vals = eq
ks = c()
for(i in 1:N){
  ks <- c(ks, coefs(roots[i],a,b,c,d,eps,disc,init_vals))
}

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

if(N > 0){
  superpos <- rep(Re(ks[1]), length(t))
  if(N > 1){
    for(i in 2:N){
      if(i == 2){
        superpos <- superpos + Re(ks[i])*exp(t*Re(roots[i]))
      } else{
        superpos <- superpos + Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))
      }
    }
  }
  superpos = c(rep(init_vals,eps^2+1), superpos)
  fig = fig %>% add_trace(y = superpos, type = "scatter", mode = "lines", name = "components sum", color = I("black"))
}

fig

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

comps = c(2,3)
if(length(comps) != 0){
  for(i in as.numeric(comps)){
    if(i == 2){
      comp = c(rep(init_vals,eps^2), NA, init_vals+Re(ks[i])*exp(t*Re(roots[i])))
      fig <- fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = "Trend")
    } else{
      comp = c(rep(init_vals,eps^2), NA, init_vals+Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))) 
      fig = fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = paste0("Cycle ", i-2))
    }
  }
}

y1 = init_vals+Re(ks[2])*exp(t[1]*Re(roots[2]))
y2 = init_vals+Mod(2*ks[3])*exp(t*Re(roots[3]))*cos(t*Im(roots[3]) + Arg(2*ks[3]))

fig <- fig %>% add_segments(x = eps, xend = eps, y = y1, yend = y2, line = list(color = 'rgb(173, 173, 173)', width = 1, dash = 'dash'),showlegend=F) 
fig %>% layout(legend = list(x = 0.8,y=0.98, font = list(size=24)))

# Figure 7 & 8 ----
eps = 6
mu = 15
m = 1
lam = 0.3
r = 1
s = 2

N = 1000
ft = 100
t <- seq(from = 0, to = ft, by = 1/eps)

a = lam*r+lam*s*mu/eps
b = -lam*s*mu/eps
c = 0.165
d = lam*s*m/eps

rho = sapply(0:N, function(x) lambertW(-b*exp(a*eps)*eps, x)/eps-a)
roots = compute_roots(rho)

eq = c/(lam*(r+s*m))
disc = 0.2*eq
init_vals = 0.2*eq
ks = c()
for(i in 1:N){
  ks <- c(ks, coefs(roots[i],a,b,c,d,eps,disc,init_vals))
}

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

if(N > 0){
  superpos <- rep(Re(ks[1]), length(t))
  if(N > 1){
    for(i in 2:N){
      if(i == 2){
        superpos <- superpos + Re(ks[i])*exp(t*Re(roots[i]))
      } else{
        superpos <- superpos + Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))
      }
    }
  }
  superpos = c(rep(init_vals,eps^2+1), superpos)
  fig = fig %>% add_trace(y = superpos, type = "scatter", mode = "lines", name = "components sum", color = I("black"))
}

fig

to_plot = data.frame(time = c(t,seq(ft,ft + eps,by=1/eps)))
fig = plot_ly(to_plot, x = ~time)

comps = c(2,3,4)
if(length(comps) != 0){
  for(i in as.numeric(comps)){
    if(i == 2){
      comp = c(rep(init_vals,eps^2), NA, Re(ks[1])+Re(ks[i])*exp(t*Re(roots[i])))
      fig <- fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = "Trend")
    } else{
      comp = c(rep(init_vals,eps^2), NA, init_vals+Mod(2*ks[i])*exp(t*Re(roots[i]))*cos(t*Im(roots[i]) + Arg(2*ks[i]))) 
      fig = fig %>% add_trace(y = comp, type = 'scatter', mode = 'lines', name = paste0("Cycle ", i-2))
    }
  }
}

y1 = Re(ks[1])+Re(ks[2])*exp(t[1]*Re(roots[2]))
y2 = init_vals+Mod(2*ks[3])*exp(t*Re(roots[3]))*cos(t*Im(roots[3]) + Arg(2*ks[3]))

fig <- fig %>% add_segments(x = eps, xend = eps, y = y1, yend = y2, line = list(color = 'rgb(173, 173, 173)', width = 1, dash = 'dash'),showlegend=F) 
fig %>% layout(legend = list(x = 0.8,y=0.9, font = list(size=24)))

# Using the step-by-step solution
tfinal = 100
times = seq(0,tfinal,by=h)
x = rep(init_vals,length(times))

start = iFT(eps,h)
x[start] = disc

for(t in times[(start):iFT(tfinal,h)]){
  denom = (1+lam*s*h*mu/(2*eps))
  sumy = sum(sapply(seq(h,eps-h,by=h), function(i) yt(m,mu,x,t-i,h)))
  x[iFT(t+h,h)] = (h*c+(1-lam*h*r-(h*lam*s*(h*m-mu))/(2*eps))*x[iFT(t,h)]-lam*s*h^2/eps*(yt(m,mu,x,t-eps,h)/2+sumy))/denom
}

x = x[1:length(x)-1]
plot(times, x, type = "l")
lines(times,superpos[1:601], col = "red")