Simulation (2018-12-15)

model.R

################################################################################
# A simulation of the following system:
#
# dx/dt = L1(z) - k1 * x
# dy/dt = L2(x) - k2 * y
# dz/dt = L3(y) - k3 * z
#
# Parameters:
#   k1 = 1.2  k2 = 1.4  k3 = 1.5
#   a1 = 2    a2 = 3    a3 = 4
# Functions:
#  L1(p) = | k1 * a1 for p <= 1
#          | 0       for p >  1
#  L2(p) = | k2 * a2 for p <= 1
#          | 0       for p >  1
#  L3(p) = | k3 * a3 for p <= 1
#          | 0       for p >  1
#
# Author: Andrey Akinshin
################################################################################

library(deSolve)
library(ggplot2)
library(tidyr)
library(rgl)

L1 <- function(p, k1, a1) ifelse(p <= 1, k1 * a1, 0)
L2 <- function(p, k2, a2) ifelse(p <= 1, k2 * a2, 0)
L3 <- function(p, k3, a3) ifelse(p <= 1, k3 * a3, 0)

model <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {

    dx <- L1(z, k1, a1) - k1 * x
    dy <- L2(x, k2, a2) - k2 * y
    dz <- L3(y, k3, a3) - k3 * z
    list(c(dx, dy, dz))
  })
}

parameters <- c(k1 = 1.2, k2 = 1.4, k3 = 1.5,
                a1 = 2  , a2 = 3,   a3 = 4)
state      <- c(x = 0, y = 0, z = 0)
times      <- seq(0, 10, by = 0.01)

out <- ode(y = state, times = times, func = model, parms = parameters)

traj <- data.frame(
  time = out[,"time"],
  x    = out[,"x"],
  y    = out[,"y"],
  z    = out[,"z"]
)

simPlot <-
  traj %>%
  gather("var", "value", 2:4) %>%
  ggplot(aes(x=time, y=value, group=var, colour=var)) +
  geom_line() +
  theme_bw() +
  xlab("Time") +
  ylab("")
simPlot
ggsave("plot-2d.png")

plot3d(traj[,c(2,3,4)], axes=T, type="l", col="red")