Hello,
I am trying to solve the following algebraic differential equation :
dy1 = h - dy3
dy2 = g - dy4
y3 = K1*y1/(y3+y4)
y4 = K2*y2/(y3+y4)
I tried using the function daspk, but it gives me the following error :
> out <- daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0)
daspk-- warning.. At T(=R1) and stepsize H (=R2) the
nonlinear solver failed to converge
repeatedly of with abs (H) = HMIN &g, 0
Warning messages:
1: In daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0) :
repeated convergence test failures on a step - inaccurate Jacobian or
preconditioner?
2: In daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0) :
Returning early. Results are accurate, as far as they go
> head(out)
time 1 2 3 4 5 6 7 8 9
[1,] 0 0.9 1.33 0 0 10 20 1e-04 0.001 0.001
[2,] 0 0.9 1.33 0 0 10 20 1e-04 0.001 0.001
I read the documentation but it only says that daspk can only solve DAE of
index one maximum, which I do not understand how to determine. I was wondering
if I had to use the function radau but I do not get how to do the M matrix? I
did not managed to get it with the example.
Could it be an error in my program? Here it is :
liquide <- function(z, C, dC, parameters) {
with(as.list(c(C,dC,parameters)),{
Ct = C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + C[7] + C[8] + C[9] + Ceau
V <- 1
K32 <- 6.54*10^2 # m^3/mol
K33 <- 1.37*10^-2 # m^3/mol
K34 <- 0.330 # sans unité
K35 <- 5.81*10^1 # sans unité
kf2 <- 1.37 # m^3/mol
kf3 <- 1.37 # m^3/mol
kf4 <- 8.68*10^-1 # m^3
kf5 <- 157.2*10^-6 # m^3
K2 <- 10^1.44*10^3 # mol/m^3
K3 <- 10^(-3.224)*10^3 # mol/m^3
Ke <- 10^-11 # mol/m^3
r1 <- kf4*C[4]/V - kf4/K34*C[5]^2/(Ct*V)
r2 <- kf3*C[1]*C[2] - kf3/K33*C[4]
r3 <- kf2*C[1]^2 - kf2/K32*C[3]
r4 <- kf5*C[3]/V - kf5/K35*C[5]*C[6]*C[8]/(Ct^2*V)
res1 <- dC[1] + r2 + r3 # dNO2/dt
res2 <- dC[2] + r2 # dNO/dt
res3 <- dC[3] - r3 + r4 # dN2O4/dt
res4 <- dC[4] - r2 + r1 # dN2O3/dt
res5 <- dC[5] -2*r1 - r4 + dC[7] # dHNO2/dt
res6 <- dC[6] - r4 + dC[8] # dHNO3/dt
res7 <- C[7] - K3*C[5]/C[9] # dNO2-/dt
res8 <- C[8] - K2*C[6]/C[9] # dNO3-/dt
res9 <- dC[9] - dC[8] - dC[7] # dCH/dt
list(c(res1, res2, res3, res4, res5, res6, res7, res8, res9))
})
}
yini <- c(0.9, 1.33, 0, 0, 10, 20, 0.001, 22.9, 23)
dyini <- c(1,1,1,1,1,1,1,1,1)
Qm <- 4 # kg/h
Ceau <- (Qm - yini[1]*0.046 - yini[2]*0.03 + yini[3]*0.092 + yini[4]*0.076 +
yini[5]*0.062 + yini[6]*0.063)/0.018 # mol/m^3
parameters <- c(Qm = Qm,
Ceau = Ceau)
liquide(20,yini, dyini,parameters)
z <- seq(0, 120, by = 1) # en seconde
library(deSolve)
out <- daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0)
head(out)
plot(out)
Thanks
Raphaëlle Carraud
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
