R help - Jul 2013 - Differential problem

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.

Does anyone know how to solve this problem?

Thank you

Raphaëlle Carraud


Raphaëlle Carraud

	[[alternative HTML version deleted]]

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.

Hello,

I am trying to solve the following algebraic differential equation :

dy1 = h - dy3
dy2 = g - dy4
y3 = K1*y1/(y3+y4)
y4 = K2*y2/(y3+y4)

I tried using the function daspk, but it gives me the following error :
> out <- daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0)
daspk--  warning.. At T(=R1) and stepsize H (=R2) the
      nonlinear solver failed to converge
      repeatedly of with abs (H) = HMIN &g, 0
Warning messages:
1: In daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0) :
  repeated convergence test failures on a step - inaccurate Jacobian or
preconditioner?
2: In daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0) :
  Returning early. Results are accurate, as far as they
go
> head(out)
     time   1    2 3 4  5  6     7     8     9
[1,]    0 0.9 1.33 0 0 10 20 1e-04 0.001 0.001
[2,]    0 0.9 1.33 0 0 10 20 1e-04 0.001 0.001

I read the documentation but it only says that daspk can only solve DAE of index
one maximum, which I do not understand how to determine. I was wondering if I
had to use the function radau but I do not get how to do the M matrix? I did not
managed to get it with the example.
Could it be an error in my program? Here it is :

liquide <- function(z, C, dC, parameters) {
  with(as.list(c(C,dC,parameters)),{

    Ct = C[1] + C[2] + C[3] + C[4] + C[5] + C[6] + C[7] + C[8] + C[9] + Ceau
    V <- 1

    K32 <- 6.54*10^2              # m^3/mol
    K33 <- 1.37*10^-2             # m^3/mol
    K34 <- 0.330                  # sans unité
    K35 <- 5.81*10^1              # sans unité

    kf2 <- 1.37                   # m^3/mol
    kf3 <- 1.37                   # m^3/mol
    kf4 <- 8.68*10^-1             # m^3
    kf5 <- 157.2*10^-6            # m^3

    K2 <- 10^1.44*10^3            # mol/m^3
    K3 <- 10^(-3.224)*10^3        # mol/m^3
    Ke <- 10^-11                  # mol/m^3

    r1 <- kf4*C[4]/V - kf4/K34*C[5]^2/(Ct*V)
    r2 <- kf3*C[1]*C[2] - kf3/K33*C[4]
    r3 <- kf2*C[1]^2 - kf2/K32*C[3]
    r4 <- kf5*C[3]/V - kf5/K35*C[5]*C[6]*C[8]/(Ct^2*V)

    res1 <- dC[1] + r2 + r3                     # dNO2/dt
    res2 <- dC[2] + r2                          # dNO/dt
    res3 <- dC[3] - r3 + r4                     # dN2O4/dt
    res4 <- dC[4] - r2 + r1                     # dN2O3/dt
    res5 <- dC[5] -2*r1 - r4 + dC[7]            # dHNO2/dt
    res6 <- dC[6] - r4 + dC[8]                  # dHNO3/dt
    res7 <- C[7] - K3*C[5]/C[9]                 # dNO2-/dt
    res8 <- C[8] - K2*C[6]/C[9]                 # dNO3-/dt
    res9 <- dC[9] - dC[8] - dC[7]               # dCH/dt

    list(c(res1, res2, res3, res4, res5, res6, res7, res8, res9))
  })
}

yini <- c(0.9, 1.33, 0, 0, 10,  20, 0.001, 22.9, 23)
dyini <- c(1,1,1,1,1,1,1,1,1)

Qm <- 4                       # kg/h
Ceau <- (Qm - yini[1]*0.046 - yini[2]*0.03 + yini[3]*0.092 + yini[4]*0.076 +
           yini[5]*0.062 + yini[6]*0.063)/0.018    # mol/m^3

parameters <- c(Qm = Qm,
                Ceau = Ceau)

liquide(20,yini, dyini,parameters)

z <- seq(0, 120, by = 1)  # en seconde

library(deSolve)
out <- daspk(y = yini, dy = dyini, times = z, res = liquide, parms = 0)
head(out)
plot(out)

Thanks

Raphaëlle Carraud

	[[alternative HTML version deleted]]