El modelo de crecimiento logístico fue publicado por primera vez por Pierre François Verhulst en 1838 y lo utilizó para describir el crecimiento auto-limitado de una población biológica. Las suposiciones básicas de este modelo son Blanchard et al., 1999:

• Si la población es pequeña la tasa de crecimiento de la población es proporcional a su tamaño

• Si la población es demasiado grande para ser soportada por su entorno y recursos, el tamaño de la población disminuye

Un modelo matemático que recoge estas observaciones está definido por el siguiente problema de valor inicial:

$\frac{dP}{dt}=aP(1-\frac{P}{K})$

$P\left(t_0\right)=P_0$

En él, a>0 y K>0 son parámetros característicos del modelo; el valor de a se conoce como la tasa intrínseca de crecimiento, el valor k es la capacidad de carga del sistema y representa la población máxima con la que el medio se satura.

Autores como Blanchard et al., 1999 afirman que existen poblaciones en las que si el número de individuos es elevado, la razón de cambio decrece y puede llegar a ser negativa y si el tamaño de la población es demasiado pequeño esta razón también decrece. Un modelo que incorpora estas observaciones es el siguiente:

$\frac{dP}{dt}=aP(\frac{P}{L}-1)(1-\frac{P}{K})$

en el cual L es una constante que se conoce como el factor de escasez y 0<L<K.

De acuerdo con Blanchard et al., 1999, la ardilla negra es un pequeño mamífero nativo de las montañas rocallosas cuya población es muy territorial y satisface el modelo de crecimiento logístico modificado; utilizaremos el paquete deSolve para analizar el comportamiento de las soluciones del modelo para lo cual, debido a su territorialidad, consideramos una taza intrínseca de crecimiento moderada con valor a=1.5, L=6 como valor del factor de escasez y K=14 como la capacidad de carga del sistema. Adicionalmente consideramos como tamaños iniciales de la población, los valores P(0)=3,P(0)=5,P(0)=7,P(0)=9,P(0)=15 y P(0)=20 que nos permiten sugerir conclusiones acerca del comportamiento de las soluciones del modelo.

La implementación en R.

En esta sección se desarrollan en R las etapas descritas en la sección 2, para resolver numéricamente la modificación de la ecuación logística con las condiciones iniciales consideradas anteriormente para lo cual utilizamos la función ode() del paquete deSolve, así mismo inicialmente se cargan las librerías internas del paquete.

Una vez cargado el paquete deSolve, las siguientes líneas muestran cómo se cargan sus librerías y se definen los parámetros del modelo.

Fig2

Elaboración propia

# Figura 2 # install.packages("deSolve") library(deSolve) # carga de las librerias # La explicitación del modelo # Vector que contiene los valores de los parámetros params1=c(a=1.5,L=6, K=14)

Aquí definimos las condiciones iniciales para el modelo, cada una de ellas como un vector.

CondIn

Elaboración propia

# Condiciones iniciales cini1=c(3) cini2=c(5) cini3=c(7) cini4=c(9) cini5=c(15) cini6=c(20)

La ecuación del modelo se incluye a través de una función con nombre Logis que calcula la razón de cambio de la variable de estado P, la devuelve como una lista y usa como argumentos los tiempos de corrida, la variable de estado y los parámetros.

Ecuaciones

Elaboración propia

# Ecuaciones # P es la variable de estado Logis= function(t, P, parms) { Salida1=with(as.list(c(P, parms)),{ dP=a*P*(P/L-1)*(1-P/K) list(c(dP))}) return(Salida1)}

El modelo se corre en el intervalo de tiempo [0,2] con paso 0.1 y se resuelve utilizando la función ode() que contiene la rutina lsoda como método de integración predeterminado y que tiene como argumentos básicos las condiciones iniciales, el tiempo de corrida, la función Logis y los parámetros.

Elaboración propia

# La carga del modelo # Tiempos tiempos1=seq(0,2,0.1) # La integración utiliza la función ode. sol1=ode(y=cini1, t=tiempos1, func=Logis, parms=params1) sol2=ode(y=cini2, t=tiempos1, func=Logis, parms=params1) sol3=ode(y=cini3, t=tiempos1, func=Logis, parms=params1) sol4=ode(y=cini4, t=tiempos1, func=Logis, parms=params1) sol5=ode(y=cini5, t=tiempos1, func=Logis, parms=params1) sol6=ode(y=cini6, t=tiempos1, func=Logis, parms=params1)

