Skip to content

Commit

Permalink
more action tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@thomasmanke
thomasmanke committed Mar 13, 2024
1 parent 64b18f7 commit e363d78
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 165 deletions.
Binary file added posts/Lorenz/lorenz.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
198 changes: 33 additions & 165 deletions posts/Planets_test/Planets_test.qmd
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
---
title: Planetary Motion (Test)
title: Lorentz (Test)
date: 2024/03/02
date-modified: last-modified
categories:
Expand All @@ -23,9 +23,6 @@ filters:

## Lorenz Attractor (1963)

This blog shows how the `deSolve` package can be used to solve non-linear ODE numerically.
Here at the example of the Lorenz Attractor:

$$
\begin{eqnarray*}
\dot{X} = a (Y - X) \\
Expand All @@ -34,63 +31,14 @@ $$
\end{eqnarray*}
$$

This set of equation is meant to model convection in a 2D fluid layer.

- $X \propto$ rate of convection, $a = \sigma \propto$ Prandtl number
- $Y \propto$ horizontal temperature variation ($b = \rho \propto$ Rayleigh number)
- $Z \propto$ vertical temperature variation ($c = \beta \propto$ layer dimensions)

For specific choices of parameters, this system is known to be chaotic: $a=10, b=28, c=8/3$



## Lorenz and climate

The app can calculate many initial conditions simultaenously.
This illustrates the *butterfly effect*: For chaotic systems, small differences in initial conditions
can have large differences after some time $t$.


