Skip to content

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

5 messages · Daniel Kaschek, Vince Pile, Karline Soetaert

Original
Vince Pile
# 
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")
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      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:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      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:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

  
  


  


      

    

    

      
      

        


  

        

      Vince Pile
    

    

      
      #
    

  


  

      
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