La función ode() tiene como salida un objeto de la clase deSolve que es una matriz que contiene los valores de la variable de estado en los tiempos proporcionados. Las primeras 10 iteraciones de este proceso, para la solución, sol1, se obtienen con:

# La representación tabular head(sol1, n=10)

Lo que genera una lista con los valores de los tiempos y el correspondiente valor de la variable de estado que se muestran en la tabla 3.

Tabla 3
Primeras 10 iteraciones

Elaboración propia

Las gráficas de las seis soluciones y algunos textos ilustrativos, se obtienen con el siguiente código y se muestran en la figura 2.

RepGraf

Elaboración propia

# La representación gráfica plot(sol1,sol2,sol3, sol4,sol5, sol6, lwd=2, xlab="t", ylab="P(t)", col=c("red", "green", "blue", "deeppink3", "brown4", "orange"), main= "Crecimiento logístico modificado") legend("topright",legend=c("P(0)=3", "P(0)=5", "P(0)=7","P(0)=9", "P(0)=15","P(0)=20"), col=c("red", "green","blue", "deeppink3", "brown4", "orange"), lwd=c(2,2), bg="aliceblue")

Figura 2
Simulaciones numéricas del modelo logístico modificado
Elaboración propia

Esta gráfica sugiere que si el tamaño de la población inicial de ardillas negras es menor que la dada por el factor de escasez la población tiende a extinguirse, que si la población de ardillas negras está entre el factor de escasez y la capacidad de soporte el tamaño de la población crece acercándose a la capacidad de carga y que si el tamaño de la población es mayor que la capacidad de carga esta disminuye acercándose a este valor.

Modificamos el código anterior para incluir las soluciones constantes, P=0, P=6 y P=14, denominadas Puntos de equilibrio, así como otras condiciones iniciales.

Fig3a

Elaboración propia

fig3b

Elaboración propia

# figura 3 # install.packages("deSolve") library(deSolve) # Solución de la ecuación logística tiempos2=seq(from = 0, to = 3, by = 0.01) param2 = c(a=1.5, L=6, K=14) logis = function(t,y,param){ salida = with(as.list(c(y,param)),{ dP=a*P*(P/L-1)*(1-P/K) list(c(dP))}) return(salida)}

oncecur=seq(0,20,2)

listadatos=lapply(oncecur,function(i){ode(y=c(P=i), times=tiempos2,func=logis,parms = param2)})

curvas= rbind(listadatos[[1]],listadatos[[2]],listadatos[[3]], listadatos[[4]],listadatos[[5]],listadatos[[6]], listadatos[[7]],listadatos[[8]],listadatos[[9]], listadatos[[10]],listadatos[[11]])

curvas2=data.frame(curvas) curvas2$group=c(rep("P(0)=0",301),rep("P(0)=2",301), rep("P(0)=4",301),rep("P(0)=6",301),rep("P(0)=8",301), rep("P(0)=10",301),rep("P(0)=12",301),rep("P(0)=14",301), rep("P(0)=16",301),rep("P(0)=18",301),rep("P(0)=20",301)) curvas2

logis2=curvas2 # Modificar los datos data_label = logis2 data_label$label = NA data_label$label[which(data_label$time == min(data_label$time))] <- data_label$group[which(data_label$time == min(data_label$time))]

install.packages("ggplot2") require(ggplot2)

colcur = c("red","#999999","#E69F00", "red", "deeppink4", "green", "#000000", "red", "#0072B2","#009E73" , "#CC79A7")

install.packages("ggrepel") require(ggrepel)

# ggplot2 traza el gráfico con etiquetas ggplot(data_label, aes(time, P, group = forcats::fct_inorder(group))) + geom_line(aes(colour= forcats::fct_inorder(group)),size=0.8) + scale_colour_manual(values=cbp2)+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_rect(colour = "black", size=0.5, fill=NA))+ geom_label_repel(aes(label = label,color=group), nudge_x = -0.15, size=2, label.size=NA, fill="honeydew2", na.rm = TRUE,segment.color=NA) + theme(legend.position = "none") + labs(x = "t", y="P(t)")

