5

He estado recientemente jugueteando con la librería odeint para resolver un sencillo sistema de ecuaciones diferenciales ordinarias.

He conseguido construir un código que parece resolver correctamente el sistema y representar R y J en el eje y, a lo largo del tiempo:

import numpy as np
from scipy.integrate import odeint 
import matplotlib.pyplot as plt

#Parameteros
a=0
b=0.4
c=-0.3
d=0

# Inicialización
tstart= 0
tstop= 50
increment = 1
t = np.arange(tstart,tstop+increment,increment)

ymin, ymax, ystep = -5,5,0.5
y = np.arange(ymin, ymax+ystep, ystep)

y0 = [3.14,-0.5]

# Función que devuelve dy/dt
def mydiff2(y, t):
  dRdt = a*y[0]+ b*y[1]
  dJdt = c*y[0]+ d*y[1]
  
  dydt= [dRdt,dJdt]
  return dydt

# Resolviendo ODE
y = odeint(mydiff2, y0, t)
print(y)
R = y[:,0]
J = y[:,1]

# Graficando los resultados
plt.plot(t,R)
plt.plot(t,J)
plt.title('Simulación con 2 variables')
plt.xlabel('t')
plt.ylabel('y')
plt.grid()
plt.axis([0, 50, -5, 5])
plt.legend(["R", "J"])
plt.show()

Sin embargo, ahora me gustaría representar los campos vectoriales para R y J en un gráfico cuyo eje y sea dy/dt y cuyo eje x sea tiempo. He estado intentando varias cosas pero no consigo representar los campos vectoriales correctamente y de manera legible.

Apreciaría ayuda con esto. Abajo uno de mis intentos.

# Campo vectorial
for y0 in y:
  line = odeint(mydiff2, y0, t)
  plt.plot(t, line, 'b')


x = np.linspace(tstart, tstop, 50)
X, Y = np.meshgrid(x, y)

U = 1
V = mydiff2(Y, None)
N = np.sqrt(U**2 + V**2)
U /= N
V /= N

plt.quiver(X, Y, U, V, angles='xy')
plt.xlabel('tiempo')
plt.ylabel('dy/dt')
plt.axis([tstart, tstop, ymin, ymax])
plt.show()

1 respuesta 1

2
+50

Tienes un sistema dinámico acoplado de dos variables, R y J.

En forma matricial podrías escribir la ecuación diferencial de la siguiente forma:

donde A=[[a,b],[c,d]], x=[[R],[J]] y x'=[[dRdt],[dJdt]]

La solución general ya existe, y es facil de encontrala con los valores propios de la matriz A. (buscar: solución de matriz fundamental)

Código de la pregunta:

# %% Parte 1
# Campo vectorial
for y0 in y:
  line = odeint(mydiff2, y0, t)
  plt.plot(t, line, 'b')

# %% Parte 2
x = np.linspace(tstart, tstop, 50)
X, Y = np.meshgrid(x, y)

U = 1
V = mydiff2(Y, None)
N = np.sqrt(U**2 + V**2)
U /= N
V /= N
# %% Parte 3
plt.quiver(X, Y, U, V, angles='xy')
plt.xlabel('tiempo')
plt.ylabel('dy/dt')
plt.axis([tstart, tstop, ymin, ymax])
plt.show()

La primera parte del código no es necesaria, pues. Para graficar el campo vectorial no es necesario resolver la ecuación diferencial, basta con tomar x'=[[dRdt],[dJdt]] el cual lo obtienes al evaluar la función diferencial de x'=[dRdt,dJdt] <=mydiff2(y, t):

La segunda parte del código no evalúa correctamente la función mydiff2() la cual debe contener como parámetro y=[[R],[J]] ('x' en formula general).

Para cada combinación de entrada (R,J) obtendrás un (dRdt,dJdt) de salida distinto. Notese que no depende del tiempo (buscar sistema lineal invariante en el tiempo)

Para obtener cada punto podrías escribir un ciclo for que ejecute la función para cada punto y obtener correspondientemente el vector del campo. En mi caso busqué la forma de hacer un producto matricial con tamaños de matriz (2x2)(2x100x100)

código:

#%%
import numpy as np
from scipy.integrate import odeint 
import matplotlib.pyplot as plt
#%%

#Parameteros
a=0
b=0.4
c=-0.3
d=0

A=np.array([[a,b],[c,d]])

