Skip to content
Prev 326678 / 398502 Next

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:

            
Show us your script where you define the function and run the DAE solver. Without that nobody can provide an answer.

Berend.