Queremos graficar el retrato de fase del siguiente sistema no lineal:

$$\ddot{y}+sin(y)=0$$

que en variables de estado se transforma en :

$$\left\{ \begin{array}{lcc} \dot{y}_{1}=y_{2}\\ \\ \dot{y_2}=-sin(y_{1}) \end{array} \right.$$



In [2]:

    
import numpy as np
import matplotlib.pyplot as plt



In [3]:

    
%matplotlib inline



In [4]:

    
def f(Y, t):
    y1, y2 = Y
    return [y2, -np.sin(y1)]



In [5]:

    
y1 = np.linspace(-2.0, 8.0, 20)
y2 = np.linspace(-2.0, 2.0, 20)



In [6]:

    
Y1, Y2 = np.meshgrid(y1, y2)



In [7]:

    
t = 0



In [8]:

    
u, v = np.zeros(Y1.shape), np.zeros(Y2.shape)



In [9]:

    
NI, NJ = Y1.shape



In [10]:

    
for i in range(NI):
    for j in range(NJ):
        x = Y1[i, j]
        y = Y2[i, j]
        yprime = f([x, y], t)
        u[i,j] = yprime[0]
        v[i,j] = yprime[1]



In [11]:

    
plt.quiver(Y1, Y2, u, v, color='r')

plt.xlabel('$y_1$')
plt.ylabel('$y_2$')
plt.xlim([-2, 8])
plt.ylim([-4, 4])
#plt.savefig('phase-portrait.png')









    Out[11]:





(-4, 4)



In [12]:

    
from scipy.integrate import odeint



In [14]:

    
for y20 in [0, 0.5, 1, 1.5, 2, 2.5]:
    tspan = np.linspace(0, 50, 200)
    y0 = [0.0, y20]
    ys = odeint(f, y0, tspan)
    plt.plot(ys[:,0], ys[:,1], 'b-') # path
    plt.plot([ys[0,0]], [ys[0,1]], 'o') # start
    plt.plot([ys[-1,0]], [ys[-1,1]], 's') # end
    
plt.xlim([-2, 8])









    Out[14]:





(-2, 8)



In [ ]: