#------------------------------------------------------
# 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
#------------------------------------------------------
# Dissociation constant
K = 1
# parameters
pars <- c(
ka = 1e6, # forward rate
r = 1,
prod = 0.1)
#-------------------------------------
# Chemical problem formulation 1: ODE
#-------------------------------------
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))
})
}
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))
#------------------------------------------------------
# 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))
#------------------------------------------------------
# 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 modelRun the code above in your browser using DataLab