Todo comenzó porque estaba resolviendo numéricamente una ecuación diferencial ordinaria (EDO para los amigos) y quería ver su comportamiento para diferentes condiciones iniciales. ¿Y cuál es la mejor herramienta para eso? El plano de fase. Así que pinté uno en Matlab y luego, como yo no soy de quedarme satisfecha al terminar algo, me pregunté: ¿cómo se hará esta historieta en Julia?
Después de hurgar en la documentación, dar a todos los botones y pintar cosas sin sentido alguno, lo conseguí con PyPlot, uno de los paquetes que se pueden usar en Julia para hacer gráficos.
Os cuento cómo calculé y cómo pinté. Todo es para Julia 1.0.
Unas notas sobre PyPlot
Antes de meternos en harina, no está de más recordar que Julia permite utilizar diferentes paquetes para dibujar, están entre otros el mítico Plots; GR, que va como un tiro y el que he usado yo esta vez: PyPlot, que es una interfaz para Julia del Matplotlib de Python.
No tengo especial predilección por este paquete, simplemente encontré un ejemplo que se parecía a lo que yo quería hacer y seguí la norma básica de la ingeniería: no reinventar la rueda. Norma especialmente certera cuando eres ingeniera mecánica.
Debo añadir también que PyPlot me funcionó directamente, probablemente porque tenía Miniconda instalado y tal vez Matplotlib también. Ni idea. Estoy todo el rato probando chismes, no sé ya ni lo que hago.
En fin, los gráficos en PyPlot funcionan casi como los de Matlab excepto por un par de cosas que a mí me dejaron loca y os pongo aquí para que no os pase lo mismo:
- Para mostrar la figura actual en una ventana necesitamos hacer
gcf()(get current figure), al menos desde un script. - PyPlot añade por defecto las figuras encima de las que ya había. Antes había opción
hold oncomo en Matlab pero se desaconseja su uso. - A veces no querremos añadir las figuras encima, para borrar anterior tenemos
clf()(clear figure)
Sabiendo esto ya podemos hacer gráficos sin control.
El plano de fase
Para que todo el rollo de la ecuación diferencial sea un poco menos abstracto os voy a dar un poco de contexto: estoy estudiando un sistema de un grado de libertad, es decir, una masa atada a una pared mediante un muelle y un amortiguador.
De esta manera, el plano de fase consiste en un gráfico en el que se representa la posición de la masa en el eje horizontal y su velocidad en el vertical y se pintan unas flechillas que indican la dirección en la que saldrá el sistema al darle unas condiciones iniciales de posición y velocidad.
Nos sirve para ver si aumenta o se reduce su velocidad, tiende hacia una posición de equilibrio o no, o para hacer bonito en un paper, la finalidad última de toda nuestra investigación.
Bien, la ecuación que describe el movimiento del sistema de un grado de libertad es la siguiente:
donde el puntico hace referencia a la derivada respecto al tiempo y dos punticos a la segunda derivada. Los valores m (masa), c (amortiguamiento) y k (rigidez) son características del sistema y son conocidos.
Voy a hacer una magia todo loca y convertir la ecuación de segundo orden en dos ecuaciones de primer orden. La primera va a ser simplemente el cambio de variable:
osease: que la velocidad es la derivada de la posición respecto al tiempo.
La segunda proviene de la ecuación original y viene a decir que la segunda derivada de la posición (la aceleración) es la primera derivada de la velocidad y por lo tanto:
Si cogemos las dos ecuaciones, reorganizamos de tal manera que las cosas con puntico estén a la izquierda y escribimos el sistema en forma matricial tenemos lo siguiente:
Todo este lío nos lleva a algo, paciencia. La cuestión es que queremos ver cómo cambia la posición y la velocidad en cada punto. ¿Y qué cacharro matemático explica cómo cambian las cosas? Exactamente, la derivada. Derivada que precisamente tenemos explícitamente representada en el sistema de ecuaciones.
¿Por qué digo eso? Si yo quiero saber cuánto valen las derivadas de mis variables en el punto (1,0) no tengo más que sustituirlo en la ecuación:
como m, c y k son valores conocidos, puedo obtener los valores de las derivadas. Sabiendo esto, puedo pintar una flechillas en el plano que indiquen hacia dónde tiende el sistema si le doy determinadas condiciones iniciales.
Antes de pintar las flechillas, os dejo aquí la función que he escrito en Julia para calcular las derivadas de la posición y el desplazamiento en diferentes puntos del plano. Como es una función genérica, he utilizado x e y para referirme a las variables y u y v para referirme a los cambios en cada punto y no para lo mismo que en el texto. Es aposta para confundir.
Perdonadme la falta de eficiencia y el probable uso de estructuras ridículas, soy novata.
""" vectorfield(R, [x1 x2 y1 y2]) Calcula el plano de fase de una ecuación diferencial Input - R: matriz de coeficientes - span: límites del campo [x1 x2 y1 y2] Output - x : posición - y : velocidad - u : vector de desplazamiento en un punto (x,y) - v : vector de velocidad en un punto (x,y) """ function vectorfield(R, span) # Span x = range(span[1],span[2],length=20) y = range(span[3],span[4],length=20) # Rejilla de puntos x = repeat(x', length(x)) # solo repite filas y = repeat(y, 1, length(y)) # Inicialización u = zeros(size(x)); v = zeros(size(y)); # Derivada en cada punto for i = 1:size(x,1)*size(x,2) der = R*[x[i]; y[i]]; u[i] = der[1]; v[i] = der[2]; end return (x,y,u,v) end
Trayectoria: Euler progresivo
Antes de que pintemos el gráfico, en plan bonus voy a contar aquí cómo he simulado el comportamiento del sistema en el tiempo con un Euler progresivo, el método más sencillo para resolver una ecuación diferencial numéricamente.
Básicamente consiste en echar una línea desde el primer punto, que conocemos al ser las condiciones iniciales, hasta un poquito más lejos. La línea en cuestión tendrá la pendiente de la derivada en el punto inicial, que conocemos ya que nos la da la ecuación diferencial. Una vez que hemos llegado al segundo punto, tiramos otra línea y así poco a poco vamos avanzando.
Es decir, en el instante de tiempo i+1 x valdrá lo que valía en el instante i más el cambio:
siendo h el paso que damos en el tiempo.
He escrito una función que me calcula la posición y velocidad del sistema en el tiempo. ¡Cuidado si le quitáis el amortiguamiento que se vuelve inestable! Temas numéricos eh, no se volvería inestable en la realidad.
""" euler(R, ic, N, h) Calcula la respuesta a una ecuación diferencial con Euler progresivo Input - R : matrices en Y(t+1) = Y(t) + h*R*Y(t) - ic : condiciones iniciales - N : número de puntos - h : paso Output - x, v : desplazamiento y velocidad del sistema en el tiempo """ function euler(R, ic, N, h) Y = zeros(2,N); Y[:,1] = ic; for t = 1:(N-1) Y[:,t+1] = Y[:,t] + h*R*Y[:,t]; end x = Y[1,:]; v = Y[2,:]; return (x, v) end
El dibujo final
Ahora que ya sabemos cuánto cambian tanto la posición y la velocidad en los puntos del plano y hemos calculado la trayectoria para un caso concreto, vamos a pintar. Para ello he usado la función quiver que, dándole la posición y los componentes en ambas direcciones, pinta unas simpáticas flechillas.
using PyPlot
# Datos del sistema
m = 1
k = 1
cc = 2*sqrt(k*m)
xi = 0.3
c = xi*cc
# Datos de la simulación
N = 1000
ti = 0
tf = 50
t = range(ti,tf,length=N)
h = t[2]-t[1]
ic = [1 0]'
R = [0 1; -k/m -c/m]
F = zeros(2,N)
# Trayectoria
x,v = euler(R, F, ic, N, h)
clf()
plot(x,v)
gcf()
# Plano de fase
span = [-1 1 -1 1]
x,y,u,v = vectorfield(R, span)
q = quiver(x,y,u,v,color="gray")
xlabel("x")
ylabel("v")
gcf()
Así, nos aparece este gráfico:
Aparte de ser muy cuco, este gráfico nos indica que si le damos cualquier tipo de condición inicial a nuestro sistema, tenderá a pararse: todas las flechitas llevan al punto (0,0). O como diría yo si me quiero hacer la interesante: las raíces de este sistema son complejas con parte real negativa.
Y hasta aquí hemos llegado, os dejo una pregunta de deberes: ¿qué pinta creéis que tendría el plano de fase si le quitásemos el amortiguamiento? Contadme vuestras sospechas.
Referencias
Os dejo con la canción más bonita del mundo.
