Skip to content

[R-sig-dyn-mod] How to change a parameter after steady?

5 messages · Jinsong Zhao, Thomas Petzoldt

Original
Jinsong Zhao
# 
Hi there,

I have established a model. I set the times to a enough long period, so 
that the model runs into steady state. After the model reaching steady, 
I hope to change a or several parameters and run the model for a fixed 
periods. I expect to see some changes of state variables or auxiliary 
variables.

How to implement it? Thanks a lot!

Best wishes,
Jinsong
Thomas Petzoldt
# 
Hi,

there are several methods to implement this, so it would help us if you
could provide a *minimal* reproducible example.

Thomas
On 4/25/2014 9:37 PM, Jinsong Zhao wrote:
Jinsong Zhao
#
On 2014/4/26 3:25, Thomas Petzoldt wrote:
Can the following example be used here? Assumed the model reached steady 
when times >= 100. Then I hope to change r to 0.12, or K to 15. How to? 
May forcings in the previous post work? Any other solution?

Thanks again.

Best,
Jinsong

library(deSolve)
model <- function (time, y, parms) {
     with(as.list(c(y, parms)), {
          dN <- r * N * (1 - N / K)
          list(dN)
     })
}
y <- c(N = 0.1)
parms <- c(r = 0.1, K = 10)
times <- seq(0, 100, 1)
out <- ode(y, times, model, parms)

plot(out)
Thomas Petzoldt
#
On 4/27/2014 8:13 AM, Jinsong Zhao wrote:
Yes, forcings can be used if your aim is to create a "long-term" 
simulation with changing parameters. Just make r and K time dependent 
forcing functions instead of fixed parameters.

Nevertheless, I still don't understand why you want to do this. If your
aim is to start a model from an equilibrium, you may consider to
identify the equilibrium point with a solver from package rootSolve and
then to use this as the initial value of your simulation.

It is also possible to run several simulations with different parameters
independently, but that's another story.

Thomas


model <- function (time, y, parms) {
     with(as.list(c(y, parms)), {
          dN <- r * N * (1 - N / K)
          list(dN)
     })
}


## find equilibrium
library(rootSolve)
y <- c(N = 0.1)
parms <- c(r = 0.1, K = 10)
times <- c(0, Inf)
equilibrium <- steady(y, times, model, parms, method="runsteady")


## change parameters and run simulation
library(deSolve)
y <- equilibrium$y
parms <- c(r = 0.12, K = 15)
times <- seq(0, 100, 1)
out <- ode(y, times, model, parms)

plot(out)
Jinsong Zhao
#
On 2014/4/27 2:08, Thomas Petzoldt wrote:
Thank you very much for the detailed explanation and code demo.

I had tried the solution provided by rootSolve. However, I was misled by 
the plot that magnified the tiny difference in numeric of the output. I 
should take care about the numeric result.

Thanks again!

Regards,
Jinsong