Hi R help,
I wanted to simulate two pool model (A&B) using deSolve package for time 0 to
12 by 1. Initial values of the state variables are A=5, B=3. The fluxes are as
follows1) Flux into A= 5 units per unit time 2) Flux from A to B= 0.33) Flux
out of A=0.1 4) Flux from B to A=0.35) Flux out of B=0.3
Here is the R code I compiled to estimate the size of A and B and graph the
output
library(deSolve)
# Define time sequence from 0 to 12 by 1
time <- seq(0,12, by=1)
# Define the function
Mod <- function (t, parms){ derivs <- function(t, state, parms){
with(as.list (c(state, parms)), { #Fluxes inA <- kinA
AtoB <- kAtoB*A Aout <- kAout*A inB <- kinB BtoA <- kBtoA*B
Bout <- kBout*B # Rate of change dA <- inA+BtoA-AtoB-Aout
dB <- inB+AtoB-BtoA-Bout return (list (c(dA, dB))) }) }
#Step 4: Define some starting values for the pools state <- c(A=5, B=3)
return (ode(y=state, times=t,func=derivs, parms=parms, method="rk4"))}
# Starting values of fluxes/parameterspars <- c(kinA=5, kAtoB=0.3, kAout=0.1,
kinB=2, kBtoA=0.3, kBout=0.3)
#Model results
Mod(time, pars)data <- data.frame(Mod(time, pars))plot(Mod(time, pars))
The result does not look right. I could not figure out where I made a mistake
in compiling code.
Any help is highly appreciated.
Wodaj
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
