I've recently been fiddling with the odeint library to solve a simple system of ordinary differential equations.
I've managed to build code that seems to correctly solve the system and plot R and J on the y-axis, over time:
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()
However, now I would like to plot the vector fields for R and J on a graph whose y-axis is dy/dt and whose x-axis is time. I have been trying various things but I am not able to represent the vector fields correctly and in a readable way.
I would appreciate help with this. Below one of my attempts.
# 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()
You have a coupled dynamical system of two variables, R and J.
In matrix form you could write the differential equation as follows:
where A=[[a,b],[c,d]], x=[[R],[J]] and x'=[[dRdt],[dJdt]]
The general solution already exists, and it is easy to find it with the eigenvalues of matrix A. (search: fundamental matrix solution)
Question code:
The first part of the code is not necessary, then. To graph the vector field it is not necessary to solve the differential equation, it is enough to take x'=[[dRdt],[dJdt]] which you obtain by evaluating the differential function of
x'=[dRdt,dJdt] <=mydiff2(y, t):
The second part of the code does not correctly evaluate the function
mydiff2()
which must contain as a parametery=[[R],[J]]
('x' in general formula).For each input combination (R,J) you will get a different output (dRdt,dJdt). Note that it does not depend on time (look for linear time-invariant system)
To get each point you could write a for loop that executes the function for each point and correspondingly get the field vector. In my case I looked for a way to make a matrix product with matrix sizes (2x2)(2x100x100)
code:
Plot: I output a "phase space" plot of x->dRdt axis and y->dJdt axis, because since the dynamical system does not depend on time you cannot have a function to plot along t. You could graph the derivative dRdt and dJdt but it would not correspond to the field but to a particular solution of the family, while the field is a general view of all the families knowing the mesh that curves the "space".
odeint
corresponds to sine wavesodeint
.