Ecuaciones Diferenciales Parciales (EDP) en R

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

Ecuaciones Diferenciales Parciales (EDP)

En contraste con las EDO, donde solo hay una variable independiente, las ecuaciones diferenciales parciales (PDE) contienen derivadas parciales con respecto a más de una variable independiente, por ejemplo t (tiempo) y x (una dimensión espacial). Para distinguir este tipo de ecuaciones de EDO, las derivadas se representan con el símbolo \(\partial\). \[ \frac{\partial y}{\partial t}=f(t,x,y,\frac{\partial y}{\partial x},p) \]

Las ecuaciones diferenciales parciales se pueden resolver subdividiendo una o más de las ecuaciones independientes continuas. variables en un número de celdas de cuadrícula, y reemplazando el derivadas por ecuaciones algebraicas discretas aproximadas, las llamadas diferencias finitas.

para casos en que varían con el tiempo, se acostumbra discretizar la(s) coordenada(s) espacial(es) solamente, mientras que el tiempo se deja en forma continua. Esto se llama el método de las líneas, y de esta manera, una PDE se traduce en una gran número de ecuaciones diferenciales ordinarias acopladas, que se pueden resolver con los habituales solucionadores de problemas de valor inicial. Esto aplica a PDE parabólicas, como la ecuación del calor, y PDE hiperbólicas, como la ecuación de onda.

para problemas invariantes en el tiempo, generalmente todas las variables independientes se discretizan y las derivadas se aproximan mediante ecuaciones algebraicas, que se resuelven mediante técnicas de búsqueda de raíces. Esta técnica se aplica a PDE elípticas.

R-package ReacTran proporciona funciones para generar diferencias finitas en una cuadrícula estructurada. Después,los casos variables en el tiempo resultantes se pueden resolver con funciones especialmente diseñadas del paquete deSolve, mientras que los casos invariantes en el tiempo se pueden resolver con métodos de resolución de raíz del paquete rootSolve.

Muchas ecuaciones diferenciales parciales se pueden resolver por aproximación numérica (diferenciación finita) después de reescribirlos como un conjunto de EDO.

Funciones tran.1D, tran.2D y tran.3D de Paquete R ReacTran implementar aproximaciones en diferencias finitas de la ecuación de transporte difusivon-advectivon que, para el caso 1-D, es:

\[ -\frac{1}{A_x}[\frac{\partial}{\partial x}A_x(-D\frac{\partial C}{\partial x})-\frac{\partial}{\partial x}(A_xCu)] \] Aquí D es el “coeficiente de difusión”, u es la “tasa de advección” y Ax es alguna propiedad (por ejemplo, el área de superficie) que depende de la variable independiente, x. Cabe señalar que la precisión de la finita las aproximaciones de diferencia no se pueden especificar en el funciones ReacTran. Depende del usuario asegurarse que las soluciones sean lo suficientemente precisas, p. incluyendo más puntos de cuadrícula.

Los modelos de reacción de difusión son una clase fundamental de modelos que describen cómo la concentración de materia, la energía, la información, etc. evoluciona en el espacio y el tiempo bajo la influencia del transporte difusivo y la transformación.

Ejemplo: EDP em una dimensión
Considere el modelo de reacción de difusión 1-D en [0,10]: \[ \frac{\partial C}{\partial t}=\frac{\partial}{\partial t}(D\frac{\partial C}{\partial x})-Q \] con C la concentración, t el tiempo, x la distancia del origen, Q, la tasa de consumo, y con condiciones de contorno (valores en los bordes del modelo): \[ \begin{eqnarray} \frac{\partial C}{\partial x}_{x=0}&=&0 \\ C_{x=10}&=&C_{ext} \end{eqnarray} \] Para resolver este modelo en R, primero el modelo 1-D Grid es definido; divide 10 cm (L) en 1000 cajas (N).

library(ReacTran)
Grid <- setup.grid.1D(N = 1000, L = 10)

