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

Vince Pile vpileggi10 at gmail.com
Tue Feb 10 20:09:46 CET 2015


Great! now I have 2 ways and can compare...Thanks very much Karline!... and
of course a great big thank you to you and your colleagues for deSolve and
rootSolve:)

On Tue, Feb 10, 2015 at 11:11 AM, Karline Soetaert <Karline.Soetaert at nioz.nl
> wrote:

> 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
>
> _______________________________________________
> 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]]



More information about the R-sig-dynamic-models mailing list