Rafael Martínez Fonseca

R-package

library(deSolve)   #PVI,ESD,EDP y EDO
library(ReacTran) #EDP
library(bvpSolve) #PVI
library(rootSolve) #EDO
library(PBSddesolve) #DDE

Ecuación Diferencial Algebraica(EDA)

Una ecuación diferencial algebraica esta conformada para una ecuación diferencial acompañada de una ecuación algebraica. \[ \begin{eqnarray} y'&=& f(t,y,p) \\ 0&=& g(t,y,p) \end{eqnarray} \] El índice de un EDA es el número de diferenciaciones necesarias hasta obtener un sistema compuesto únicamente por EDO. La función daspk del paquete deSolve resuelve (relativamente simple) EDA de índice como máximo 1, mientras que la función radau resuelve EDA de índice hasta 3.

Ejemplo:

El llamado “problema de Rober” describe una reacción autocatalítica entre tres especies químicas, \(y1\), \(y2\) e \(y3\). El problema puede ser formulado ya sea como una EDO o como la siguiente EDA: \[ \begin{eqnarray} y_1'&=& -0.04y_1+10^4y_2y_3 \\ y_2'&=& 0.04y_1-10^4y_2y_3 -310^7y_2^2\\ 1&=& y_1+y_2+y_3 \end{eqnarray} \]

donde las dos primeras ecuaciones son diferenciales ecuaciones que especifican la dinámica de la química especies y1 e y2, mientras que la tercera ecuación algebraica asegura que la concentración suma de los tres la especie permanece 1.

El DAE tiene que ser especificado por la función residual en lugar de las tasas de cambio.(Como en EDO).

\[ \begin{eqnarray} r_1'&=& -y_1'-0.04y_1+10^4y_2y_3 \\ r_2'&=& -y_2'+0.04y_1-10^4y_2y_3 -310^7y_2^2\\ r_3&=& -1+y_1+y_2+y_3 \end{eqnarray} \] Implementación en R:

edafun<-function(t, y, dy, parms) {
res1 <- - dy[1]-0.04*y[1]+1e4*y[2]*y[3]
res2 <- - dy[2]+0.04*y[1]-1e4*y[2]*y[3]-3e7*y[2]^2
res3 <- y[1]+y[2]+y[3]-1
list(c(res1, res2, res3),
error = as.vector(y[1]+y[2]+y[3])- 1)
}
yini <- c(y1 = 1, y2 = 0, y3 = 0)
dyini <- c(-0.04, 0.04, 0)
times <- 10 ^ seq(-6,6,0.1)

Los argumentos de entrada de la función daefun son los tiempo actual(t), los valores de las variables de estado y sus derivadas (y, dy) y los parámetros (parms).Devuelve los residuales, concatenados y una salida. variable, el error en la ecuación algebraica. Este último se añade para comprobar la exactitud de los resultados. Para DAE resueltos con daspk, tanto las variables de estado como sus derivados deben inicializarse. Aquí nos aseguramos de que las condiciones iniciales para y obedezcan la restricción algebraica, mientras que también la condición inicial de las derivadas es consistente con la dinámica.

library(deSolve)
print(system.time(out <-daspk(y = yini,dy = dyini,times = times,res = edafun,
parms = NULL)))
   user  system elapsed 
   0.08    0.01    0.10

Se puede utilizar un método de trazado S3 para trazar todas las variables a la vez:

plot(out, ylab = "conc.", xlab = "time",
type = "l", lwd = 2, log = "x")
mtext("IVP DAE", side = 3, outer = TRUE,line = -1)

Hay un cambio inicial muy rápido en las concentraciones, principalmente debido a la rápida reacción entre y1 y y2 y entre y2. Después de eso, la reacción lenta de y1 con y2 hace que el sistema cambie mucho más suavemente. Esto es típico de los problemas rígidos.

Solución del problema DAE para las sustancias y1,y2,y3; error de balance de masa: desviación de la suma total de uno.