Dynamical System Case Study 2 (Piecewise linear LLL-system)

Andrey Akinshin · 2022-07-17

We consider the following dynamical system:

$$ \begin{cases} \dot{x}_1 = L(a_1, k_1, x_3) - k_1 x_1,\\ \dot{x}_2 = L(a_2, k_2, x_1) - k_2 x_2,\\ \dot{x}_3 = L(a_3, k_3, x_2) - k_3 x_3, \end{cases} $$

where $L$ is a piecewise linear function:

$$ L(a, k, x) = \begin{cases} ak & \quad \textrm{for}\; 0 \leq x \leq 1,\\ 0 & \quad \textrm{for}\; 1 < x. \end{cases} $$

In this case study, we build a Shiny application that draws 3D phase portraits of this system for various sets of input parameters.

Shiny app parameters

The application allows controlling the following parameters:

$a_1$, $a_2$, $a_3$, $k_1$, $k_2$, $k_3$: the dynamical system parameter.
$\textrm{TotalTime}$: the total time simulation.
$\textrm{SkipTime}$: the initial time interval that will not be presented on the plot.
$N$: the number of simulated trajectories.
$\textrm{Seed}$: the randomization seed that controls the initial positions of all the trajectories.

We uniformly generate the initial point for each simulated trajectory from the following area:

$$ [0.5; (2+a_1)/3] \times [0.5; (2+a_2)/3] \times [0.5; (2+a_3)/3]. $$

Source code

The source code is also available on GitHub: https://github.com/AndreyAkinshin/pwLLL

core.R:

library(deSolve)

simulate <- function(context) {
  with(context, {
    L <- function(a, k, x) ifelse(x < 1, a * k, 0)

    model <- function(t, state, parms) {
      x1 <- as.numeric(state[1])
      x2 <- as.numeric(state[2])
      x3 <- as.numeric(state[3])
      list(c(
        dx1 = L(a1, k1, x3) - k1 * x1,
        dx2 = L(a2, k2, x1) - k2 * x2,
        dx3 = L(a3, k3, x2) - k3 * x3
      ))
    }

    times.step <- 0.01
    times.skip <- min(TotalTime, SkipTime)
    times <- seq(0, TotalTime, by = 0.01)

    simulate <- function(index) {
      get.start <- function(a) runif(1, 0.5, (2 + a) / 3)
      start <- c(x1 = get.start(a1), x2 = get.start(a2), x3 = get.start(a2))
      traj <- as.data.frame(lsoda(start, times, func = model, parms = 0))
      traj$name <- paste0("t", index)
      traj <- traj[round(times.skip / times.step):nrow(traj),]
      traj
    }
    set.seed(Seed)
    traj <- do.call("rbind", lapply(1:N, function(i) simulate(i)))

    list(
      context = context,
      model = model,
      traj = traj
    )
  })
}

randomCodeName <- function() {
  consonants <- c("b", "c", "d", "f", "g", "h", "j", "k", "l", "m", "n",
                  "p", "q", "r", "s", "t", "v", "w", "x", "y", "z")
  vowels <- c("a", "e", "i", "o", "u")
  paste0(sample(consonants, 1), sample(vowels, 1),
         sample(consonants, 1), sample(vowels, 1),
         sample(consonants, 1))
}

server.R:

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

source("core.R")