El proceso anterior genera la gráfica que se muestra en la figura 3. Observe que todas las soluciones cercanas a la solución P=6 se ”alejan" de ella y que, en general, las soluciones "tienden" hacia la solución de equilibrio P=0 o P=14, es decir, P=6 es una fuente y P=0, P=14 son sumideros.

Figura 3
Simulaciones numéricas del modelo logístico modificado
Elaboración propia

El sistema de Chen

El concepto de atractor extraño surge a raíz de los trabajos de E. Lorenz, quien en 1963, escribió un notable artículo en el cual en un intento, por estudiar la posibilidad de predecir el comportamiento del flujo atmosférico a largo plazo, es decir, conocido el estado del flujo atmosférico en un instante dado determinar su comportamiento en tiempos futuros distantes, simuló experimentalmente el proceso convectivo atmosférico y observó que el flujo resultante podía ser muy irregular con trayectorias no periódicas inestables. Matemáticamente este proceso puede ser descrito utilizando un modelo de Rayleigh - Benard el cual es un sistema de ecuaciones en derivadas parciales que se deduce de las ecuaciones de Navier - Stokes para el que no se conocen soluciones exactas a excepción de las triviales. Suponiendo una expansión de Fourier de cierto tipo Lorenz reduce el sistema de ecuaciones en derivadas parciales a un sistema autónomo de ecuaciones diferenciales ordinarias que se conoce como el sistema de Lorenz Mosquera, 1992.

El Sistema de Lorenz posee dependencia sensitiva con respecto a las condiciones iniciales, consecuentemente exhibe comportamiento caótico y muestra lo que se conoce como un atractor extraño que corresponde, en su diagrama de fase, a una estructura geométrica con forma poco usual en la que las las trayectorias del sistema nunca se cortan y son líneas de longitud infinita encerradas en un área finita, que describen órbitas no periódicas Mosquera, 1992.

A partir del trabajo de Lorenz se han propuesto varios sistemas de ecuaciones diferenciales que contienen atractores extraños, entre ellos, el atractor de Rossller, el atractor de Chen, el atractor de Chua, el atractor de Liu y otros. En esta sección implementaremos en R el sistema de ecuaciones diferenciales ordinarias que produce el atractor de Chen.

El sistema de Chen Sooraksa y Chen, 2018 fue propuesto en 1999 y puede considerarse como un modelo de clima controlado bajo la perspectiva anticontrol. Desde el punto de vista matemático, es el dual del sistema de Lorenz a través de la inversión en el tiempo y corresponde a un sistema autónomo de ecuaciones diferenciales ordinarias definido por

$\frac{dx}{dt}=-ax+ay$

$\frac{dy}{dt}=(c-a)x+cy-xz$

$\frac{dz}{dt}=xy-bz$

Las variables x y z representan, el promedio espacial de la velocidad hidrodinámica y la diferencia de temperatura entre la corriente ascendente y la descendente y la variable y el gradiente de temperatura. a, b y c son parámetros característicos del sistema y dependiendo de sus valores se tienen diversas situaciones; por ejemplo, de acuerdo a Sooraksa y Chen, 2018 con los valores de los parámetros a=35, b=3 y c=28 y las condiciones iniciales x(0)=-3, y(0)=2, z(0)=20 el sistema presenta comportamiento impredecible y se obtiene un conjunto de estructura geométrica caótica que se conoce cómo el atractor de Chen.

La implementación en R.

El sistema de Chen presenta una variada dinámica, sin embargo, presentamos el caso más importante de esta, su atractor extraño.

El siguiente código implementa la descripción realizada en la sección 2, para resolver numéricamente el Sistema de Chen, para lo cual nuevamente se utiliza la función ode() del paquete deSolve. Siguiendo a Sooraksa y Chen, 2018 consideramos los valores de los parámetros a=35, b=3 y c=28 para los cuales se produce el atractor de Chen.

Con base en lo expuesto en K. Soetaert et al., 2016 para construir un gráfico de los resultados numéricos obtenidos es posible utilizar un método diseñado explícitamente para objetos de la clase deSolve, que ordena las figuras en dos filas y dos columnas. En este caso en ellas se presentan las gráficas de las variables dependientes contra la independiente y puesto que queda un espacio en el margen inferior derecho en él, se muestra el gráfico del promedio espacial de la velocidad hidrodinámica contra la diferencia de temperatura entre la corriente ascendente y la descendente, que en este caso corresponde a la proyección del atractor de Chen en el plano xz.

