Differential problem
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.