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

