R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("odesolve-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('odesolve')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ccl4data"
> 
> ### * ccl4data
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ccl4data
> ### Title: Closed chamber study of CCl4 metabolism by rats.
> ### Aliases: ccl4data
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(ccl4data)
> 
> 
> 
> cleanEx(); ..nameEx <- "ccl4data.avg"
> 
> ### * ccl4data.avg
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ccl4data.avg
> ### Title: Closed chamber study of CCl4 metabolism by rats.
> ### Aliases: ccl4data.avg
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(ccl4data.avg)
> 
> 
> 
> cleanEx(); ..nameEx <- "lsoda"
> 
> ### * lsoda
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lsoda
> ### Title: Solve System of ODE (ordinary differential equation)s.
> ### Aliases: lsoda
> ### Keywords: math
> 
> ### ** Examples
> 
> ### lsexamp -- example from lsoda source code
> 
> ## names makes this easier to read, but may slow down execution.
> parms <- c(k1=0.04, k2=1e4, k3=3e7)
> my.atol <- c(1e-6,  1e-10,  1e-6)
> times <- c(0,4 * 10^(-1:10))
> lsexamp <- function(t, y, p)
+   {
+     yd1 <- -p["k1"] * y[1] + p["k2"] * y[2]*y[3]
+     yd3 <- p["k3"] * y[2]^2
+     list(c(yd1,-yd1-yd3,yd3),c(massbalance=sum(y)))
+   }
> exampjac <- function(t, y, p)
+   {
+     c(-p["k1"],  p["k1"],  0,
+   
+         p["k2"]*y[3],
+       - p["k2"]*y[3] - 2*p["k3"]*y[2],
+                        2*p["k3"]*y[2],
+   
+       p["k2"]*y[2],  -p["k2"]*y[2],  0
+       )
+   }
>   
> require(odesolve)
[1] TRUE
> ## measure speed (here and below)
> system.time( 
+ out <- lsoda(c(1,0,0),times,lsexamp, parms, rtol=1e-4, atol= my.atol)
+ )  
[1] 0.05 0.00 0.05 0.00 0.00
> out
       time            1            2          3 massbalance
 [1,] 0e+00 1.000000e+00 0.000000e+00 0.00000000           1
 [2,] 4e-01 9.851712e-01 3.386380e-05 0.01479493           1
 [3,] 4e+00 9.055333e-01 2.240655e-05 0.09444430           1
 [4,] 4e+01 7.158403e-01 9.186334e-06 0.28415047           1
 [5,] 4e+02 4.505250e-01 3.222964e-06 0.54947175           1
 [6,] 4e+03 1.831976e-01 8.941774e-07 0.81680155           1
 [7,] 4e+04 3.898729e-02 1.621940e-07 0.96101254           1
 [8,] 4e+05 4.936362e-03 1.984221e-08 0.99506362           1
 [9,] 4e+06 5.161832e-04 2.065786e-09 0.99948381           1
[10,] 4e+07 5.179810e-05 2.072029e-10 0.99994820           1
[11,] 4e+08 5.283548e-06 2.113430e-11 0.99999472           1
[12,] 4e+09 4.658902e-07 1.863562e-12 0.99999953           1
[13,] 4e+10 1.424583e-08 5.698339e-14 0.99999999           1
attr(,"istate")
[1] 2
> 
> ## This is what the authors of lsoda got for the example:
> 
> ## the output of this program (on a cdc-7600 in single precision)
> ## is as follows..
> ##
> ##   at t =  4.0000e-01   y =  9.851712e-01  3.386380e-05  1.479493e-02
> ##   at t =  4.0000e+00   y =  9.055333e-01  2.240655e-05  9.444430e-02
> ##   at t =  4.0000e+01   y =  7.158403e-01  9.186334e-06  2.841505e-01
> ##   at t =  4.0000e+02   y =  4.505250e-01  3.222964e-06  5.494717e-01
> ##   at t =  4.0000e+03   y =  1.831975e-01  8.941774e-07  8.168016e-01
> ##   at t =  4.0000e+04   y =  3.898730e-02  1.621940e-07  9.610125e-01
> ##   at t =  4.0000e+05   y =  4.936363e-03  1.984221e-08  9.950636e-01
> ##   at t =  4.0000e+06   y =  5.161831e-04  2.065786e-09  9.994838e-01
> ##   at t =  4.0000e+07   y =  5.179817e-05  2.072032e-10  9.999482e-01
> ##   at t =  4.0000e+08   y =  5.283401e-06  2.113371e-11  9.999947e-01
> ##   at t =  4.0000e+09   y =  4.659031e-07  1.863613e-12  9.999995e-01
> ##   at t =  4.0000e+10   y =  1.404280e-08  5.617126e-14  1.000000e+00
> 
> ## Using the analytic jacobian speeds up execution a little :
> 
> system.time( 
+ outJ <- lsoda(c(1,0,0),times,lsexamp, parms, rtol=1e-4, atol= my.atol,
+               jac = exampjac)
+ )
[1] 0.04 0.00 0.04 0.00 0.00
>   
> all.equal(out, outJ) # TRUE
[1] TRUE
> 
> ## Example for using hmax
> 
> ## Parameters for steady state conditions
> parms  <- c(a=0.0, b=0.0, c=0.1, d=0.1, e=0.1, f=0.1, g=0.0)
> 
> ## A simple resource limited Lotka-Volterra-Model
> ## Note passing parameters through using a closure
> lvmodel <- with(as.list(parms), function(t, x, parms) {
+     import <- sigimp(t)
+     ds <- import - b*x["s"]*x["p"] + g*x["k"]
+     dp <- c*x["s"]*x["p"]  - d*x["k"]*x["p"]
+     dk <- e*x["p"]*x["k"]  - f*x["k"]
+     res<-c(ds, dp, dk)
+     list(res)
+   })
> 
> ## vector of timesteps
> times  <- seq(0, 100, length=101)
> 
> ## external signal with rectangle impulse
> signal <- as.data.frame(list(times = times,
+                             import = rep(0,length(times))))
> 
> signal$import[signal$times >= 10 & signal$times <=11] <- 0.2
> 
> sigimp <- approxfun(signal$times, signal$import, rule=2)
> 
> ## Start values for steady state
> y<-xstart <- c(s=1, p=1, k=1)
> 
> ## LSODA (default step size)
> out2 <- as.data.frame(lsoda(xstart, times, lvmodel, parms))
> 
> ## LSODA: with fixed maximum time step
> out3 <- as.data.frame(lsoda(xstart, times, lvmodel, parms, hmax=1))
> 
> par(mfrow=c(2,2))
> plot (out2$time, out2$s,  type="l",   ylim=c(0,3))
> lines(out3$time, out3$s, col="green", lty="dotted")
> 
> plot (out2$time, out2$p, type="l",    ylim=c(0,3))
> lines(out3$time, out3$p, col="green", lty="dotted")
> 
> plot (out2$time, out2$k, type="l",    ylim=c(0,3))
> lines(out3$time, out3$k, col="green", lty="dotted")
> 
> plot (out2$p, out2$k, type="l", ylim=range(out2$k,out3$k))
> lines(out3$p, out3$k, col="green", lty="dotted")
> 
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "rk4"
> 
> ### * rk4
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rk4
> ### Title: Solve System of ODE (ordinary differential equation)s by
> ###   classical Runge-Kutta 4th order integration.
> ### Aliases: rk4
> ### Keywords: math
> 
> ### ** Examples
> 
> ## A simple resource limited Lotka-Volterra-Model
> lvmodel <- function(t, x, parms) {
+   s <- x[1] # substrate
+   p <- x[2] # producer
+   k <- x[3] # consumer
+   with(as.list(parms),{
+     import <- approx(signal$times, signal$import, t)$y
+     ds <- import - b*s*p + g*k
+     dp <- c*s*p  - d*k*p
+     dk <- e*p*k  - f*k
+     res<-c(ds, dp, dk)
+     list(res)
+   })
+ }
> 
> ## vector of timesteps
> times  <- seq(0, 100, length=101)
> 
> ## external signal with rectangle impulse
> signal <- as.data.frame(list(times = times,
+                             import = rep(0,length(times))))
> 
> signal$import[signal$times >= 10 & signal$times <=11] <- 0.2
> 
> ## Parameters for steady state conditions
> parms  <- c(a=0.0, b=0.0, c=0.1, d=0.1, e=0.1, f=0.1, g=0.0)
> 
> ## Start values for steady state
> y<-xstart <- c(s=1, p=1, k=1)
> 
> ## Classical RK4 with fixed time step
> out1  <- as.data.frame(rk4(xstart, times, lvmodel, parms))
> 
> ## LSODA (default step size)
> out2 <- as.data.frame(lsoda(xstart, times, lvmodel, parms))
> 
> ## LSODA: with fixed maximum time step
> out3 <- as.data.frame(lsoda(xstart, times, lvmodel, parms, hmax=1))
> 
> par(mfrow=c(2,2))
> plot (out1$time, out1$s,  type="l",   ylim=c(0,3))
> lines(out2$time, out2$s, col="red",   lty="dotted")
> lines(out3$time, out3$s, col="green", lty="dotted")
> 
> plot (out1$time, out1$p, type="l",    ylim=c(0,3))
> lines(out2$time, out2$p, col="red",   lty="dotted")
> lines(out3$time, out3$p, col="green", lty="dotted")
> 
> plot (out1$time, out1$k, type="l",    ylim=c(0,3))
> lines(out2$time, out2$k, col="red",   lty="dotted")
> lines(out3$time, out3$k, col="green", lty="dotted")
> 
> plot (out1$p, out1$k, type="l")
> lines(out2$p, out2$k, col="red",   lty="dotted")
> lines(out3$p, out3$k, col="green", lty="dotted")
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.98 0.02 1.01 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')