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

Mostramos por ejemplo las funciones del paquete desolve.

ode(y, times, func, parms, method = c("lsoda",
"lsode", "lsodes", "lsodar",
"vode", "daspk", "euler", "rk4",
"ode23", "ode45", "radau", "bdf",
"bdf_d", "adams", "impAdams",
"impAdams_d"), ...)

Ecuaciones Diferenciales Ordinarias (EDO)

Las ecuaciones diferenciales ordinarias describen el cambio de una variable de estado \(y\) como una función \(f\) de una independiente variable \(t\) (por ejemplo, tiempo o espacio), de sí mismo y, opcionalmente, un conjunto de otras variables p, a menudo llamadas parámetros:

\[ y^\prime=\frac{dy}{dx}=f(t,y,p) \] en ocasiones para dar solución a una EDO se necesitan valores extras, por lo general son condiciones iniciales y/o los llamados valores en la frontera que se describen mas adelante.

Problemas de Valor Inicial (PVI)

Si en la ecuación diferencial hay una condición extra se les llama problemas de valor inicial por las condiciones iniciales corresponden a información de la función en un punto.(ejemplo: \(y(x=0)=0\)).

Existen dos algoritmos numéricos para estos problemas, el llamado Runge-Kutta y el multipaso lineal. Finalmemte contamos dos importantes familias, La Adams familia y las fórmula de diferenciación hacia atrás(BDF:backward diffentiation formulae).

Otra distinción importante es entre explícito y métodos implícitos, donde los últimos métodos pueden resolver una clase particular de ecuaciones (llamadas “rígidas” ecuaciones) donde los métodos explícitos tienen problemas con estabilidad y eficiencia. La rigidez ocurre, por ejemplo, si un problema tiene componentes con diferentes tasas de variación según la variable independiente. Muy a menudo habrá una compensación entre el uso de métodos explícitos que requieren poco trabajo por paso de integración y métodos implícitos que son capaces de dar pasos de integración más grandes, pero necesita (mucho) más trabajo para un paso.

En R, los problema de valor inicial se resuelven con la función del paquete desolve, cual implementación nos permite resolver muchas EDO del código vode,las ecuaciones diferenciales algebraicas se resuelven con daspk,todas estas funciones pertenecen a los métodos multipasos y comprenden las metodos Adams de diferenciación hacia atrás. Los primeros métodos son explícitos, los últimos implícitos. Además, este paquete contiene una implementación de-novo de un solucionador de Runge-Kutta bastante general. Finalmente, el método implícito de Runge Kutta radau ha sido añadido recientemente.

Ejemplo:Una EDO de valor inicial

Considere la famosa ecuación de van der Pol, que describe un oscilador no conservativo con amortiguamiento no lineal y que fue desarrollado originalmente para circuitos eléctricos que emplean tubos de vacío. La oscilación se describe mediante una EDO de segundo orden:

\[ z''-\mu(1-z^2)z'+z=0 \] Tal sistema se puede reescribir rutinariamente como un sistema de dos EDO de 1er orden, si sustituimos \(y_1=z\). \[ y_1=z \xrightarrow{Derivando} y_1'=z'=y_2 \xrightarrow{Derivando} z''=y_2' \] Escribimos nuestro sistema de EDO.

\[ \begin{eqnarray} y_1'&=& y_2 \\ y_2'&=& \mu(1-y_1^2)y_2-y_1 \end{eqnarray} \] Hay un parámetro, \(\mu\), y dos diferenciales variables, \(y_1\) e \(y_2\) con valores iniciales (en \(t = 0\)): \[ \begin{eqnarray} y1_{(t=0)}&=& 2 \\ y2_{(t=0)}&=& 0 \end{eqnarray} \] La ecuación de van der Pol se usa a menudo como prueba problema para los solucionadores de ODE, ya que, para \(\mu\) grandes, su dinámica consiste en partes donde la solución cambia muy lentamente, alternando con regiones de muy agudo cambios. Esta “rigidez” hace que la ecuación sea bastante desafiante de resolver.

En R, este modelo se implementa como una función de van der Pol(vdpol) cuyas entradas son el tiempo actual (t), los valores de las variables de estado (y) y los parámetros (mu); la función devuelve una lista con como primer elemento el derivados, concatenados.

vdpol <- function(t,y,mu){list(c(y[2],mu*(1-y[1]^2)*y[2]-y[1]))} #Escritura del sistema

Después de definir la condición inicial del estado variables (c_inicial), el modelo se resuelve y la salida escrito en puntos de tiempo (tiempos) seleccionados, usando la función de integración ode de Solve. La rutina predeterminada lsoda, que ode invoca automáticamente cambia entre métodos rígidos y no rígidos, dependiendo del problema. Ejecutamos el modelo para una rigidez típica (mu = 1000) y no rígido (mu = 1) situación:

library(deSolve)
c_inicial <- c(y1 = 2, y2 = 0)   #Valores iniciales 
rigido <- ode(y = c_inicial, func = vdpol,times = 0:3000, parms = 1000)
no_rigido <- ode(y = c_inicial, func = vdpol,times = seq(0, 30, by = 0.01),parms = 1)

El modelo devuelve una matriz, de clase deSolve, con en su primera columna los valores de tiempo, seguido de los valores de las variables de estado:

head(rigido, n = 3)
     time       y1            y2
[1,]    0 2.000000  0.0000000000
[2,]    1 1.999333 -0.0006670373
[3,]    2 1.998666 -0.0006674088
head(no_rigido, n = 3)
     time       y1          y2
