## =======================================================================
## Coupled chemical reactions including an equilibrium
## modeled as (1) an ODE and (2) as a DAE
##
## The model describes three chemical species A,B,D:
## subjected to equilibrium reaction D <- > A + B
## D is produced at a constant rate, prod
## B is consumed at 1s-t order rate, r
## Chemical problem formulation 1: ODE
## =======================================================================
## Dissociation constant
K <- 1
## parameters
pars <- c(
ka = 1e6, # forward rate
r = 1,
prod = 0.1)
Fun_ODE <- function (t, y, pars)
{
with (as.list(c(y, pars)), {
ra <- ka*D # forward rate
rb <- ka/K *A*B # backward rate
## rates of changes
dD <- -ra + rb + prod
dA <- ra - rb
dB <- ra - rb - r*B
return(list(dy = c(dA, dB, dD),
CONC = A+B+D))
})
}
## =======================================================================
## Chemical problem formulation 2: DAE
## 1. get rid of the fast reactions ra and rb by taking
## linear combinations : dD+dA = prod (res1) and
## dB-dA = -r*B (res2)
## 2. In addition, the equilibrium condition (eq) reads:
## as ra = rb : ka*D = ka/K*A*B = > K*D = A*B
## =======================================================================
Res_DAE <- function (t, y, yprime, pars)
{
with (as.list(c(y, yprime, pars)), {
## residuals of lumped rates of changes
res1 <- -dD - dA + prod
res2 <- -dB + dA - r*B
## and the equilibrium equation
eq <- K*D - A*B
return(list(c(res1, res2, eq),
CONC = A+B+D))
})
}
## =======================================================================
## Chemical problem formulation 3: Mass * Func
## Based on the DAE formulation
## =======================================================================
Mass_FUN <- function (t, y, pars)
{
with (as.list(c(y, pars)), {
## as above, but without the
f1 <- prod
f2 <- - r*B
## and the equilibrium equation
f3 <- K*D - A*B
return(list(c(f1, f2, f3),
CONC = A+B+D))
})
}
Mass <- matrix(nr=3, nc=3, byrow = TRUE,
data=c(1, 0, 1, # dA + 0 + dB
-1, 1, 0, # -dA + dB +0
0, 0, 0)) # algebraic
times <- seq(0, 100, by = 2)
## Initial conc; D is in equilibrium with A,B
y <- c(A = 2, B = 3, D = 2*3/K)
## ODE model solved with daspk
ODE <- as.data.frame(daspk(y = y, times = times, func = Fun_ODE,
parms = pars, atol = 1e-10, rtol = 1e-10))
## Initial rate of change
dy <- c(dA = 0, dB = 0, dD = 0)
## DAE model solved with daspk
DAE <- as.data.frame(daspk(y = y, dy = dy, times = times,
res = Res_DAE, parms = pars, atol = 1e-10, rtol = 1e-10))
MASS<- daspk(y=y, times=times, func = Mass_FUN, parms = pars, mass = Mass)
## ================
## plotting output
## ================
opa <- par(mfrow = c(2,2))
for (i in 2:5)
{
plot(ODE$time,ODE[,i],xlab = "time",
ylab = "conc",main = names(ODE)[i],type = "l")
points(DAE$time,DAE[,i],col = "red")
}
legend("bottomright",lty = c(1,NA),pch = c(NA,1),
col = c("black","red"),legend = c("ODE","DAE"))
# difference between both implementations:
max(abs(ODE-DAE))
par(mfrow = opa)
## =======================================================================
## same DAE model, now with the Jacobian
## =======================================================================
jacres_DAE <- function (t, y, yprime, pars, cj)
{
with (as.list(c(y, yprime, pars)), {
## res1 = -dD - dA + prod
PD[1,1] <- -1*cj # d(res1)/d(A)-cj*d(res1)/d(dA)
PD[1,2] <- 0 # d(res1)/d(B)-cj*d(res1)/d(dB)
PD[1,3] <- -1*cj # d(res1)/d(D)-cj*d(res1)/d(dD)
## res2 = -dB + dA - r*B
PD[2,1] <- 1*cj
PD[2,2] <- -r -1*cj
PD[2,3] <- 0
## eq = K*D - A*B
PD[3,1] <- -B
PD[3,2] <- -A
PD[3,3] <- K
return(PD)
})
}
PD <- matrix(nc = 3, nr = 3, 0)
DAE2 <- as.data.frame(daspk(y = y, dy = dy, times = times,
res = Res_DAE, jacres = jacres_DAE, jactype = "fullusr",
parms = pars, atol = 1e-10, rtol = 1e-10))
max(abs(DAE-DAE2))
## See \dynload subdirectory for a FORTRAN implementation of this model
## =======================================================================
## The chemical model as a DLL, with production a forcing function
## =======================================================================
times <- seq(0, 100, by = 2)
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),0.1*(1+runif(11)*1)))
ODE_dll <- as.data.frame(daspk(y=y,dy=dy,times=times,res="chemres",
dllname="deSolve", initfunc="initparms",
initforc="initforcs", parms=pars, forcings=prod,
atol=1e-10,rtol=1e-10,nout=2, outnames=c("CONC","Prod")))
opa <- par(mfrow = c(1,2))
plot(ODE_dll$time,ODE_dll$Prod,xlab = "time",
ylab = "/day",main = "production rate",type = "l")
plot(ODE_dll$time,ODE_dll$D,xlab = "time",
ylab = "conc",main = "D",type = "l")
par(mfrow = opa)Run the code above in your browser using DataLab