inst/doc/compiledCode.R

In deSolve: Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')

### R code from vignette source 'compiledCode.Rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
library("deSolve")
options(prompt = "R> ")
options(width=70)


###################################################
### code chunk number 2: the_Rmodel
###################################################
model <- function(t, Y, parameters) {
  with (as.list(parameters),{

    dy1 = -k1*Y[1] + k2*Y[2]*Y[3]
    dy3 = k3*Y[2]*Y[2]
    dy2 = -dy1 - dy3

    list(c(dy1, dy2, dy3))
  })
}


###################################################
### code chunk number 3: Jacobian_in_R
###################################################
jac <- function (t, Y, parameters) {
  with (as.list(parameters),{

    PD[1,1] <- -k1
    PD[1,2] <- k2*Y[3]
    PD[1,3] <- k2*Y[2]
    PD[2,1] <- k1
    PD[2,3] <- -PD[1,3]
    PD[3,2] <- k3*Y[2]
    PD[2,2] <- -PD[1,2] - PD[3,2]

    return(PD)
  })
}


###################################################
### code chunk number 4: Run_Rmodel
###################################################
parms <- c(k1 = 0.04, k2 = 1e4, k3=3e7)
Y     <- c(1.0, 0.0, 0.0)
times <- c(0, 0.4*10^(0:11))
PD    <- matrix(nrow = 3, ncol = 3, data = 0)
out   <- ode(Y, times, model, parms = parms, jacfunc = jac)


###################################################
### code chunk number 5: compile_DLLmodel_F (eval = FALSE)
###################################################
##   system("R CMD SHLIB mymod.f")


###################################################
### code chunk number 6: compile_DLLmodel_C (eval = FALSE)
###################################################
##   system("R CMD SHLIB mymod.c")


###################################################
### code chunk number 7: compiledCode.Rnw:724-766
###################################################
caraxisfun <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    yb <- r * sin(w * t)
    xb <- sqrt(L * L - yb * yb)
    Ll <- sqrt(xl^2 + yl^2)
    Lr <- sqrt((xr - xb)^2 + (yr - yb)^2)

    dxl <- ul; dyl <- vl; dxr <- ur; dyr <- vr

    dul  <- (L0-Ll) * xl/Ll      + 2 * lam2 * (xl-xr) + lam1*xb
    dvl  <- (L0-Ll) * yl/Ll      + 2 * lam2 * (yl-yr) + lam1*yb - k * g

    dur  <- (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
    dvr  <- (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k * g

    c1   <- xb * xl + yb * yl
    c2   <- (xl - xr)^2 + (yl - yr)^2 - L * L

    list(c(dxl, dyl, dxr, dyr, dul, dvl, dur, dvr, c1, c2))
  })
}
eps <- 0.01; M <- 10; k <- M * eps^2/2;
L <- 1; L0 <- 0.5; r <- 0.1; w <- 10; g <- 1

parameter <- c(eps = eps, M = M, k = k, L = L, L0 = L0,
               r = r, w = w, g = g)

yini <- c(xl = 0, yl = L0, xr = L, yr = L0, ul = -L0/L, vl = 0,
          ur = -L0/L, vr = 0, lam1 = 0, lam2 = 0)

# the mass matrix
Mass      <- diag(nrow = 10, 1)
Mass[5,5] <- Mass[6,6] <- Mass[7,7] <- Mass[8,8] <- M * eps * eps/2
Mass[9,9] <- Mass[10,10] <- 0
Mass

# index of the variables: 4 of index 1, 4 of index 2, 2 of index 3
index <- c(4, 4, 2)

times <- seq(0, 3, by = 0.01)
out <- radau(y = yini, mass = Mass, times = times, func = caraxisfun,
        parms = parameter, nind = index)


###################################################
### code chunk number 8: caraxis
###################################################
plot(out, which = 1:4, type = "l", lwd = 2)


###################################################
### code chunk number 9: figcaraxis
###################################################
plot(out, which = 1:4, type = "l", lwd = 2)


###################################################
### code chunk number 10: compiledCode.Rnw:949-978
###################################################
## the model, 5 state variables
f1 <- function  (t, y, parms) {
  ydot <- vector(len = 5)

  ydot[1] <-  0.1*y[1] -0.2*y[2]
  ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3]
  ydot[3] <-           -0.3*y[2] +0.1*y[3] -0.2*y[4]
  ydot[4] <-                     -0.3*y[3] +0.1*y[4] -0.2*y[5]
  ydot[5] <-                               -0.3*y[4] +0.1*y[5]

  return(list(ydot))
}

## the Jacobian, written in banded form
bandjac <- function  (t, y, parms) {
  jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE,
                data = c( 0  , -0.2, -0.2, -0.2, -0.2,
                          0.1,  0.1,  0.1,  0.1,  0.1,
                         -0.3, -0.3, -0.3, -0.3,    0))
  return(jac)
}

## initial conditions and output times
yini  <- 1:5
times <- 1:20

## stiff method, user-generated banded Jacobian
out  <- lsode(yini, times, f1, parms = 0, jactype = "bandusr",
              jacfunc = bandjac, bandup = 1, banddown = 1)


###################################################
### code chunk number 11: compiledCode.Rnw:1061-1072
###################################################
## Parameter values and initial conditions
Parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32,
        16.17, 281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671,
        0.4027255, 1000, 0.02, 1.0, 3.8)

yini <- c( AI=21, AAM=0, AT=0, AF=0, AL=0, CLT=0,  AM=0  )

