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

[Vince Pile]

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