Dynamical System Case Study 1 (symmetric 3d system)



Let’s consider the following dynamical system:

\[\begin{cases} \dot{x}_1 = f(x_3) - x_1,\\ \dot{x}_2 = f(x_1) - x_2,\\ \dot{x}_3 = f(x_2) - x_3, \end{cases} \]

where \(f(x) = \alpha / (1+x^m)\) is a Hill function. In this case study, we explore the phase portrait of this system for \(\alpha = 18,\; m = 3\).

Preparation

First of all, let’s find the stationary point of this system. Since in such point \(\dot{x}_1=\dot{x}_2=\dot{x}_3=0\), we have the following equation:

\[x_1 = f(f(f((x_1)))). \]

Since \(x_1\) is monotonically increasing and \(f(f(f(x_1))))\) is monotonically decreasing, this equation has exactly one solution. It’s easy to see that this solution is \(x_0=x_1=x_2=x_3=2\).

Next, we build the corresponding linearization matrix:

\[M = \begin{bmatrix} -1 & 0 & f'(x_0)\\ f'(x_0) & -1 & 0\\ 0 & f'(x_0) & -1 \end{bmatrix}. \]

Since \(f'(x) = -\alpha m x^{m - 1} / (x ^ m + 1)^2\), we have:

\[f'(x_0) = -18 * 3 * 2^2 / (2^3+1)^2 = -216 / 81 = -8/3 \approx -2.66667. \]

Thus,

\[M = \begin{bmatrix} -1 & 0 & -8/3\\ -8/3 & -1 & 0\\ 0 & 8/3 & -1 \end{bmatrix}. \]

Now we can get the eigenvalues \(\lambda_i\) and eigenvectors \(e_i\):

\[\lambda_1 \approx -3.666667,\quad \lambda_2 \approx 0.333333+2.309401i,\quad \lambda_3 \approx 0.333333-2.309401i, \]

\[\begin{bmatrix} e_1 \\ e_2 \\ e_3 \end{bmatrix} = \begin{bmatrix} -0.5773503, & -0.5773503, & -0.5773503\\ -0.2886751+0.5i, & -0.2886751-0.5i, & 0.5773503\\ -0.2886751-0.5i, & -0.2886751+0.5i, & 0.5773503 \end{bmatrix}. \]

The most expressive phase portrait could be obtained using a projection on \((e_2, e_3)\) or \((\Re(e_2), \Im(e_2)) = (\Re(e_3), -\Im(e_3))\).

Phase portraits

In order to explore phase portraits, we perform several simulation studies. In each simulation, we generate several trajectories from random start points and draw three plots:

  1. A projection of simulated trajectories on \((\Re(e_2), \Im(e_2))\). Each trajectory has own color from a rainbow palette. The limit cycle is shown using the black color.
  2. A projection of simulated trajectories on \((\Re(e_2), \Im(e_2))\). For each trajectory point, the color defines the current speed (\(\sqrt{\dot{x}_1^2+\dot{x}_2^2+\dot{x}_3^2}\)). The limit cycle is shown using the black color. The beginning of each trajectory is drop in order to provide a more consistent picture.
  3. A 3D projection of simulation trajectories on \((e1_, \Re(e_2), \Im(e_2))\). These visualizations are heavy, so you should click on “Show 3D plot” in order to explore the plot.

Simulation 1 (10 trajectories)

In this simulation, we generate 10 trajectories from random start points around the stationary point uniformly taken in \((x_1, x_2, x_3)\).




Simulation 2 (100 trajectories)

In this simulation, we generate 100 trajectories from random start points around the stationary point uniformly taken in \((x_1, x_2, x_3)\).




Simulation 3 (500 trajectories)

In this simulation, we generate 500 trajectories from random start points around the stationary point uniformly taken in \((e_1, e_2, e_3)\).




Source code

library(deSolve)
library(ggplot2)
library(dplyr)
library(tidyr)
library(plotly)

simParams <- list(
  totalTime = 50,
  step = 0.01,
  startThreshold = 3,
  endThreshold = 45
)

hypot <- function(x, y) sqrt(x^2 + y^2)
fp <- list(a = 18, m = 3)
f <- function(x, a = fp$a, m = fp$m) a / (1 + x ^ m)
df <- function(x, a = fp$a, m = fp$m) -a * m * x ^ (m - 1) / (x ^ m + 1)^2

model <- function(t, x, parms) {
  with(as.list(c(parms, x)), {
    dx1 <- f(x3) - x1
    dx2 <- f(x1) - x2
    dx3 <- f(x2) - x3
    list(c(dx1, dx2, dx3))
  })
}
calc.speed <- function(x) sqrt(sum(unlist(model(0, x, 0))^2))

x0 <- uniroot(function(x) f(f(f(x)))-x, c(0, fp$a), tol = 1e-9)$root
M <- matrix(c(
  -1, 0, df(x0),
  df(x0), -1, 0,
  0, df(x0), -1
), nrow = 3, ncol = 3, byrow = TRUE)
eigenM <- tryCatch(eigen(M), error = function(e) list(values = rep(NA, 6), vectors = matrix(rep(NA, 36), ncol = 6)))
Ml <- eigenM$values
Me <- eigenM$vectors

simulate <- function(name, start) {
  times  <- seq(0, simParams$totalTime, by = simParams$step)
  traj <- data.frame(lsoda(start, times, model, c()))
  traj$speed <- sapply(1:nrow(traj), function(index) calc.speed(traj[index, 2:4]))

  proj <- function(e){
    dotProduct <- rowSums(t(t((traj[,2:4]) - x0) * e))
    l <- sqrt(sum(e^2))
    dotProduct / l
  }
  traj$e1 <- proj(Re(Me[,1]))
  traj$e2 <- proj(Re(Me[,2]))
  traj$e3 <- proj(Im(Me[,2]))
  traj$name <- name
  traj
}
generate.input <- function(type, maxTraj = -1) {
  if (type == 1) {
    if (maxTraj <= 0)
      maxTraj <- 1
    return(data.frame(
      x1 = abs(x0 + rnorm(maxTraj) * 3),
      x2 = abs(x0 + rnorm(maxTraj) * 3),
      x3 = abs(x0 + rnorm(maxTraj) * 3)
    ))
  }
  if (type == 2) {
    input <- expand.grid(e2 = seq(-5, 8, by = 0.2), e3 = seq(-7.5, 7, by = 0.2))
    noise <- runif(nrow(input), -5, 5)
    ue1 <- Re(Me[,1])
    ue2 <- Re(Me[,2])
    ue3 <- Im(Me[,2])
    input$x1 <- noise * ue1[1] + input$e2 * ue2[1] + input$e3 * ue3[1]
    input$x2 <- noise * ue1[2] + input$e2 * ue2[2] + input$e3 * ue3[2]
    input$x3 <- noise * ue1[3] + input$e2 * ue2[3] + input$e3 * ue3[3]
    input <- input[input$x1 > 0 & input$x2 > 0 & input$x3 > 0,]
    if (maxTraj > 0 & maxTraj < nrow(input))
      input <- input[sample(1:nrow(input), maxTraj),]
    return(input)
  }
  stop("Unknown type")
}
simulate2 <- function(input) {
  get.input <- function(i) c(x1 = input[i, "x1"], x2 = input[i, "x2"], x3 = input[i, "x3"])
  do.call("rbind", lapply(1:nrow(input), function(i) simulate(paste0("t", i), get.input(i))))
}

draw.simple <- function(traj) {
  ggplot(traj, aes(e2, e3, col = name)) +
    geom_path(aes(group = name)) +
    geom_point(aes(x, y), data.frame(x = 0, y = 0), col = "red") +
    geom_path(aes(e2, e3, group = name), traj[traj$time > simParams$endThreshold,], col = "black", size = 2, alpha = 0.1) +
    scale_color_manual(values = sample(rainbow(length(unique(traj$name))))) +
    theme(legend.position="none") +
    labs(x = "Re(e2)", y = "Im(e2)")
}
draw.speed <- function(traj) {
  ggplot(traj[traj$time > simParams$startThreshold,], aes(e2, e3, col = speed)) +
    geom_path(aes(group = name)) +
    geom_point(aes(x, y), data.frame(x = 0, y = 0), col = "red") +
    geom_path(aes(e2, e3, group = name), traj[traj$time > simParams$endThreshold,], col = "black", size = 2, alpha = 0.1) +
    scale_color_gradientn(colours = rainbow(7)) +
    labs(x = "Re(e2)", y = "Im(e2)")
}
draw.3d <- function(traj) {
  traj3d <- traj[traj$time > simParams$startThreshold,c("time", "e1", "e2", "e3", "name")]
  traj3d[,2:4] <- round(traj3d[,2:4], 3) # Compressing
  plot_ly(
    traj3d, x = ~e1, y = ~e2, z = ~e3,
    type = 'scatter3d', mode = 'lines', opacity = 1, color = ~name,
    colors = sample(rainbow(length(unique(traj$name))))) %>%
    layout(showlegend = FALSE)
}


# Simulation 1
set.seed(1)
traj1 <- simulate2(generate.input(1, 10))
# draw.simple(traj1)
# draw.speed(traj1)
# draw.3d(traj1)

# Simulation 2
set.seed(2)
traj2 <- simulate2(generate.input(1, 100))
# draw.simple(traj2)
# draw.speed(traj2)
# draw.3d(traj2)

# Simulation 3
set.seed(3)
traj3 <- simulate2(generate.input(2))
# draw.simple(traj3)
# draw.speed(traj3)
# draw.3d(traj3)