Objetivo

En este documento modelamos cómo la cafeína disminuye en el organismo usando ecuaciones diferenciales ordinarias (EDO) resueltas con el paquete deSolve.

Modelo 1: Eliminación simple \(dX/dt = -kX\)

Definición de la EDO y función para deSolve

Preparación del entorno

Explicación. Aquí se instala y se carga el paquete deSolve, que sirve para resolver ecuaciones diferenciales numéricamente en R.

library(deSolve)

Definición de la EDO y función para deSolve

Explicación.

Convierte el estado en lista para acceder a X y al mismo tiempo crea un vector dxdt para almacenar la derivada.

Donde X es la cantidad de cafeína y k la constante de eliminación.

Devuelve la derivada como lista.

ed.sol = function(t, state, parms)
{
  with(as.list(state),
       {
         dxdt = rep(0, length(state))
         dxdt = -k*X
         return(list(dxdt))
       })
}

Parámetros, condición inicial y malla temporal

Explicación. X = 70 cantidad inicial de cafeína (mg).

k = 0.2 constante de eliminación (1/hora).

t = seq(0,30,1) el tiempo de simulación (0 a 30 horas en pasos de 1).

init = c(X=X) condición inicial empaquetada para ode.

X = 70
k = 0.2
t = seq(0,30,1)
init = c(X=X)

Integración numérica con ode()

Explicación. Aquí se llama a la función ode() para integrar la ecuación:

-y = init estado inicial.

-times = t los tiempos en que queremos la solución.

-func = ed.sol la función que define la EDO.

-parms = NULL no hay parámetros adicionales.

El resultado ) es una tabla con tiempo y valores de X. head() muestra las primeras filas.

miOutput = ode(y=init, times = t, func = ed.sol, parms = NULL)
head(miOutput)
##      time        X
## [1,]    0 70.00000
## [2,]    1 57.31122
## [3,]    2 46.92244
## [4,]    3 38.41684
## [5,]    4 31.45300
## [6,]    5 25.75162

Gráfica del modelo 1

Explicación. Graficamos \(X(t)\) vs tiempo para visualizar la decadencia exponencial. Etiquetamos ejes con unidades.

plot(miOutput, type = "l", xlab = "Tiempo(hr)", ylab = "Cafeína(mg)")