# [R] Differential problem

Raphaëlle Carraud raphaelle.carraud at oc-metalchem.com
Thu Jul 11 13:53:29 CEST 2013

```Sorry for the bug, I had eliminated some lines to avoid making the program too big. Here is the version that works :

reaction<-function(z, state, dval, parameters) {
with(as.list(c(state)),{
# rate of change

Tr <- 273+90
P <- 0.98*10^5

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             #dHNO2/dz
res6 <- r5 -dI                   #dHNO3/dz
res7 <- -r5 - r4 - dG             #dH2O/dz
res8 <- -r1/2 - dH                #dO2/dz

list(c(res1, res2, res3, res4, res5, res6, res7, res8))
})   # end with(as.list ...
}

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)  # en seconde

library(deSolve)
#out <- ode(y = state, times = z, func = reaction, parameters)

out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0)
plot(out)

I obtain the following message:

> library(deSolve)
> #out <- ode(y = state, times = z, func = reaction, parameters)
>
> out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = 0)

I tried adding the dval and parameters as you said:
with(as.list(c(state,dval,parameters)),{

I get the following message:

Warning messages:
1: In daspk(y = state, dy = dval, times = z, res = reaction, parms = 0) :
matrix of partial derivatives is singular with direct method-some equations redundant
2: In daspk(y = state, dy = dval, times = z, res = reaction, parms = 0) :
Returning early. Results are accurate, as far as they go

For the calling of the daspk function, I followed the documentation, where you have the same inversion:

daefun <- function(t, y, dy, parameters) {
+ res1 <- dy[1] + y[1] - y[2]
+ res2 <- y[2] * y[1] - t
+
+ list(c(res1, res2))
+ }
> library(deSolve)
> yini <- c(1, 0)
> dyini <- c(1, 0)
> times <- seq(0, 10, 0.1)
> ## solver
> system.time(out <- daspk(y = yini, dy = dyini, times = times, res = daefun, parms = 0))

Is it wrong? When I modify the order, I obtain again that object dA is not found, so I guessed the doc was right.

-----Message d'origine-----
De : Berend Hasselman [mailto:bhh at xs4all.nl]
Envoyé : jeudi 11 juillet 2013 12:33
À : Raphaëlle Carraud
Cc : r-help at r-project.org
Objet : Re: [R] Differential problem

On 11-07-2013, at 12:05, Raphaëlle Carraud <raphaelle.carraud at oc-metalchem.com> wrote:

> 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)
> plot(out)
>

When I run this I get this error message:

And shouldn't the first line of the reaction function be this

with(as.list(c(state,dval,parameters)),{

with(as.list(c(state)),{

The call of daspk also seems incorrect; shouldn't it be

out <- daspk(y = state, dy = dval, times = z, res = reaction, parms = parameters)

And finally: where does variable P come from? Not defined anywhere!

Berend

> -----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]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help