fig4

Elaboración propia

# Figura 4 # El Sistema de Chen # Las ecuaciones del modelo Chen = function (t, y, parms) { with(as.list(y), { dX = -a * X+a*Y dY = (c-a)*X+c*Y - X*Z dZ =X*Y - b*Z list(c(dX, dY, dZ)) })}

# Valores de los parámetros y condiciones iniciales param=c(a = 35, b = 3, c = 28) Eini = c(X = -3, Y = 2, Z = 20)

# la integración se realiza por 100 veces cada 0.001 vez library(deSolve) tiempos = seq(from = 0, to = 100, by = 0.001) salida = ode(y = Eini, times = tiempos, func = Chen, parms = param)

# La gráfica del modelo plot(salida,xlab="tiempo", col= "blue", lwd = 1) plot(salida[,"X"], salida[,"Z"], col= "red", type = "l", xlab = "X", ylab = "Z", main = "Proyección plano XZ")

La gráfica resultante de este proceso se muestra en la figura 4.

Figura 4
Simulaciones numéricas del sistema de Chen
Elaboración propia

Las siguientes líneas de R, generan una gráfica 3D que ilustra el denominado Atractor de Chen y que se muestra en la figura 5.

Figura 5
Una vista 3d del atractor de Chen
Elaboración propia

# Figura 5 # gráfica 3D library(scatterplot3d) scatterplot3d(salida[,-1], type = "l", lwd = 1, xlab = "X", ylab = "Y", zlab = "Z", main = "El atractor de Chen", col.main="red",color=("darkorchid"),box=FALSE, angle=45)

En esta sección se ha presentado únicamente una muestra del variado comportamiento que presenta este sistema para diversos valores de sus parámetros, en Sooraksa y Chen, 2018 puede encontrarse una descripción más completa de este sistema.

El paquete phaseR

El paquete phaseR fue desarrollado, alrededor del año 2012, por Michael J. Grayling, Gerhard Burger, Stephen P. Ellner y John M. Guckenheimer y es un paquete R para el análisis cualitativo de sistemas autónomos de ecuaciones diferenciales ordinarias en una y dos dimensiones, para lo cual utiliza el método del plano fase. En la tabla 4 se muestra una descripción de las funciones incorporadas en este paquete Grayling, 2022.

Tabla 4
Descripción de las funciones del paquete phaseR

Elaboración propia

Su instalación se realiza de la manera usual, internamente el paquete utiliza la función ode() y el propósito central de esta sección es la de utilizar algunas de las funciones descritas en la tabla 4 para ilustrar las potencialidades del paquete para tratar desde el punto de vista cualitativo, problemas de valor inicial definidos para sistemas autónomos en $ℝ$ y en ${ℝ}^2$

De acuerdo a Grayling, 2022 presentamos la sintaxis básica de cada una de las funciones del paquete phaseR() que utilizaremos en este documento.

Para la función flowField la sintaxis es:

flowField(func, xlim, ylim, parameters, points, system, add,...)

y una descripción de sus argumentos se presenta en la tabla 5

Tabla 5
Los argumentos de flowField() y su descripción

Elaboración propia

Para la función nullclines la sintaxis es:

nullclines(func, xlim, ylim, parameters, points, system, add,...)

en ella se observa que su sintaxis coincide con la de la función flowFields en consecuencia la interpretación de sus argumentos es exactamente la misma.

Para la función numericalSolution la sintaxis es:

numericalSolution(func, y0, tlim, type, parameters, col,...)

y la descripción de sus argumentos se presenta en la tabla 6

Tabla 6
Los argumentos de numericalSolution() y su descripción

Elaboración propia

Para la función trajectory la sintaxis es:

trajectory(func, y0, t.lim, parameters, system, col,...)

y una descripción de sus argumentos se presenta en la tabla 7.

Los argumentos trajectory() y su descripción


Elaboración propia

Para la función stability la sintaxis es:

stability(func, ystar, parameters, system, summary,...)

y una descripción de sus argumentos se presenta en la tabla 8

Tabla 8
Los argumentos de stability() y su descripción

Elaboración propia

