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:b...@xs4all.nl]
Envoyé : jeudi 11 juillet 2013 11:18
À : Raphaëlle Carraud
Cc : r-help@r-project.org
Objet : Re: [R] Differential problem
On 11-07-2013, at 09:13, Raphaëlle Carraud <raphaelle.carr...@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@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.
______________________________________________
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.