## the rate of change
DLLfunc(y = yini, dllname = "deSolve", func = "derivsccl4",
        initfunc = "initccl4", parms = Parms, times = 1,
        nout = 3, outnames = c("DOSE", "MASS", "CP")  )


###################################################
### code chunk number 12: compiledCode.Rnw:1083-1099
###################################################
pars <- c(K = 1, ka   = 1e6, r    = 1)

## Initial conc; D is in equilibrium with A,B
y     <- c(A = 2, B = 3, D = 2*3/pars["K"])

## Initial rate of change
dy  <- c(dA = 0, dB = 0, dD = 0)

## production increases with time
prod <- matrix(nc=2,data=c(seq(0,100,by=10),seq(0.1,0.5,len=11)))

DLLres(y=y,dy=dy,times=5,res="chemres",
          dllname="deSolve", initfunc="initparms",
          initforc="initforcs", parms=pars, forcings=prod,
          nout=2, outnames=c("CONC","Prod"))



###################################################
### code chunk number 13: compiledCode.Rnw:1267-1275
###################################################
Flux <- matrix(ncol=2,byrow=TRUE,data=c(
  1, 0.654, 11, 0.167,   21, 0.060, 41, 0.070, 73,0.277, 83,0.186,
  93,0.140,103, 0.255,  113, 0.231,123, 0.309,133,1.127,143,1.923,
  153,1.091,163,1.001,  173, 1.691,183, 1.404,194,1.226,204,0.767,
  214, 0.893,224,0.737, 234,0.772,244, 0.726,254,0.624,264,0.439,
  274,0.168,284 ,0.280, 294,0.202,304, 0.193,315,0.286,325,0.599,
  335, 1.889,345, 0.996,355,0.681,365,1.135))
head(Flux)


###################################################
### code chunk number 14: compiledCode.Rnw:1280-1281
###################################################
parms <- 0.01


###################################################
### code chunk number 15: compiledCode.Rnw:1287-1289
###################################################
meanDepo <- mean(approx(Flux[,1],Flux[,2], xout=seq(1,365,by=1))$y)
Yini <- c(y=meanDepo/parms)


###################################################
### code chunk number 16: compiledCode.Rnw:1305-1312
###################################################
times <- 1:365
out <- ode(y=Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc="scocforc",  forcings=Flux,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation","Depo"))
head(out)


###################################################
### code chunk number 17: compiledCode.Rnw:1318-1324
###################################################
fcontrol <- list(method="constant")
out2 <- ode(y=Yini, times, func = "scocder",
    parms = parms, dllname = "deSolve",
    initforc="scocforc",  forcings=Flux, fcontrol=fcontrol,
    initfunc = "scocpar", nout = 2,
    outnames = c("Mineralisation","Depo"))


###################################################
### code chunk number 18: scoc
###################################################
par (mfrow=c(1,2))
plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='linear'")
lines(out[,"time"], out[,"Mineralisation"], lwd=2, col="blue")
legend("topleft",lwd=1:2,col=c("red","blue"), c("Flux","Mineralisation"))

plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='constant'")
lines(out2[,"time"], out2[,"Mineralisation"], lwd=2, col="blue")


###################################################
### code chunk number 19: figscoc
###################################################
par (mfrow=c(1,2))
plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='linear'")
lines(out[,"time"], out[,"Mineralisation"], lwd=2, col="blue")
legend("topleft",lwd=1:2,col=c("red","blue"), c("Flux","Mineralisation"))

plot(out, which = "Depo", col="red",
  xlab="days", ylab="mmol C/m2/ d", main="method='constant'")
lines(out2[,"time"], out2[,"Mineralisation"], lwd=2, col="blue")


###################################################
### code chunk number 20: compiledCode.Rnw:1359-1391
###################################################
SPCmod <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    import <- input(t)
    dS <- import - b*S*P + g*C    # substrate
    dP <- c*S*P  - d*C*P          # producer
    dC <- e*P*C  - f*C            # consumer
    res <- c(dS, dP, dC)
    list(res, signal = import)
  })
}

## The parameters
parms  <- c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)

## vector of timesteps
times  <- seq(0, 100, by=0.1)

## external signal with several rectangle impulses
signal <- as.data.frame(list(times = times,
                 import = rep(0, length(times))))

signal$import <- ifelse((trunc(signal$times) %% 2 == 0), 0, 1)
sigimp <- approxfun(signal$times, signal$import, rule = 2)

## Start values for steady state
xstart <- c(S = 1, P = 1, C = 1)

## Solve model
print (system.time(
out <- ode(y = xstart,times = times,
           func = SPCmod, parms, input = sigimp)
))


###################################################
### code chunk number 21: lv
###################################################
plot(out)


###################################################
### code chunk number 22: figlv
###################################################
plot(out)


###################################################
### code chunk number 23: compiledCode.Rnw:1510-1513
###################################################
eventdata <- data.frame(var=rep("C",10),time=seq(10,100,10),value=rep(0.5,10),
  method=rep("multiply",10))
eventdata


###################################################
### code chunk number 24: compiledCode.Rnw:1600-1618
###################################################
derivs <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
    if (t < tau)
      ytau <- c(1, 1)
    else
      ytau <- lagvalue(t - tau, c(1, 2))

    dN <- f * N - g * N * P
    dP <- e * g * ytau[1] * ytau[2] - m * P
    list(c(dN, dP), tau=ytau[1], tau=ytau[2])
  })
}

yinit <- c(N = 1, P = 1)
times <- seq(0, 500)
parms <- c(f = 0.1, g = 0.2, e = 0.1, m = 0.1, tau = .2)
yout <- dede(y = yinit, times = times, func = derivs, parms = parms)
head(yout)