## The app
```{shinylive-r}
#| standalone: true
#| viewerHeight: 600
library(shiny)
library(dplyr)
library(tidyr)
library(deSolve)
library(plotly)
library(DT)
library(shinyalert)
### Define Planetary data #####
planets = data.frame(
name=c("Mercury", "Venus", "Earth", "Mars"),
a=c(0.386, 0.724, 1.00, 1.52), # semi-major axis in AU (for Earth: a=1.00)
e=c(0.2056, 0.068, 0.0167, 0.0934), # eccentricity
w=c(0.241, 0.615, 1.00, 1.00), # orbital velocities - will be overwritten by Kepler
show=c(T, T, T, T) # include/exclude in plot
)
# semi-minor axis b from eccentricity: b = a*(1-e^2)
# orbital velocity w from Kepler's law: T ~ sqrt(a^3)
planets = planets %>% mutate(b=a*(1 - e^2), w = 1/sqrt(a^3))
rownames(planets)=planets$name # assign rownames
### Define Ellipses #####
createEllipses <- function(Tmax=360, dT=20,params=planets){
t = seq(from=0,to=Tmax, by=dT) # timesteps in days
d = data.frame(x=0,y=0,day=t, name="Sun") # sun in center (0,0) for all t
for (name in rownames(params)) {
if (!params[name,"show"]) { next } # skip if planet unselected
a = params[name, "a"]
b = params[name, "b"]
w = params[name, "w"]
phi <- 2*pi*w*t/365.242 # angles for all time t
d <- rbind(d, data.frame(x = a*cos(phi), y = b*sin(phi), day=t, name=name))
}
return(d)
}
library(tidyverse)
# The lorenz functions in a for deSolve can understand
Lorenz<-function(t, state, parameters) {
Expand All @@ -103,130 +51,50 @@ Lorenz<-function(t, state, parameters) {
})
}
#solveIni <- function(state, ts=1:20, ps=parameters) {
# # solve ODE and convert to data frame
# df = data.frame()
# df
# ode(y = state, times = ts, func=Lorenz, parms=ps) %>% as.data.frame()
#}
# state2states <- function(state, N = 5, sd = 0.1){
# # input: initial state = named vector of length p (number of variables)
# # output: matrix of N row vectors of perturbed initial states
# p <- length(state)
# M <- matrix(rnorm(N * p, mean = 1, sd = sd), nrow = N)
# states <- M*state
# colnames(states) = names(state)
# return(states)
# }
ui <- fluidPage(
titlePanel('Planetary Motion'),
ui <- fluidPage(
titlePanel("Lorenz Attractor"),
fluidRow(
column(2,
numericInput('Tmax', 'Tmax', min=10, max=1000,step=1, value=365),
),
column(2,
numericInput('dT', 'dT', min=1, max=30,step=1, value=10.0),
),
column(2,
numericInput('ms', 'ms per frame', min=0, step=20, value=100),
),
column(6,
selectInput('center', 'choose center', c("Sun",rownames(planets))),
),
column(3, sliderInput("param_a", label = "a:", min = -10, max = 20, value = 10, step=0.001)),
column(3, sliderInput("param_b", label = "b:", min = -10, max = 100, value = 28, step=0.001)),
column(3, sliderInput("param_c", label = "c:", min = -10, max = 20, value = 8./3., step=0.001)),
column(3, sliderInput("time", label = "Time Range:", min = 1, max = 200, value = c(1,100), step=1)),
column(12, actionButton("update", "Update Parameters") ),
column(3, sliderInput("ini_n", label = "ini_n:", min = 1, max = 1000, value = 50, step=1)),
column(3, sliderInput("ini_sd", label = "ini_sd:", min = 0, max = 10, value = 1, step=0.01)),
column(6, sliderInput("ms", label = "ms / frame", min = 0, max = 500, value = 200, step=1)),
tabsetPanel(
tabPanel("Orbits", plotlyOutput('planets',width = "600px", height = "300px")),
tabPanel("Parameters", DTOutput("params", width="50%")),
tabPanel("Session Info", verbatimTextOutput("session_info"))
# tabPanel("Lorenz Attractor",plotlyOutput("animated",width="600px",height="600px")),
tabPanel("(X,Y,Z) vs Time",plotOutput("time_series")),
tabPanel("Session Info",verbatimTextOutput("session_info"))
)
)
)
######
server<-function(input,output,session){
options(warn = -1)
rv <- reactiveValues(
ellipses = createEllipses(360, 20),
planets = planets
)
# observe update button
observeEvent(input$update, {
n_frames <- input$Tmax / input$dT
msg1 <- paste("Number of frames = ", n_frames)
msg2 <- 'This may take some time'
if (n_frames > 80) {
shinyalert(msg1, msg2, type = "info")
}
rv$ellipses <- createEllipses(input$Tmax, input$dT, rv$planets)
})
# observe editing of parameters
observeEvent(input$params_cell_edit, {
info = input$params_cell_edit
str(info)
planets <<- editData(planets, input$params_cell_edit, 'params')
rv$planets <- planets
# ensure that logical stays logical
rv$planets <- rv$planets %>% mutate(show=as.logical(show))
})
### recenter ellipses if center is changed by user
recenter_ellipses <- reactive({
center <- input$center
df_r <- rv$ellipses
# n_obj = number of objects in df_r (planets + sun)
# n_obj may change in response to parameter editing --> input$params_cell_edit
n_obj <- df_r %>% select(name) %>% n_distinct()
df_c <- df_r %>% filter(name==center) # get all coordinates for new center
# expand df_c n_obj times (to same length as df_r)
# here is a strong assumption that the times will match
# could be generalized by proper pivoting
df_c <- df_c[rep(seq_len(nrow(df_c)), n_obj), ]
df_r[,1:2] <- df_r[,1:2] - df_c[,1:2] # recenter
df_r
})
### render plot: careful not to use too many frames (>100) --> very slow
output$planets<-renderPlotly({
# heliocentric view
helio <- rv$ellipses %>%
plot_ly(x = ~x, y = ~y, color = ~name, type='scatter', mode='lines') %>%
add_trace(mode='markers', size=18, frame = ~day, showlegend=FALSE)
# other view
other <- recenter_ellipses() %>%
plot_ly(x = ~x, y = ~y, color = ~name, type='scatter', mode='lines', showlegend=FALSE) %>%
add_trace(mode='markers', size=18, frame = ~day, showlegend=FALSE)
subplot(helio, other) %>% animation_opts(frame=input$ms, transition=input$ms, redraw=FALSE)
})
### parameter table - editable
output$params <- renderDT(rv$planets, selection='none', editable='cell')
### sessionInfo()
output$session_info <- renderPrint( sessionInfo() )
}
vals <- reactiveValues()
observe({
parameters <- c(a = input$param_a, b = input$param_b, c = input$param_c)
times <- seq(input$time[[1]], input$time[[2]], by = 0.01)
state <- c(X = 1, Y = 1, Z = 1) # single initial state
# single solution - high time resolution
vals$out <- ode(y = state, times = times, func = Lorenz, parms = parameters) %>%
as.data.frame()
})
output$time_series <- renderPlot({
vals$out %>%
pivot_longer(-time, names_to= "coord" , values_to = "value") %>%
ggplot(aes(x = time, y = value)) +
geom_path() +
facet_wrap(~coord)
})
output$session_info <- renderPrint( sessionInfo() )
}
shinyApp(ui = ui, server = server)
```
Expand Down

0 comments on commit e363d78

Please sign in to comment.