[1,] 0.00 2.000000  0.00000000
[2,] 0.01 1.999901 -0.01970323
[3,] 0.02 1.999608 -0.03882242

Las cifras se generan utilizando el método de trazado S3 para objetos de la clase deSolve:

plot(rigido, type = "l", which = "y1",lwd = 2,col='blue', ylab = "y",main = "PVI EDO, Rigido")

plot(no_rigido, type = "l", which = "y1",lwd = 2,col='blue', ylab = "y",main = "PVI EDO,No Rigido")

Solución de la ecuación de van der Pol, una ecuación diferencial ordinaria de valor inicial, no rígida caso, \(\mu = 1\).

Solver No Rigido Rigido
deSolve 0.37 271.19
Isoda 0.26 0.23
adams 0.13 616.13
bdf 0.15 0.22
radau 0.53 0.72

Comparación de solucionadores para un rígido y un parametrización no rígida de la ecuación de van der Pol (tiempo en segundos, valores medios de diez simulaciones en una CPU AMD AM2 X2 3000).

Una comparación de tiempos para dos solucionadores explícitos, el método de Runge-Kutta (oda23) y el de adams método, con el solucionador multipaso implícito (bdf,fórmula de diferenciación hacia atrás) muestra una clara ventaja para este último en el caso rígido (Figura 1). los El solucionador predeterminado (lsoda) no es necesariamente el más rápido, pero muestra un comportamiento robusto debido a la detección automática de rigidez. Utiliza el Adams multipaso explícito método para el caso no rígido y el método BDF para el caso rígido. La precisión es comparable para todos solucionadores con atol =``rtol```= 10^{−6} , el valor por defecto.

Problemas de Valores en la Frontera (PVF)

Si las condiciones adicionales se especifican en diferentes valores de la variable independiente, las ecuaciones diferenciales se denominan problemas de valores en la frontera (PVF).

El paquete bvpSolve implementa tres métodos para resolver problemas de valores en la frontera. El método de solución más simple es el único método de disparo, que combina la integración del problema de valor inicial con un algoritmo de búsqueda de raíces no lineal. Dos métodos de solución más estables implementan un código Runge Kutta (MIRK) monoimplícito, basado en el código FORTRAN twpbvpC y la colocación método, basado en el código FORTRAN colnew. Algunos problemas de valores en la frontera también se puede resolver con funciones de los paquetes Re acTran y rootSolve.

Ejemplo: \[ \xi y''-y=-(\xi \pi^2+1)\cos(\pi x) \quad x \in[-1,1] \] y las condicionales: \[ \begin{eqnarray} y{(x=-1)}&=& 0 \\ y{(x=+1)}&=& 0 \end{eqnarray} \] La ecuación de segundo orden primero se reescribe como dos ecuaciones de primer orden: \[ \begin{eqnarray} y_1'&=& y_2 \\ y_2'&=& \frac{1}{\xi}[y_1-(\xi \pi^2+1)\cos(\pi x)] \end{eqnarray} \] Implementación en R.

Problema <- function(x, y, xi) {
list(c(y[2],1/xi * (y[1] - (xi*pi*pi+1) * cos(pi*x))))
  }

Con valores decrecientes de \(\xi\), este problema se convierte en cada vez más difícil de resolver. Usamos tres valores de \(\xi\), y resolvemos el problema con el disparo,MIRK y el método de colocación. Observe cómo las condiciones iniciales yini y las condiciones al final del intervalo de integración yend se especifican, donde NA indica que el valor no es conocido. La variable independiente se llama aquí x.

library(bvpSolve)

x <- seq(-1, 1, by = 0.01)
shoot <- bvpshoot(yini = c(0, NA),yend = c(0, NA), x = x, parms = 0.01,func = Problema)

twp <- bvptwp(yini = c(0, NA), yend = c(0,NA), x = x, parms = 0.0025,func = Problema)

coll <- bvpcol(yini = c(0, NA),yend = c(0, NA), x = x, parms = 1e-04,func = Problema)

La aproximación numérica generada por bvptwp está muy cerca de la solución analítica,para \(\xi= 0.0025\).

xi <- 0.0025
analytic <- cos(pi*x)+exp((x-1)/sqrt(xi))+exp(-(x+1)/sqrt(xi))
max(abs(analytic-twp[,2]))
[1] 7.788209e-10

Se observa una discrepancia baja similar (\(4·10^{−11}\)) para el \(\xi= 0.0001\) resuelto por bvpcol; el tiroteo método es considerablemente menos preciso (\(1.4·10^{−5}\)), aunque la misma tolerancia (atol = \(·10^{−8}\)) se utilizó para todas las carreras. La gráfica muestra cómo la forma de la solución se ve afectada por el parámetro \(\xi\), haciéndose cada vez más más empinada cerca de los límites, y por lo tanto más y más difícil de resolver, a medida que \(\xi\) se hace más pequeño.

plot(shoot[, 1], shoot[, 2], type = "l", lwd = 2,
ylim = c(-1, 1), col = "blue",
xlab = "x", ylab = "y", main = "BVP ODE")
#
lines(twp[, 1], twp[, 2], col = "red", lwd = 2)
#
lines(coll[, 1], coll[, 2], col = "green", lwd = 2)
#
legend("topright", 
legend = c("0.01", "0.0025","0.0001"), col = c("blue", "red", "green"),
title = expression(xi), lwd = 2)

Solución del problema BVP ODE, para diferentes valores de parámetro \(\xi\).