La ecuación del modelo incluye un término de transporte, aproximado por la función ReacTran tran.1D y un término de consumo (Q). La condición límite aguas abajo, prescrita como una concentración (C.down) debe especificarse, el gradiente cero en el límite de corriente arriba es el valor predeterminado:

pde1D <-function(t, C, parms) {
tran <- tran.1D(C = C, D = D,
C.down = Cext, dx = Grid)$dC
list(tran - Q) # return value: rate of change
}

Los parametros del modelo son:

D <- 1 # diffusion constant
Q <- 1 # uptake rate
Cext <- 20

En una primera aplicación, el modelo se resuelve para estado estacionario, que recupera la condición en la que el las concentraciones son invariantes:

\[ 0=\frac{\partial}{\partial x}(D \frac{\partial C}{\partial x})-Q \] En R, las condiciones de estado estacionario se pueden estimar usando funciones del paquete rootSolve que implementan entre otros, un algoritmo de Newton-Raphson (Presione et al., 2007). Para modelos unidimensionales, constante.1D es más eficiente. La “conjetura” inicial del estado estacionario la solución (y) no es importante; aquí tomamos simplemente N números al azar. El argumento nspec = 1 informa al solver que solo se describe un componente.

Aunque un sistema de 1000 ecuaciones necesita ser resuelto, esto toma solo una fracción de segundo:

library(rootSolve)
print(system.time(
std <- steady.1D(y = runif(Grid$N),
func = pde1D, parms = NULL, nspec = 1)
))

   user  system elapsed 
   0.04    0.02    0.06 

Los valores de las variables de estado (y) se trazan contra la distancia, en medio de las celdas de la cuadrícula (Cuadrícula$x.mid).

plot (Grid$x.mid, std$y,
      type = "l",lwd = 2,
      main = "steady-state PDE",
      xlab = "x",
      ylab = "C",
      col = "blue")

Solución de estado estacionario del modelo de reacción de difusión 1-D.

La solución analítica se compara bien con la aproximación numérica:

analytical <- Q/2/D*(Grid$x.mid^2 - 10^2) + Cext
max(abs(analytical - std$y))

[1] 1.250004e-05

A continuación, el modelo se ejecuta dinámicamente durante 100 veces unidades usando la función deSolve ode.1D, y comenzando con una concentración uniforme:

require(deSolve)
times <- seq(0, 100, by = 1)
system.time(
out <- ode.1D(y = rep(1, Grid$N),
times = times, func = pde1D,
parms = NULL, nspec = 1)
)

   user  system elapsed 
   0.34    0.26    0.61 

Aquí, out es una matriz, cuya primera columna contiene los tiempos de salida, y las siguientes columnas los valores de las variables de estado en los diferentes recuadros; imprimimos el primeras columnas de las últimas tres filas de esta matriz:

tail(out[, 1:4], n = 3)

       time         1         2         3
 [99,]   98 -27.55783 -27.55773 -27.55754
[100,]   99 -27.61735 -27.61725 -27.61706
[101,]  100 -27.67542 -27.67532 -27.67513

Trazamos el resultado usando un color azul-amarillo-rojo esquema, y usando la imagen del método S3 de deSolve. La figura 6 muestra que, a medida que pasa el tiempo, se desarrollan gradientes desde la distribución uniforme, hasta que el sistema casi alcanza el estado estacionario al final de la simulación.

image(out, xlab = "tiempo, dias",ylab = "Distancia, cm",main = "EDP", add.contour = TRUE)

Solución dinámica del modelo de reacción de difusión 1-D.

Cabe señalar que el modelo de estado estacionario es efectivamente un problema de valores en la frontera, mientras que el modelo transitorio es un prototipo de una ecuación diferencial parcial “parabólica” . Mientras que R también puede resolver los otros dos principales clases de PDE, es decir, del tipo “hiperbólica” y “elíptica”, está mucho más allá del alcance de este documento para elaborar sobre eso.