Differential problem
Ok, now it's good (Pt was in my workplace so it worked for me, I am not used to R using these value to make the program run so I hadn't looked...)
reaction<-function(z, state, dval, parameters) {
with(as.list(c(state,dval,parameters)),{
# rate of change
Tr <- 273+90
P <- 0.98*10^5
K2 <- 10^(2993/Tr-9.226)*(10^-3)
K3 <- 10^(2072/Tr-7.234)*(10^-3)
K4 <- 10^(-20.83/Tr-0.5012)
K5 <- 10^(-965.5/Tr-1.481)
k1 <- (10^(652.1/Tr-0.7356))*(8.314*Tr/P)^2
kf2 <- 1.4*10^-33*(Tr/300)^(-3.8)*6.022*10^23*10^-6
kb2 <- kf2/K2*P/(8.314*Tr)
kf3 <- 3.1*10^-34*(Tr/300)^(-7.7)*10^(-6)*6.022*10^23
kb3 <- kf3/K3*P/(8.314*Tr)
kf4 <- 41
kf5 <- 0.25
r1 <- k1*A^2*H
r4 <- kf4*D*G - kf4/K4*E^2
r5 <- kf5*C*G - kf5/K5*E*I
res1 <- -dA + dB + 2*dC - 2*r1 - 2*r5 #
res2 <- dA + dD + r1 + r4 #
res3 <- K2 - C/B^2 #
res4 <- K3 - D/(A*B) #
res5 <- r5 + 2*r4 - dE #dHNO2/dz
res6 <- r5 -dI #dHNO3/dz
res7 <- -r5 - r4 - dG #dH2O/dz
res8 <- -r1/2 - dH #dO2/dz
list(c(res1, res2, res3, res4, res5, res6, res7, res8))
}) # end with(as.list ...
}
xi <- c(0.3, #x_NO
0.1, #x_NO2
0, #x_N2O4
0, #x_N2O3
0.05, #x_HNO2
0.05, #x_HNO3
0.2, #x_H2O
0.3) #x_O2
Pt <- 0.98*10^5
state <- c(A = xi[1]*Pt,
B = xi[2]*Pt,
C = xi[3]*Pt,
D = xi[4]*Pt,
E = xi[5]*Pt,
I = xi[6]*Pt,
G = xi[7]*Pt,
H = xi[8]*Pt)
dval <- c(dA = 1,
dB = 1,
dC = 0.5,
dD = 0.2,
dE = 0,
dI = 0,
dG = 0,
dH = 0)
parameters <- c(Pt = 0.98*10^5)
z <- seq(0, 1, by = 0.01) # en seconde
library(deSolve)
#out <- ode(y = state, times = z, func = reaction, parameters)
out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0)
head(out)
plot(out)
-----Message d'origine-----
De?: Berend Hasselman [mailto:bhh at xs4all.nl]
Envoy??: jeudi 11 juillet 2013 14:13
??: Rapha?lle Carraud
Cc?: r-help at r-project.org
Objet?: Re: [R] Differential problem
On 11-07-2013, at 13:53, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Sorry for the bug, I had eliminated some lines to avoid making the program too big. Here is the version that works :
reaction<-function(z, state, dval, parameters) {
with(as.list(c(state)),{
# rate of change
Tr <- 273+90
P <- 0.98*10^5
K2 <- 10^(2993/Tr-9.226)*(10^-3)
K3 <- 10^(2072/Tr-7.234)*(10^-3)
K4 <- 10^(-20.83/Tr-0.5012)
K5 <- 10^(-965.5/Tr-1.481)
k1 <- (10^(652.1/Tr-0.7356))*(8.314*Tr/P)^2
kf2 <- 1.4*10^-33*(Tr/300)^(-3.8)*6.022*10^23*10^-6
kb2 <- kf2/K2*P/(8.314*Tr)
kf3 <- 3.1*10^-34*(Tr/300)^(-7.7)*10^(-6)*6.022*10^23
kb3 <- kf3/K3*P/(8.314*Tr)
kf4 <- 41
kf5 <- 0.25
r1 <- k1*A^2*H
r4 <- kf4*D*G - kf4/K4*E^2
r5 <- kf5*C*G - kf5/K5*E*I
res1 <- -dA + dB + 2*dC - 2*r1 - 2*r5 #
res2 <- dA + dD + r1 + r4 #
res3 <- K2 - C/B^2 #
res4 <- K3 - D/(A*B) #
res5 <- r5 + 2*r4 - dE #dHNO2/dz
res6 <- r5 -dI #dHNO3/dz
res7 <- -r5 - r4 - dG #dH2O/dz
res8 <- -r1/2 - dH #dO2/dz
list(c(res1, res2, res3, res4, res5, res6, res7, res8))
}) # end with(as.list ...
}
xi <- c(0.3, #x_NO
0.1, #x_NO2
0, #x_N2O4
0, #x_N2O3
0.05, #x_HNO2
0.05, #x_HNO3
0.2, #x_H2O
0.3) #x_O2
state <- c(A = xi[1]*Pt,
B = xi[2]*Pt,
C = xi[3]*Pt,
D = xi[4]*Pt,
E = xi[5]*Pt,
I = xi[6]*Pt,
G = xi[7]*Pt,
H = xi[8]*Pt)
dval <- c(dA = 1,
dB = 1,
dC = 0.5,
dD = 0.2,
dE = 0,
dI = 0,
dG = 0,
dH = 0)
parameters <- c(Pt = 0.98*10^5)
Doesn't run. Since variable Pt is not defined when you calculate vector state. So define Pt <- .. before xi as in the original example. In the function reaction isn't variable P just Pt from the parameter vector? If so then either do P <- Pt or just use Pt directly (but see next remark).
z <- seq(0, 1, by = 0.01) # en seconde library(deSolve) #out <- ode(y = state, times = z, func = reaction, parameters) out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0) head(out) plot(out) I obtain the following message:
library(deSolve) #out <- ode(y = state, times = z, func = reaction, parameters) out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0)
Error in eval(expr, envir, enclos) : object 'dA' not found.
Obviously since dA isn't defined outside of dval when the functions was defined. And do parms=parameters if you want to use Pt (as I told you in the previous post).
I tried adding the dval and parameters as you said:
with(as.list(c(state,dval,parameters)),{
I get the following message:
Warning messages:
1: In daspk(y = state, dy = dval, times = z, res = reaction, parms = 0) :
matrix of partial derivatives is singular with direct method-some
equations redundant
2: In daspk(y = state, dy = dval, times = z, res = reaction, parms = 0) :
Returning early. Results are accurate, as far as they go
These are warning messages. You get plots. So now is the time to start looking at initial values etc. Since I know next to nothing about DAE's you are on your own here unless someone else comes up with suggestions.
For the calling of the daspk function, I followed the documentation, where you have the same inversion:
daefun <- function(t, y, dy, parameters) {
+ res1 <- dy[1] + y[1] - y[2]
+ res2 <- y[2] * y[1] - t
+
+ list(c(res1, res2))
+ }
library(deSolve) yini <- c(1, 0) dyini <- c(1, 0) times <- seq(0, 10, 0.1) ## solver system.time(out <- daspk(y = yini, dy = dyini, times = times, res = daefun, parms = 0))
Is it wrong? When I modify the order, I obtain again that object dA is not found, so I guessed the doc was right.
Of course see above. Berend