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

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

--
Daniel Kaschek
Institute of Physics
Freiburg University

Room 210
Phone: +49 761 2038531