Cada una de estas funciones posee argumentos adicionales que no son necesarios para nuestros objetivos sin embargo una descripción de estos puede verse en Grayling, 2022.

El siguiente ejemplo muestra de manera sintética cómo utilizar algunas de estas funciones para realizar un análisis cualitativo del siguiente problema.

Ejemplo 2. Utilizar R para generar el campo direccional, construir las isoclinas nulas, trayectorias con diversas condiciones iniciales y un análisis de estabilidad para los puntos de equilibrio de ecuación diferencial

$\frac{dP}{dt}=aP(\frac{P}{L}-1)(1-\frac{P}{K})$

$P\left(0\right)=P_0$

Este es el modelo de la ecuación logística modificada considerada en la subsección 3.1 por lo que como valores de los parámetros consideramos los seleccionados en aquella sección. El código R que genera lo solicitado, es el siguiente

fig6

Elaboración propia

# figura 6 # install.packages("phaseR") library(phaseR) # Análisis cualitativo de la ecuación logística modificada tiempos2=seq(from = 0, to = 5, by = 0.01) param2 = c(a=1.5, L=6, K=14) logis.mod = function(t,P,param){ salida = with(as.list(c(P,param)),{ dP=a*P*(P/L-1)*(1-P/K) list(c(dP))}) return(salida)}

# El campo vectorial logis.mod.flowField = flowField(logis.mod, xlim = c(0, 5), ylim = c(-1, 20), parameters = param2, points = 21, system = "one.dim", add = FALSE, xlab = "t", ylab = "P(t)") grid()

# Las isoclinas nulas logis.mod.nullclines = nullclines(logis.mod, xlim = c(0, 5), ylim = c(-1, 20), parameters = param2, system = "one.dim")

# Trayectorias con diversas condiciones iniciales logis.mod.trajectory = trajectory(logis.mod, y0 = c(1.5, 2.5,7.5, 11.5,15, 18 ), tlim = c(0, 5), parameters = param2, system = "one.dim", col = rep("red", 6))

# Estabilidad de los puntos de equilibrio logis.mod.Estab1=stability(logis.mod,ystar = c(0),parameters = param2, system = "one.dim") logis.mod.Estab2=stability(logis.mod,ystar = c(6),parameters = param2, system = "one.dim") logis.mod.Estab3=stability(logis.mod,ystar = c(14),parameters = param2, system = "one.dim")

El código anterior genera la gráfica de la figura 6 que muestra el campo direccional, las isoclinas nulas y trayectorias para algunas condiciones iniciales.

Figura 6
Campo vectorial, isoclinas nulas, trayectorias de la ecuación logística modificada
Elaboración propia

En cuanto a estabilidad de los puntos de equilibrio los resultados son los siguientes

discriminant = -1.5, classification = Stable discriminant = 0.85714, classification = Unstable discriminant = -2, classification = Stable

y que por tanto clasifican a P=6 como una el punto de equilibrio inestable y a P=0 y P=14 como puntos de equilibrio asintóticamente estables.

La Ley de Hooke fue formulada en 1768 por Robert Hooke y básicamente establece que:

"La fuerza que devuelve un resorte a su posición de equilibrio es proporcional a la distancia que se desplaza de esta posición" (Cerón & Guerrero, 2008).

Cuando este enunciado se expresa matemáticamente produce una ecuación diferencial lineal de segundo orden que se estudia, usualmente, en un curso básico de ecuaciones diferenciales ordinarias. En este aparte se analiza la ecuación diferencial no lineal

$m\frac{d^{2}x}{d{t}^2}+kx+a{x}^3=0$

que describe el movimiento de un cuerpo de masa m que cuelga de un resorte y fue propuesta en Cerón & Guerrero, 2008 como un modelo que generaliza la ley de Hooke. Físicamente x representa el desplazamiento de un cuerpo de masa m en el extremo del resorte no lineal sin fricción, k es un parámetro positivo que se denomina elasticidad del resorte y cuya rigidez $a\ne 0$ varía con el desplazamiento. Si a>0 la rigidez aumenta con el desplazamiento y se tiene un resorte duro, pero para a<0 la rigidez disminuye con el desplazamiento y se tiene un resorte suave (Cerón & Guerrero, 2008).

La utilización de los paquetes deSolve y PhaseR para analizar esta ecuación diferencial requiere transformarla en un sistema dos-dimensional de ecuaciones diferenciales; para lo cual utilizamos el cambio de variable $\frac{dx}{dt}=y$ con lo que se obtiene el sistema no lineal bidimensional

$\frac{dx}{dt}=y$

$\frac{dy}{dt}=-\frac{k}{m}x-\frac{a}{m}{x}^3$

en el cual, llamando $c^2=\frac{k}{m}$ y $b=\frac{a}{2k}\ne 0$ se puede escribir como

$\frac{dx}{dt}=y$

$\frac{dy}{dt}=-c^2(x+2b{x}^3)$

El resorte duro

Cerón at el (2008) demuestran que, en el caso del resorte duro, el comportamiento de las trayectorias es independiente del valor de b, por ello para considerar este caso, consideramos los valores $c^2=1$ y $b=2$ e implementamos el siguiente código, que utiliza la función ode y genera las gráficas del desplazamiento y la velocidad en ventanas diferentes.

fig7

Elaboración propia

# Figura 7 library(deSolve) #Definimos las condiciones iniciales Cini = c(X=1, Y=0) # Definimos los tiempos en los cuales queremos construir la solución tiempos = seq(from=0, to=8, by=0.01) # Definimos los valores de los parámetros del modelo params = c(c=1, b=2) # Definimos la función que expresa nuestro sistema de ecuaciones hooke1= function(t,y,parms){ Salida = with(as.list(c(y,parms)),{ dX = Y dY = -c^2*(X+2*b*X^3) list(c(dX,dY)) }) return(Salida) } # Se llama la función ode para construir la solución Salida= ode(y = Cini, times = tiempos, func = hooke1, parms = params) # Se construye la gráfica de la solución plot(Salida,xlab="tiempo", col="red",type="l")

Este proceso que proporciona la imagen que se muestra en la figura 7

Figura 7
A la izquierda, posición y a la derecha, velocidad
Elaboración propia

Este gráfico no muestra toda la interacción del sistema, consecuentemente construimos un gráfico más acorde a nuestras necesidades para lo cual generamos, en primer lugar, el gráfico del desplazamiento del resorte y en la misma ventana el de su velocidad.

fig8

Elaboración propia

# Figura 8 par(mfrow=c(1,1)) # Posicion X contra t plot(Salida[,"time"],Salida[,"X"], ylim=c(-2, max(Salida[ ,"X"], 3)), type="l", col="red", lwd=2, xlab="Tiempo", ylab="", main="Un modelo que generaliza la ley de Hooke") grid() # velocidad Y contra t lines(Salida[,"time"], Salida[,"Y"], type="l", col="blue", lwd=2) legend("topright", legend=c("Desplazamiento", "Velocidad"), col=c("red", "blue"), lwd=c(2,2), bg="snow1")

Este proceso que genera la gráfica de la figura 8, la cual sugiere que las soluciones de la ecuación diferencial tienen carácter periódico cómo se verifica teóricamente en Cerón & Guerrero, 2008.

Figura 8
Desplazamiento y velocidad contra tiempo
Elaboración propia

Una manera diferente de visualizar la dinámica de un sistema autónomo bidimensional, es utilizar el gráfico de las trayectorias en el plano fase, es decir, un gráfico del desplazamiento contra la velocidad, lo cual en R se puede obtener utilizando algunas de las funciones del paquete phaseR, ilustramos este hecho a continuación, para lo cual en primera instancia instalamos el correspondiente paquete y activamos sus librerías.

install.package{"phaseR"} library(phaseR)

En el caso del resorte duro se demuestra en Cerón & Guerrero, 2008 que el sistema tiene un único punto de equilibrio en el punto (0,0) y que las trayectorias del sistema corresponden a órbitas cerradas para todo b>0, consecuente con este hecho para el análisis de la ecuación seguimos considerando como valores de los parámetros c=1 y b=2.

El siguiente código nos ilustra, de una manera diferente, el comportamiento del sistema en el plano fase, para el cual graficamos trayectorias con diferentes condiciones iniciales, las isoclinas nulas y el campo direccional.

fig9

Elaboración propia

# Figura 9 library(phaseR) Hooke=function(t, y, param){ # Variables x =y[1] y =y[2] # Parámetros c =param[1] b =param[2] # Ecuaciones diferenciales dy =numeric(2) dy[1]=y dy[2]= -c^2*(x+2*b*x^3) list(dy)} # Campo vectorial Hooke.flowField= flowField(Hooke, xlim = c(-3, 5), ylim = c(-12, 12.5), param = c(1, 2), points =21 , add = FALSE) # Isoclinas nulas Hooke.nullclines = nullclines(Hooke, xlim = c(-5, 5), ylim = c(-12, 12.5), param = c(1, 2), points = 500) Cini = matrix(c(1, 2, 1.5, 2, 2.5, 2, 3,5), ncol = 2, nrow = 4, byrow = TRUE) # Curvas solución Hooke.trajectory = trajectory(Hooke, y0 = Cini, tlim= c(0,10), param = c(1, 2), col = rep("red", 4))

En este proceso produce la gráfica que ilustra la figura 9 que verifica que para esta condición inicial y estos valores de los parámetros el sistema posee órbitas cerradas.

Figura 9
Campo direccional, isoclinas nulas y trayectorias, resorte duro
Elaboración propia

El siguiente script utiliza la función numericalSolution para graficar en una misma ventana, el desplazamiento y la velocidad contra el tiempo. La gráfica resultante se muestra en la figura 10.

fig10

Elaboración propia

# Figura 10 Hooke.num.sol= numericalSolution(Hooke, y0 = c(1, 2), tlim=c(0,10), type = "one", parameters = c(1, 2), col = c("red", "green"), ylab = "X, Y", ylim = c(-3, 3)) legend("topright",legend=c("X", "Y"),col=c("red", "green"), lwd=c(2,2), bg="aliceblue")

Figura 10
Desplazamiento y Velocidad contra el tiempo, resorte duro.
Elaboración propia

Para finalizar esta sección utilizamos la función stability, de este paquete, para determinar, la estabilidad del punto de equilibrio en el origen, lo que obtenemos con la siguiente línea.

Hooke.estabilidad.d=stability(Hooke,ystar=c(0,0),parameters=c(1,2))

R genera como respuesta la línea siguiente, que nos dice que el punto de equilibrio (0,0) corresponde a un centro como se verifica teóricamente en Cerón & Guerrero, 2008.

tr = 0, Delta = 1, discriminant = -4, classification = Centre

El resorte suave.

En este caso el sistema posee tres puntos de equilibrio, $\left(0,0\right)$, $(\frac{1}{\sqrt{-2b}},0)$ y $(\frac{-1}{\sqrt{-2b}},0)$ y puesto que en Cerón at el (2008) se demuestra que el sistema posee órbitas cerradas en la región

$D=\{(x,y)||x|<\frac{1}{\sqrt{-2b}},y\in ℝ\}$

para el análisis de la ecuación consideramos como valores de los parámetros c=1 y b=-2 e inicialmente una condición inicial en el interior de esta región. En el primer script utilizamos, únicamente con fines ilustrativos, la función rk la cual usa directamente la rutina de integración Runge-Kutta5 como un método de integración por defecto y proporciona, en ventanas diferentes, las gráficas del desplazamiento y la velocidad contra el tiempo.

fig11

Elaboración propia

# Figura 11 library(deSolve) Hooke2 = function(t, x, parms) { Salida1=with(as.list(c(parms, x)),{ dX = Y dY = -c^2*(X+2*b*X^3) list(c(dX, dY)) }) return(Salida1) } tiempos2 = seq(0, 20, 0.01) params = c(c = 1, b= -2) Cini2 =c(X = 0.3, Y= 0.2) # La función rk usa ode45 como un método por defecto Solucion= rk(y = Cini2, times = tiempos2, func = Hooke2, parms = params) plot(Solucion, xlab="tiempo",col="red",type="l")

Este código produce la gráfica que se ilustra en la figura 11. Nótese que hemos utilizado como condición inicial el punto (0.3, 0.2) como valores de los parámetros c=1 y b=-2 y que esta gráfica sugiere la existencia en el sistema de trayectorias periódicas.

Figura 11
Posición y Velocidad, resorte suave, Cini= (0.3, 0.2)
Elaboración propia

Si conservamos los valores de los parámetros, modificamos la condición inicial a (0.7, 0.5), el tiempo al intervalo [0, 1.2] y corremos nuevamente el código anterior, obtenemos la gráfica que se ilustra en la figura 12, la cual nos pone de manifiesto la posibilidad de la existencia de trayectorias no periódicas.

Figura 12
Posición y Velocidad, resorte suave, Cini=(0.7, 0.5)
Elaboración propia

Si en el caso del resorte suave, utilizamos las funciones trajectory, flowField y nullclines para generar trayectorias con diversas condiciones iniciales, el campo vectorial y las isoclinas nulas, el código que genera la figura 13, es el siguiente.

fig13

Elaboración propia

# Figura 13 library(phaseR) Hooke=function(t, y, parameters){ #variables x =y[1] y =y[2] # parametros c =parameters[1] b =parameters[2] # definimos las ecuaciones diferenciales dy =numeric(2) dy[1]=y dy[2]= -c^2*(x+2*b*x^3) list(dy)} Hooke.flowField= flowField(Hooke, xlim = c(-2, 2), ylim = c(-2, 2), parameters = c(1, -2), points = 21, add = FALSE) Hooke.nullclines= nullclines(Hooke, xlim = c(-2, 2), ylim = c(-2, 2), parameters = c(1, -2), points = 500) Cini=matrix(c(0.1,0, 0.3,0, 0.2,0.1, 0.3, 0.2, -0.2, -0.3,0.2, 0.3,-1,1, 1,-1,-0.8,0.8,0.8,-0.8, -1,1.5,1,-1.5,-1.5,1.5,1.5,-1.5), ncol = 2, nrow = 14, byrow = TRUE) Hooke.trajectory = trajectory(Hooke, y0 = Cini, tlim= c(0,10), parameters = c(1, -2), col = rep("mediumvioletred", 6))

Figura 13
Campo direccional, isoclinas nulas y trayectorias, resorte suave.
Elaboración propia

El siguiente script, genera la gráfica de la figura 14, y codifica dos situaciones. En primer lugar, el desplazamiento y la velocidad contra el tiempo, para la condición inicial (0.3, 0.2) la cual está en una región en la cual existen órbitas periódicas y en segundo lugar, el desplazamiento y la velocidad contra el tiempo, para la condición inicial (0.7, 0.5) que está en una región en la que no existen órbitas periódicas. El script utiliza la función numericalSolution del paquete phaseR.

fig14

Elaboración propia

# figura 14 par(mfrow=c(1,2)) Hooke.num.sol= numericalSolution(Hooke, y0 = c(0.3, 0.2), tlim=c(0,15),type = "one", parameters = c(1, -2), col = c("red", "green"), ylab = "X, Y", xlab = "tiempo", ylim = c(-0.5, 0.5)) legend("topright",legend=c("X", "Y"),col=c("red", "green"), lwd=c(2,2), bg="aliceblue")

Hooke.num.sol= numericalSolution(Hooke, y0 = c(0.7, 0.5), tlim=c(0,1.5),type = "one", parameters = c(1, -2), col = c("red", "green"), ylab = "X, Y", xlab = "tiempo", ylim = c(0, 10)) legend("topright",legend=c("X", "Y"),col=c("red", "green"), lwd=c(2,2), bg="aliceblue")

Figura 14
Desplazamiento y Velocidad contra el tiempo. Resorte suave: izquierda Cini=(0.3, 0.2), derecha Cini=(0.7, 0.5)
Elaboración propia

Estas gráficas que nos ilustran, en el caso del resorte suave, la diversidad de comportamiento que pueden tener las trayectorias del modelo que generaliza la ley de Hooke.

Cerramos esta sección utilizando la función stability para determinar la estabilidad de los puntos de equilibrio (0, 0), (1, 0) y (-1, 0).

Hooke.Estab1=stability(Hooke,ystar = c(0, 0),parameters = c(1, -2)) Hooke.Estab2=stability(Hooke,ystar = c(1, 0),parameters = c(1, -2)) Hooke.Estab3=stability(Hooke,ystar = c(-1,0),parameters = c(1, -2)) líneas de código que dan como resultado

tr = 0, Delta = 1, discriminant = -4, classification = Centre tr = 0, Delta = -11, discriminant = 44, classification = Saddle tr = 0, Delta = -11, discriminant = 44, classification = Saddle

y que por tanto clasifican, el punto de equilibrio (0, 0) como un centro y los puntos de equilibrio (1, 0) y (-1, 0) como puntos silla.