Diagrama de fases en R

Un diagrama de fases permite visualizar las soluciones de un sistema dinámico.

En particular puede ser útil para ver si existen soluciones periódicas o para visualizar el tipo de estabilidad de los puntos de equilibrio del sistema.

Diagrama de fases de un modelo de competencia de 2 especies junto con las curvas isóclinas del sistema y 4 trayectorias solución

En R se puede utilizar el paquete "pplane" para realizar un diagrama de fases (ver al final de este post sobre cómo instalar este paquete).

Librerías necesarias

Una vez instalado el paquete "pplane", se lo debe cargar mediante el comando: library(pplane).

También se debe cargar el paquete "deSolve" que permite resolver ecuaciones diferenciales: library(deSolve). Este último paquete "deSolve" está disponible para instalar directamente desde la línea de comandos de R.

Sistema dinámico

El sistema dinámico que voy a tomar como ejemplo es un modelo de competencia de dos especies $x_{1}(t)$ y $x_{2}(t)$ por los recursos de cierto ambiente. El mismo se puede describir mediante el siguiente sistema de ecuaciones diferenciales:

$\frac{dx_1}{dt} = r_1 x_1 ( 1 - \frac{x_1}{k_1} - \alpha_{12} \frac{x_2}{k_1} )$

$\frac{dx_2}{dt} = r_2 x_2 ( 1 - \frac{x_2}{k_2} - \alpha_{21} \frac{x_1}{k_2} )$

Donde $r_1$, $r_2$, $k_1$, $k_2$, $\alpha_{12}$ y $\alpha_{21}$ son parámetros del sistema.

Parámetros

Si en el instante inicial las poblaciones toman valores $x_{10}$ y $x_{20}$ y se deja evolucionar el sistema, la situación final será una de las siguientes:

• La especie 1 se extingue y la especie 2 alcanza una población igual a $k_2$
• La especie 2 se extingue y la especie 1 alcanza una población igual a $k_1$
• Ambas especies coexisten alcanzando poblaciones respectivas $x_1=\frac{\alpha_{12}k_2-k_1}{\alpha_{12}\alpha_{21}-1}$ y $x_2=\frac{\alpha_{21}k_1-k_2}{\alpha_{12}\alpha_{21}-1}$

Cuál de estas situaciones ocurrirá dependerá del valor utilizado para los 6 parámetros del modelo (en el caso en que ambos puntos de extinción son estables también dependerá de las condiciones iniciales).

En este ejemplo en particular utilizo los siguientes valores:

• r1=r2=1
• k1=400 y k2=300
• $\alpha_{12}=0.9$ y $\alpha_{21} = 0.4$

Para estos valores se verifica:

1. $\alpha_{12} < \frac{k_1}{k_2}=1.333$
2. $\alpha_{21} < \frac{k_2}{k_1}=0.75$

Para este modelo, esto implica que el punto de coexistencia es un punto de equilibrio estable y los restantes son puntos de equilibrio inestables. Es decir que el resultado final será que las especies van a coexistir.

Código utilizado

