[R-sig-dyn-mod] Which parameter identifies the 'time' to reach steady state with runsteady?

[Daniel Kaschek]

Dear Vince,

runsteady returns the steady state values. I think, it is not supposed
to return the time when steady state is "reached". The problem about
"reached" is that this will never happen. However, you can ask when the
solution is less than e.g. 5% away from the steady state solution.

I modified the last part of your code to plot the solution divided by
the steady state. In the plot you see that OM reaches the steady state
after approximately t=30 whereas O2 ans HS follow at t=50.

You can determine this time numerically by the following line (see code
below for the definition of plotOut):

apply(plotOut, 2, function(col) times[abs(1-col)<0.05][1])

Hope, this helps.
Cheers,
Daniel
 

# time
times = seq(0, 100, 0.05)

# dynamic soln
dynOut <- lsoda(y = y, func = model, time = times, parms = pars)

# steady state soln
ssOut <- runsteady(y = y, func = model, parms = pars, times = c(0,
1000))

# plot
head(dynOut)
plotOut <- t(t(dynOut[,c("OM", "O2", "HS")])/ssOut$y[c("OM", "O2",
"HS")])


matplot(times, plotOut, type="l", lty=1:4, col=1:4, lwd=3, ylim=c(0.9,
1.1))
abline(h = c(0.95, 1.05), lty=2, col="gray")
legend("topleft", legend=c("OM", "O2", "HS"), lty=1:4, col=1:4, lwd=3)


On Mo, 2015-02-09 at 21:51 -0500, Vince Pile wrote:
> Hello,
> New to working with deSolve and rootSolve and ran the example under
> runsteady and tried to determine the time to reach steady state but despite
> all my attempts to manipulate the 'istate' or 'rstate' values to match the
> plot output could not work it out.
> Thanks in advance for your help!
> Vince
> ## =============================================
> ## A simple sediment biogeochemical model
> ## =============================================
> 
> model<-function(t, y, pars) {
> 
>   with (as.list(c(y, pars)),{
> 
>     Min       = r*OM
>     oxicmin   = Min*(O2/(O2+ks))
>     anoxicmin = Min*(1-O2/(O2+ks))* SO4/(SO4+ks2)
> 
>     dOM  = Flux - oxicmin - anoxicmin
>     dO2  = -oxicmin      -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2)
>     dSO4 = -0.5*anoxicmin  +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4)
>     dHS  = 0.5*anoxicmin   -rox*HS*(O2/(O2+ks)) + D*(BHS-HS)
> 
>     list(c(dOM, dO2, dSO4, dHS), SumS = SO4+HS)
>   })
> }
> 
> # parameter values
> pars <- c(D = 1, Flux = 100, r = 0.1, rox = 1,
>           ks = 1, ks2 = 1, BO2 = 100, BSO4 = 10000, BHS = 0)
> # initial conditions
> y <- c(OM = 1, O2 = 1, SO4 = 1, HS = 1)
> 
> # time
> times = seq(0, 40, 0.05)
> 
> # dynamic soln
> dynOut <- lsoda(y = y, func = model, time = times, parms = pars, verbose =
> TRUE)
> 
> # steady state soln
> ssOut <- runsteady(y = y, func = model, parms = pars, times = c(0, 1000),
> verbose=TRUE)
> 
> # plot
> head(dynOut)
> matplot(dynOut[,1], dynOut[,c("OM", "O2", "HS")], type="l", lty=1:4,
> col=1:4, lwd=3)
> legend("topleft", legend=c("OM", "O2", "HS"), lty=1:4, col=1:4, lwd=3)
> abline(h=seq(0,900,100), v=seq(0,40,5), lty=2, col="grey")
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

-- 
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531