Skip to content

[R-sig-dyn-mod] Help with function

2 messages · Austin Mullings, Daniel Kaschek

Original
Austin Mullings
# 
I am trying to get this equation to work with deSolve, and can't seem to
get it to work. Does someone know how one would write the function for this?

dY1(t)/dt = a - bY3(t) Y4(t) - cY1(t),
dY2(t)/dt = Y5(t) [dY1(t)/dt] [1 + Y6(t)] - eY2(t) - {f Y3(t) Y4(t) Y1(t) -
g},
dY3(t)/dt = h - iY2(t),
dY4(t)/dt = Y6(t) [dY1(t)/dt] - jY4(t) + k {f Y3(t) Y4(t) Y1(t) - g} + d,
dY5(t)/dt = 1.0 - Y5(t) [Y1(t) + Y4(t) + 1],
dY6(t)/dt = 1.0 - Y6(t) [Y1(t) + 1].
Daniel Kaschek
# 
Hi Austin,

one possibility to do this is the code below. In this example code, I
randomly initialize your parameters and the initial state. Instead of
using the indexes to access parameter values or states you can use the
names as well (see ?ode, first example). However, in my experience, the
function with() slows down your code considerably. 

Best wishes,
Daniel


library(deSolve)

myfn <- function(t, y, p) {
  
  dy <- numeric(6)
  dy[1] <- p[1] - p[2]*y[4] - p[3]*y[1]
  dy[2] <- y[5]*dy[1]*(1+y[6]) - p[5]*y[2] - (p[6]*y[3]*y[4]*y[1] - p[7])
  dy[3] <- p[8] - p[9]*y[2]
  dy[4] <- y[6]*dy[1] - p[10]*y[4] + p[11]*(p[6]*y[3]*y[4]*y[1] - p[7]) + p[4]
  dy[5] <- 1.0 - y[5]*(y[1] + y[4] + 1)
  dy[6] <- 1.0 - y[6]*(y[1] + 1)
  
  return(list(dy))
  
}

pars <- runif(11)
yini <- runif(6)
times <- seq(0, 10, .01)


out <- ode(yini, times, myfn, pars)
plot(out)
On Mo, 2015-02-16 at 13:15 -0700, Austin Mullings wrote: