[R] Differential problem

2013-07-26 Thread Raphaëlle Carraud
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   12 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]]

__
R-help@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.


[R] Ode error message

2013-07-16 Thread Raphaëlle Carraud
Hello,

I am creating a program with R to solve a differential equation system. 
However, I get the following message I do not understand :

> out <- ode(y = state, times = z, func = liquide, parms = 0, atol = 0)
DLSODA-  EWT(I1) is R1 .le. 0.0
In above message, I1 = 1
 
In above message, R1 = 0
 
Error in lsoda(y, times, func, parms, ...) : 
  illegal input detected before taking any integration steps - see written 
message

or this one when I tried modifying the atoll value :

> out <- ode(y = state, times = z, func = liquide, parms = 0, atol = 10^-14)
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
  such that in the machine, T + H = T on the next step
 (H = step size). Solver will continue anyway.
In above message, R1 = 0, R2 = 0
 
DINTDY-  T (=R1) illegal
In above message, R1 = 1
 
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R1 = 0, R2 = 0
 
DINTDY-  T (=R1) illegal
In above message, R1 = 2
 
  T not in interval TCUR - HU (= R1) to TCUR (=R2)
In above message, R1 = 0, R2 = 0
 
DLSODA-  Trouble in DINTDY.  ITASK = I1, TOUT = R1
In above message, I1 = 1
 
In above message, R1 = 2

Here is my program. I also tried changing the initial values but it does not 
work well.

liquide <- function(z, state, parameters) {
  with(as.list(c(state,parameters)),{
# rate of change
  
Tr <- 273+90
C <- CA + CB + CC + CD + CE + CI + CG + CJ + CK + CH

K32 <- 6.54*10^4  
K33 <- 1.37*10^4  
K34 <- 330 
K35 <- 5.81*10^4  
kf2 <- 1.37*10^3  
kf3 <- 1.37*10^3 
kf4 <- 8.68*10^5  
kf5 <- 157.2

K2 <- 10^1.37
K3 <- 10^(-3.35) 

r1 <- kf4*CD - kf4/K34*CE^2
r2 <- kf3*CA*CB - kf3/K33*CD
r3 <- kf2*CA^2 - kf2/K32*CC
r4 <- kf5*CC - kf5/K35*CE*CI^2 


dCA <- -r2  # dNO/dt
dCB <- -r3 - r2   # dNO2/dt
dCC <- r3/2 - r4 # dN2O4/dt
dCD <- r2 - r1# dN2O3/dt
dCE <- 2*r1 + r4# dHNO2/dt
dCI <- r4 # dHNO3/dt
dCG <- -r4 - r1   # dH2O/dz
dCH <- (dCE + dCI)/((K2 + K3)*(CE + CI))  # dH/dz
dCJ <-  (CH*dCI - CI*dCH)/(K3*CH^2)  # dNO3-/dz
dCK <-  (CH*dCE - CE*dCH)/(K2*CH^2)# dNO2-/dz

   

list(c(dCA, dCB, dCC, dCD, dCE, dCI, dCG, dCH, dCJ, dCK))
  })   # end with(as.list ...
}


Ti <- 273+90   # K
Ct <- 5100   # mol/m^3

state <-c(CA = 0,   # mol/m^3 NO2
  CB = 0,   # mol/m^3 NO
  CC = 0,# mol/m^3 N2O4
  CD = 0,  # mol/m^3 N2O3
  CE = 50,  # mol/m^3 HNO2
  CI = 50, # mol/m^3 HNO3
  CG = 5000, # mol/m^3 H2O
  CH = 0, # mol/m^3 H+
  CJ = 0,# mol/m^3 NO3-
  CK = 0) # mol/m^3 NO2-
  
parameters <- c(Ct = 5100)

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

library(deSolve)
out <- ode(y = state, times = z, func = liquide, parms = 0, atol = 10^-14)
head(out)
plot(out)

Thank you

__
R-help@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.


Re: [R] Differential problem

2013-07-11 Thread Raphaëlle Carraud
Ok, now it's good (Pt was in my workplace so it worked for me, I am not used to 
R using these value to make the program run so I hadn't looked...)

reaction<-function(z, state, dval, parameters) {
  with(as.list(c(state,dval,parameters)),{
# 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

Pt <- 0.98*10^5
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)
head(out)
plot(out)


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


On 11-07-2013, at 13:53, Raphaëlle Carraud  
wrote:

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

Doesn't run.
Since variable Pt is not defined when you calculate vector state. So define Pt 
<- .. before xi as in the original example.

In the function reaction isn't variable P just Pt from the parameter vector?
If so then either do P <- Pt or just use Pt directly (but see next remark).


> z <- seq(0, 1, by = 0.01)  # en seconde
> 
> librar

Re: [R] Differential problem

2013-07-11 Thread Raphaëlle Carraud
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)
head(out)
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)
Error in eval(expr, envir, enclos) : object 'dA' not found.

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:b...@xs4all.nl] 
Envoyé : jeudi 11 juillet 2013 12:33
À : Raphaëlle Carraud
Cc : r-help@r-project.org
Objet : Re: [R] Differential problem


On 11-07-2013, at 12:05, Raphaëlle Carraud  
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 <

Re: [R] Differential problem

2013-07-11 Thread Raphaëlle Carraud
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:b...@xs4all.nl] 
Envoyé : jeudi 11 juillet 2013 11:18
À : Raphaëlle Carraud
Cc : r-help@r-project.org
Objet : Re: [R] Differential problem


On 11-07-2013, at 09:13, Raphaëlle Carraud  
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@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.

__
R-help@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.


[R] Differential problem

2013-07-11 Thread Raphaëlle Carraud
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]]

__
R-help@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.


[R] Recherche de fonction

2013-07-10 Thread Raphaëlle Carraud
Bonjour,

Je souhaite  résoudre le couple d'équation différentielles suivant :

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

en ayant connaissance des valeurs initiales de dA, dB, dC, dE, dI, dG, dH, r1, 
r2, r4, r5, K2, K3, A, B, C et D.

J'ai essayé plusieurs fonctions mais comme je ne peux pas lui faire calculer 
une des dérivée de laquelle découlerait les autre, il n'arrive pas à me fournir 
la solution.
Je n'ai pas vu d'exemple qui pourrai s'assimiler à celui-ci dans la 
documentation.

Est-il possible de résoudre ce problème sur R ?

Merci

Cordialement,
Raphaëlle Carraud

[[alternative HTML version deleted]]

__
R-help@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.