shinyServer(function(input, output, session) {

  output$plot3D <- renderPlotly({
    req(input$a1)
    req(input$a2)
    req(input$a3)
    req(input$k1)
    req(input$k2)
    req(input$k3)
    req(input$TotalTime)
    req(input$SkipTime)
    req(input$N)
    req(input$Seed)

    context <- list(
      a1 = input$a1,
      a2 = input$a2,
      a3 = input$a3,
      k1 = input$k1,
      k2 = input$k2,
      k3 = input$k3,
      TotalTime = input$TotalTime,
      SkipTime = input$SkipTime,
      N = input$N,
      Seed= input$Seed
    )
    res <- simulate(context)
    traj <- res$traj
    plot_ly(
      traj, x = ~x1, y = ~x2, z = ~x3,
      type = 'scatter3d', mode = 'lines', opacity = 1, color = ~name,
      colors = sample(rainbow(length(unique(traj$name))))) %>%
      layout(showlegend = FALSE, width = 900, height = 600)
  })

  output$htmlOutput <- renderUI({
    context <- reactiveValuesToList(input)
    rnd <- function(x) { if (is.complex(x) & abs(Im(x)) < 1e-9) round(Re(x), 6) else round(x, 6) }
    fmt <- function(x) { if (is.na(x)) "         ?" else formatC(x, digits = 6, width = 10, format = "f") }
    hr <- "-----------------------------------------------------------------<br>"
    text <- with(context, {
      paste0(
        "<pre>",
        "System:<br>",
        "dx1 = L(a1, k1, x3) - k1 * x1<br>",
        "dx2 = L(a2, k2, x1) - k2 * x2<br>",
        "dx3 = L(a3, k3, x2) - k3 * x3<br>",
        "L(a, k, x) = { a * k (for 0≤x≤1)<br>",
        "             { 0     (for 1&lt;x)<br>",
        "<br>",
        "Parameters:<br>",
        "a1 = ", a1, "<br>",
        "a2 = ", a2, "<br>",
        "a3 = ", a3, "<br>",
        "k1 = ", k1, "<br>",
        "k2 = ", k2, "<br>",
        "k3 = ", k3, "<br>",
        "TotalTime = ", TotalTime, "<br>",
        "SkipTime = ", SkipTime, "<br>",
        "N = ", N, "<br>",
        "Seed = ", Seed, "<br>",
        "</pre>"
        )
    })
    HTML(text)
  })

  output$save_state <- downloadHandler(
    filename = function() {
      paste0("data-", Sys.Date(), "-", randomCodeName(), ".yaml")
    },
    content = function(file) {
      data <- list(
        a1 = input$a1,
        a2 = input$a2,
        a3 = input$a3,
        k1 = input$k1,
        k2 = input$k2,
        k3 = input$k3,
        TotalTime = input$TotalTime,
        SkipTime = input$SkipTime,
        N = input$N,
        Seed = input$Seed
      )
      data <- data[names(data) != "restore_state"]
      list.save(data, file)
    }
  )

  loadedData <- reactive({
    list.load(input$restore_state$datapath)
  })
  observe({
    data <- loadedData()
    updateSliderInput(session, "a1", value = data$a1)
    updateSliderInput(session, "a2", value = data$a2)
    updateSliderInput(session, "a3", value = data$a3)
    updateSliderInput(session, "k1", value = data$k1)
    updateSliderInput(session, "k2", value = data$k2)
    updateSliderInput(session, "k3", value = data$k3)
    updateSliderInput(session, "TotalTime", value = data$TotalTime)
    updateSliderInput(session, "SkipTime", value = data$SkipTime)
    updateSliderInput(session, "N", value = data$N)
    updateSliderInput(session, "Seed", value = data$Seed)
  })
})

ui.R:

library(shiny)
library(plotly)

slider <- function(title, min, max, value = min) {
  sliderInput(title, title, min, max, value, 0.1)
}
sliderA <- function(title, value = 10) slider(title, 0, 20, value)
sliderK <- function(title, value = 1) slider(title, 0, 20, value)

shinyUI(fluidPage(

  titlePanel("Piecewise linear 3D system (LLL)"),

  sidebarLayout(
    sidebarPanel(
      div(style = "overflow-y:scroll; max-height: 700px; position:relative;",
      h2("Parameters"),

      sliderA("a1", 2.5),
      sliderA("a2", 3),
      sliderA("a3", 5),

      sliderK("k1", 1),
      sliderK("k2", 2.5),
      sliderK("k3", 4),

      hr(),

      sliderInput("TotalTime", "TotalTime", 1, 100, 10, 1),
      sliderInput("SkipTime", "SkipTime", 0, 20, 0.2, 0.1),
      sliderInput("N", "N (number of trajectories)", 1, 100, 10, 1),
      sliderInput("Seed", "Randomization seed", 1, 100, 1, 1),

      br(),
      br(),
      downloadButton("save_state", "Save parameters to file"),
      br(),
      br(),
      fileInput("restore_state", "Load parameters from file",
                placeholder = ".yaml file")
    )),

    mainPanel(
      plotlyOutput("plot3D", width = "900px", height = "600px"),
      htmlOutput("htmlOutput")
    )
  )
))