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

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)
}
)

list.load(input$restore_state$datapath)
})
observe({
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(),
br(),
br(),
placeholder = ".yaml file")
)),

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


The source code of this post and all the relevant files are available on GitHub.
Share: