An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20130711/0a8b3275/attachment.pl>
Differential problem
7 messages · Raphaëlle Carraud, Berend Hasselman
On 11-07-2013, at 09:13, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Hello, I have the following differential equation system to solve, the variables I wish to obtain being A, B, C, D, E, I, G and H. 0 = -dA + dB + 2*dC - 2*r1 - 2*r5 0 = dA + dD + r1 + r4 0 = K2 - C/B^2 0 = K3 - D/(A*B) 0 = r5 + 2*r4 - dE 0 = r5 -dI 0 = -r5 - r4 - dG 0 = -r1/2 - dH r1, r4 and r5 are variables expressed as functions of A, B, C, D, I, G and H, calculated previously in the integrated function. K2 and K3 are parameters. As I have two algebraic equations, I think this system is a DAE (Algebraic differential equation). I found in the package deSolve two functions that resolve DAE but I didn't manage to obtain a solution; it says that the variable dA cannot be found.
Show us your script where you define the function and run the DAE solver. Without that nobody can provide an answer. Berend.
Does anyone know how to solve this problem? Thank you Rapha?lle Carraud Rapha?lle Carraud [[alternative HTML version deleted]]
______________________________________________ R-help at 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.
Sorry,
Here is the program I have until now:
reaction<-function(z, state, dval, parameters) {
with(as.list(c(state)),{
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
res6 <- r5 -dI
res7 <- -r5 - r4 - dG
res8 <- -r1/2 - dH
list(c(res1, res2, res3, res4, res5, res6, res7, res8))
}) # end with(as.list ...
}
Ti <- 273+90 #K
Pt <- 0.98*10^5 #Pa
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)
z <- seq(0, 1, by = 0.01)
library(deSolve)
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 11:18
??: Rapha?lle Carraud
Cc?: r-help at r-project.org
Objet?: Re: [R] Differential problem
On 11-07-2013, at 09:13, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Hello, I have the following differential equation system to solve, the variables I wish to obtain being A, B, C, D, E, I, G and H. 0 = -dA + dB + 2*dC - 2*r1 - 2*r5 0 = dA + dD + r1 + r4 0 = K2 - C/B^2 0 = K3 - D/(A*B) 0 = r5 + 2*r4 - dE 0 = r5 -dI 0 = -r5 - r4 - dG 0 = -r1/2 - dH r1, r4 and r5 are variables expressed as functions of A, B, C, D, I, G and H, calculated previously in the integrated function. K2 and K3 are parameters. As I have two algebraic equations, I think this system is a DAE (Algebraic differential equation). I found in the package deSolve two functions that resolve DAE but I didn't manage to obtain a solution; it says that the variable dA cannot be found.
Show us your script where you define the function and run the DAE solver. Without that nobody can provide an answer. Berend.
Does anyone know how to solve this problem? Thank you Rapha?lle Carraud Rapha?lle Carraud [[alternative HTML version deleted]]
______________________________________________ R-help at 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.
On 11-07-2013, at 12:05, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Sorry,
Here is the program I have until now:
reaction<-function(z, state, dval, parameters) {
with(as.list(c(state)),{
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
res6 <- r5 -dI
res7 <- -r5 - r4 - dG
res8 <- -r1/2 - dH
list(c(res1, res2, res3, res4, res5, res6, res7, res8))
}) # end with(as.list ...
}
Ti <- 273+90 #K
Pt <- 0.98*10^5 #Pa
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)
z <- seq(0, 1, by = 0.01)
library(deSolve)
out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0)
head(out)
plot(out)
When I run this I get this error message:
Error in eval(expr, envir, enclos) : object 'Tr' not found
And shouldn't the first line of the reaction function be this
with(as.list(c(state,dval,parameters)),{
in stead of this
with(as.list(c(state)),{
The call of daspk also seems incorrect; shouldn't it be
out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = parameters)
And finally: where does variable P come from? Not defined anywhere!
Berend
-----Message d'origine----- De : Berend Hasselman [mailto:bhh at xs4all.nl] Envoy? : jeudi 11 juillet 2013 11:18 ? : Rapha?lle Carraud Cc : r-help at r-project.org Objet : Re: [R] Differential problem On 11-07-2013, at 09:13, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Hello, I have the following differential equation system to solve, the variables I wish to obtain being A, B, C, D, E, I, G and H. 0 = -dA + dB + 2*dC - 2*r1 - 2*r5 0 = dA + dD + r1 + r4 0 = K2 - C/B^2 0 = K3 - D/(A*B) 0 = r5 + 2*r4 - dE 0 = r5 -dI 0 = -r5 - r4 - dG 0 = -r1/2 - dH r1, r4 and r5 are variables expressed as functions of A, B, C, D, I, G and H, calculated previously in the integrated function. K2 and K3 are parameters. As I have two algebraic equations, I think this system is a DAE (Algebraic differential equation). I found in the package deSolve two functions that resolve DAE but I didn't manage to obtain a solution; it says that the variable dA cannot be found.
Show us your script where you define the function and run the DAE solver. Without that nobody can provide an answer. Berend.
Does anyone know how to solve this problem? Thank you Rapha?lle Carraud Rapha?lle Carraud [[alternative HTML version deleted]]
______________________________________________ R-help at 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.
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)
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.
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
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. -----Message d'origine----- De?: Berend Hasselman [mailto:bhh at xs4all.nl] Envoy??: jeudi 11 juillet 2013 12:33 ??: Rapha?lle Carraud Cc?: r-help at r-project.org Objet?: Re: [R] Differential problem
On 11-07-2013, at 12:05, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Sorry,
Here is the program I have until now:
reaction<-function(z, state, dval, parameters) {
with(as.list(c(state)),{
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
res6 <- r5 -dI
res7 <- -r5 - r4 - dG
res8 <- -r1/2 - dH
list(c(res1, res2, res3, res4, res5, res6, res7, res8))
}) # end with(as.list ...
}
Ti <- 273+90 #K
Pt <- 0.98*10^5 #Pa
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)
z <- seq(0, 1, by = 0.01)
library(deSolve)
out <- daspk(y = state, dy = dval, times = z, res = reaction, parms =
0)
head(out)
plot(out)
When I run this I get this error message:
Error in eval(expr, envir, enclos) : object 'Tr' not found
And shouldn't the first line of the reaction function be this
with(as.list(c(state,dval,parameters)),{
in stead of this
with(as.list(c(state)),{
The call of daspk also seems incorrect; shouldn't it be
out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = parameters)
And finally: where does variable P come from? Not defined anywhere!
Berend
-----Message d'origine----- De : Berend Hasselman [mailto:bhh at xs4all.nl] Envoy? : jeudi 11 juillet 2013 11:18 ? : Rapha?lle Carraud Cc : r-help at r-project.org Objet : Re: [R] Differential problem On 11-07-2013, at 09:13, Rapha?lle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:
Hello, I have the following differential equation system to solve, the variables I wish to obtain being A, B, C, D, E, I, G and H. 0 = -dA + dB + 2*dC - 2*r1 - 2*r5 0 = dA + dD + r1 + r4 0 = K2 - C/B^2 0 = K3 - D/(A*B) 0 = r5 + 2*r4 - dE 0 = r5 -dI 0 = -r5 - r4 - dG 0 = -r1/2 - dH r1, r4 and r5 are variables expressed as functions of A, B, C, D, I, G and H, calculated previously in the integrated function. K2 and K3 are parameters. As I have two algebraic equations, I think this system is a DAE (Algebraic differential equation). I found in the package deSolve two functions that resolve DAE but I didn't manage to obtain a solution; it says that the variable dA cannot be found.
Show us your script where you define the function and run the DAE solver. Without that nobody can provide an answer. Berend.
Does anyone know how to solve this problem? Thank you Rapha?lle Carraud Rapha?lle Carraud [[alternative HTML version deleted]]
______________________________________________ R-help at 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.
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
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