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:
x'=A*x
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.
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:
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.