La definición del modelo en R y los comandos utilizados para dibujar el plano de fases, algunas trayectorias y las curvas isóclinas se encuentran en el siguiente código que explico en las siguientes secciones:

	########
	# mayo 2013
	# plano de fases en R - modelo de competencia logístico
	# CharlieXRay en blogspot.com
	########
	########
	# librerias necesarias
	########
	library(pplane) # para dibujar el plano de fases
	library(deSolve) # para resolver el sistema de ecuaciones diferenciales
	########
	# defino el sistema dinámico a utilizar
	########
	modelo <- function(parametros){
	function(x1,x2){
	with(as.list(parametros),{ # con los componentes r1, k1, etc que conforman el array parametros
	dx1 <- x1*( 1 - (x1/k1) - a12*(x2/k1) )
	dx2 <- x2*( 1 - (x2/k2) - a21*(x1/k2) )
	return( c(dx1,dx2) ) # retorna las derivadas
	})
	}
	}
	########
	# parametros a utilizar
	########
	r1 <- 1
	r2 <- 1
	k1 <- 400
	k2 <- 300
	#######
	# parametros de interaccion de las especies
	#######
	a12 <- 0.9 # 0.9 < 400/300 = 1.333
	a21 <- 0.4 # 0.4 < 300/400 = 0.75
	parametros <- c(r1=r1,r2=r2,k1=k1,k2=k2,a12=a12,a21=a21)
	#######
	# plano de fases
	#######
	plano_fases <- phaseArrows( modelo(parametros), xlims=c(0,k1), ylims=c(0,k2) )
	#######
	# dibujo trayectoria en el plano de fases
	#######
	tmax = 30 # trayectorias entre t=0 y t=tmax
	drawTrajectories( modelo(parametros), tEnd=tmax, x0=c(x=100,y=80) ) # trayectoria con condicion inicial (100,80)
	#######
	# permito dibujar otras n=3 trayectorias indicando el inicio de cada trayectoria con un el mouse
	#######
	n <- 3 # numero de trayectorias a dibujar
	drawTrajectories( modelo(parametros), tEnd=tmax, x0=locator(n,"p") )
	#######
	# curvas isoclinas (nullclines)
	#######
	phaseNullclines(plano_fases$z, plano_fases$xy, add=TRUE)

view raw plano_fases_logistico.R hosted with ❤ by GitHub

Plano de fases

En la primera parte del código se define el modelo dinámico a utilizar.

Una vez que se tiene definido el modelo dinámico y se le han asignado valores a sus parámetros, ya se puede graficar el diagrama de fases. Para eso se utiliza el comando phaseArrows() del paquete "pplane".

Trayectorias

Al gráfico anterior se le puede añadir una trayectoria solución del sistema. Para esto se utiliza el comando drawTrajectories() del paquete "pplane". En el código de ejemplo la población inicial es (100,80).

También es posible indicar con un clik del mouse el punto inicial de la trayectoria mediante el uso del comando drawTrajectories(). Más en general se pueden indicar $n$ puntos iniciales para que dibuje $n$ trayectorias.

Curvas Isóclinas

Las curvas isóclinas (nullclines en inglés) del sistema dinámico son aquellas donde alguna componente del vector de derivadas se anula. Es decir donde $\frac{dx_1}{dt}=0$ o $\frac{dx_2}{dt}=0$.

Es posible agregarlas al plano de fases mediante el comando phaseNullclines(). Notar que en el código al usar este comando utilizo la variable "plano_fases" que definí al dibujar el diagrama de fases con el comando phaseArrows().

En este modelo hay 4 isóclinas que son las dos rectas de pendiente negativa que se ven en el diagrama y los ejes $x_1=0$ y $x_2=0$. El color de las rectas indica si corresponden a la condición $\frac{dx_1}{dt}=0$ o $\frac{dx_2}{dt}=0$.

Apéndice: Instalación del paquete pplane en R

Ejecuto los comandos en una línea de comandos de Debian GNU/Linux.

# indica comando ejecutado como super usuario
$ indica comando ejecutado como usuario normal

Instalo primero un cliente de "subversion":

# aptitude install subversion

Descargo los archivos del paquete "pplane" desde la web del proyecto en R-Forge mediante el siguiente comando:

$ svn checkout svn://r-forge.r-project.org/svnroot/pplane/

Esto crea una carpeta de nombre "pplane" con los archivos del proyecto en el directorio donde se ejecutó el comando.

Finalmente, en el mismo directorio donde se ejecutó el comando anterior, ejecuto el siguiente comando para instalar el paquete en R:

$ R CMD INSTALL --build pplane/pkg/pplane

Si el comando anterior no logra instalar "pplane" e indica que faltan dependencias (en mi caso faltaba el paquete "abind"), hay que entrar a R, instalar estas dependencias con el comando install.packages() y luego volver a ejecutar el comando que instala "pplane".

Referencias

• "Cannot install R-forge package using install.packages"
• Paquete "pplane" en R-Forge