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()