# Función que devuelve dy/dt
def mydiff2(y, t):
    y=np.array(y)
    [dRdt,dJdt] = np.dot(y.T,A.T).T
    dydt= [dRdt,dJdt]
    return dydt
#%%
# 
Rx = np.linspace(-5,5,100)
Jy = np.linspace(-5,5,100)
RX, JY = np.meshgrid(Rx,Jy)
z=np.array(mydiff2(np.array([RX,JY]),0))
plt.quiver(RX,JY,*z,scale=90)
plt.show()

Gráfica: Espacio de fases Presento como salida una gráfica del "espacio de fases" de eje x->dRdt y eje y->dJdt, porque dado que el sistema dinámico no depende del tiempo no puedes tener una función para graficar a lo largo de t. Podrías graficar la derivada dRdt y dJdt pero no correspondería al campo sino a una solución particular de la familia, mientras el campo es una vista general a todas las familias conociendo la malla que curva el "espacio".

  • En la gráfica podrás observar que es un flujo cortante circular.
  • La solución que encontraste con odeint corresponde a ondas senusoidales
  • la solución analítica es con los valores propios: λ=+-0.34641 i, de forma que es un sistema marginalmente estable (oscila perpetuamente), confirmando el espacio de fases circular y la solución encontrada con odeint.
5
  • Básicamente si tomas las condiciones iniciales y0=[3.14,-0.5], sabiendo que el campo son círculos o los valores propios... puedes saber que la solución es un vector girando de radio = 3.1796 girando a -0.3464 radianes por segundo. Commented el 20 jul. 2021 a las 9:04
  • Muchas gracias @Alfredo Maussa, una respuesta muy muy didáctica e instructiva. Unas preguntas: 1. ¿Es la solución encontrada con odeint lo que llamaríamos solución numérica porque computamos y graficamos un rango de valores? 2. ¿Podrías explicar cómo obtienes los valores propios (eigenvalues) y eigenvectors con Python? 3. ¿Por qué dices que sabiendo las condiciones iniciales podemos saber que la solución analítica es un vector giranto con radio 3.17 a -0.3464 rad/sg. Puedes explicar esto con los cálculos en Python? 4. ¿Cómo podemos saber los puntos fijos de este sistema y si son estables o no?
    – pyring
    Commented el 20 jul. 2021 a las 11:18
  • 1
    1) Sí es una solución numérica porque se utilizan iteraciones con un algoritmo (computable) para aproximar la solución a cada 'paso' de las iteraciones, se debe limitar un rango ya que el computador tiene memoria limitada y no queremos que repita el ciclo iterativo indefinidamente. 2) Aplicando álgebra lineal, en python hay modulos programado, por ejemplo numpy tiene eigenvalue, featurevector = np.linalg.eig(A). 3) Es más una interpretación analítica del espacio de fases (en teoría de control, es equivalente al campo de las variables), (continuo...) Commented el 20 jul. 2021 a las 22:12
  • 1
    3) Te invito a ver este video de 26s: youtu.be/_L0aHxhdeYE El sistema que presentas es similar, de hecho, puedes asumir que R es la velocidad (eje y en el video) y J es el ángulo (eje x en el video). Puedes ver como gira el punto (v,theta) tal cual como lo haría (R,J) en el campo graficado en la respuesta. Al saber que gira circularmente porque el campo es circular... sabemos que la maginitud del vector inicial se conserva radio=sqrt(3.14**2+0.5**2)=3.17, la velocidad de giro viene de que la solución general es de la forma exp(vectorpropio_imaginario) la cual es un 'fasor' (googlear) Commented el 20 jul. 2021 a las 22:19
  • 1
    4) al igual que '3)' también es algo analítico, un punto fijo o estable es aquel que tiende a quedarse quieto, matemáticamente se traduce en que el cambio es 0 (osea las derivadas son 0) gráficamente los puntos estacionarios es donde no hay vectores dibujados (porque son de tamaño nulo), en este caso sería solamente (0,0) lo que significa que el pendulo está quieto/fijo. por otro lado si tuviese cualquier velocidad o ángulo distinto de (0,0) el sistema se movería sí o sí por lo tanto no sería un punto fíjo / estable. en python podrías resolver para dJdt=0 y dRdt =0,lo cual deberia resultar 0,0 Commented el 20 jul. 2021 a las 22:24

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.