Hi,
I've solved a simple differential equation describing the degradation of
amino acid carbon (THAA-C) using deSolve.
Code is a follows:
# Input of model parameters, a and b describes form of curve, i is apparent
initial age of Org. C.
parameters <- c(a = a, b = b, i=i)
# Initial value of the model, G
state = c(G = G)
#specifies the function degradation as a function of time (t), inital state
(state) and parameters (given parameters in model)
degradation = function (t, state, parameters) {
#Specifies that state an parameters is inputted as lists, hence order in the
lists is important
with(as.list(c(state, parameters)),
#differential equation to be solved
{dG = (-a*(t+i)^b)*(G)
list(c(dG))
})
}
# Makes sequence of times to be solved for. Start at t0 and proceed to iteno
in steps of 1 day
times = seq(t0, iteno, by = 1)
#calls solver, to solve diffrential equation and saves results as variable
THAAC
THAAC = ode (y = state, times = times, func = degradation, parms parameters)
This all works fine, and output THAAC contains carbon concentrations of G at
all times t. However, i would like the model to give me the rates, i.e. dG
for all times. And possibly also the value of -a*(t+i)^b for all times.
Does anyone know if this is possible?
Kristoffer Piil
--
View this message in context:
http://r.789695.n4.nabble.com/deSolve-output-tp3738970p3738970.html
Sent from the R help mailing list archive at Nabble.com.
Hi,
Try with the following ODE function. This should give you an extra
column with the derivative of G in your THAAC matrix.
degradation = function (t, state, parameters) {
?with(as.list(c(state, parameters)),
? ?{dG = (-a*(t+i)^b)*(G)
? ?list(c(dG),dG=dG)
? ?})
}
Any additional variables that you want to output from your ODE system
function need to be a separate level of the output list. For instance:
# Example taken from the deSolve vignette
parameters <- c(a = -8/3,b = -10, c = 28)
state <- c(X = 1,Y = 1, Z = 1)
Lorenz<-function(t, state, parameters) {
?with(as.list(c(state, parameters)),{
?# rate of change
?dX <- a*X + Y*Z
?dY <- b * (Y-Z)
?dZ <- -X*Y + c*Y - Z
?# return the rate of change
?list(c(dX, dY, dZ),dX=dX, dY=dY,Dummy=-X/Y)
?}) # end with(as.list ...
}
times <- seq(0, 100, by = 0.01)
library(deSolve)
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
head(out)
Sebastien