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

[Karline Soetaert]

In addition to that:
You can ask the time at which steady-state is reached via the "attributes" function.
In your case you would need to write:
attributes(ssOut)$time


Karline

-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Vince Pile
Sent: dinsdag 10 februari 2015 13:50
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Which parameter identifies the 'time' to reach steady state with runsteady?

Yes, your insight and addition to the code solved my burning problem and saved me hours of tedious work.
Thanks very much Daniel!!
Double Cheers,
Vince

On Tue, Feb 10, 2015 at 3:57 AM, Daniel Kaschek < daniel.kaschek at physik.uni-freiburg.de> wrote:

> 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
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

	[[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