import numpy as np                       # llama a la libreria numpy (operaciones numericas)
import matplotlib.pyplot as plt          # llama a la libreria grafica

# Este codigo muestra los problenas de estabilidad de la integracion por Euler explicito.
#
# Vamos a definir una funcion maricial que realiza la siguente cuenta
#
# [0    1].[  y(t)  ]
# [-w^2 0] [dy(t)/dt]
#

tita0=0.1
omega2=2.0

n=80
h=(2*np.pi)/n
f=1.0/(1.0+h*h*omega2)                          # factor para Euler implicito

n=2000

y=np.linspace(tita0,tita0,n)
ydot=np.linspace(0,0,n)

yimp=np.linspace(tita0,tita0,n)
ydotimp=np.linspace(0,0,n)

i=1


while i<n:
	yimp[i]=f*yimp[i-1]+f*h*ydotimp[i-1]
	ydotimp[i]=-f*omega2*h*yimp[i-1]+f*ydotimp[i-1]
	y[i]=y[i-1]+h*ydot[i-1]
	ydot[i]=-omega2*h*y[i-1]+f*ydot[i-1]
	i=i+1

plt.figure()                                    # inicio nuevo grafico
#plt.plot(yimp,ydotimp,color='steelblue')       # grafico de diagrama de fases
#plt.plot(y,ydot,color='darkorange')            # grafico de diagrama de fases
plt.plot(y,color='darkorange')                  # grafico y vs. t e ydot vs. t
plt.plot(yimp,color='steelblue')                # grafico y vs. t e ydot vs. t
plt.show()                                      # mostrar en pantalla