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]]