# Importamos librerias necesarias
# SciPy, SciPy,integrate: librerías para solución por Runge-Kutta 4,5 de la ecuación.
# Matplotlib: para graficar función solución (theta(t))
# Numpy: Matemáticas y trabajo matricial.
import scipy, scipy.integrate
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
 
# Se define ecuación diferencial como sistema de dos
# ecuaciones de primer orden.
def dydt(t,y,L):
    theta,thetaPrima = y
    thetaPrimaPrima = -(g/L) * scipy.sin(theta)
    return [thetaPrima,thetaPrimaPrima]
 
# Definimos unas constantes a manera de ejemplo.
# g: aceleración de la gravedad.
# L: longitud del péndulo
# theta0: ángulo desde el que se suelta el péndulo.
# thetaPrima0: velocidad inicial con la que se suelta.
 
g = 9.8
L = 1.0
theta0 = scipy.pi/2 # Ángulo inicial: 90 grados.
thetaPrima0 = 0.0
 
solucionador = scipy.integrate.ode(dydt)
# 'dopri5' es para indicarle al solucionador que resuelva
# la ecuación por el método de Runge-Kutta 4,5
solucionador.set_integrator('dopri5')
# .set_f_params es la función para pasarle
# los argumentos adicionales a la ecuación diferencial.
# Para nuestro caso, la longitud del péndulo
solucionador.set_f_params(L)
# Le indicamos al solucionador el valor inicial del problema.
# theta0: ángulo desde el que se suelta el péndulo.
# thetaPrima0: velocidad inicial con la que se suelta.
solucionador.set_initial_value([theta0,thetaPrima0] , 0)
 
thetas = [] # Crear vector vacío de ordenadas.
dt = 0.01 # Pequeño intervalo delta de t.
while solucionador.successful() and solucionador.t < 10:
    solucionador.integrate(solucionador.t+dt)
    thetas.append(solucionador.y[0]*180/scipy.pi) # Llenar el vector de ordenadas con valor de theta
 
# Creamos vector de abscisas de tiempo.
t=[]
for i in range(len(thetas)):
    t.append(i*dt)
 
fig = plt.figure() # Creamos una nuevo gráfico con Matplotlib.
plt.plot(t,thetas) # Graficamos el vector de abscisas vs el vector de ordenadas.
plt.xlabel('tiempo (s)') # Título eje x
plt.ylabel('Amplitud theta(t)') # Título eje y
plt.show() # Mostramos el gráfico en pantalla.

IyBJbXBvcnRhbW9zIGxpYnJlcmlhcyBuZWNlc2FyaWFzCiMgU2NpUHksIFNjaVB5LGludGVncmF0ZTogbGlicmVyw61hcyBwYXJhIHNvbHVjacOzbiBwb3IgUnVuZ2UtS3V0dGEgNCw1IGRlIGxhIGVjdWFjacOzbi4KIyBNYXRwbG90bGliOiBwYXJhIGdyYWZpY2FyIGZ1bmNpw7NuIHNvbHVjacOzbiAodGhldGEodCkpCiMgTnVtcHk6IE1hdGVtw6F0aWNhcyB5IHRyYWJham8gbWF0cmljaWFsLgppbXBvcnQgc2NpcHksIHNjaXB5LmludGVncmF0ZQppbXBvcnQgbWF0cGxvdGxpYiBhcyBtcGwKaW1wb3J0IG1hdHBsb3RsaWIucHlwbG90IGFzIHBsdAppbXBvcnQgbnVtcHkgYXMgbnAKIAojIFNlIGRlZmluZSBlY3VhY2nDs24gZGlmZXJlbmNpYWwgY29tbyBzaXN0ZW1hIGRlIGRvcwojIGVjdWFjaW9uZXMgZGUgcHJpbWVyIG9yZGVuLgpkZWYgZHlkdCh0LHksTCk6CiAgICB0aGV0YSx0aGV0YVByaW1hID0geQogICAgdGhldGFQcmltYVByaW1hID0gLShnL0wpICogc2NpcHkuc2luKHRoZXRhKQogICAgcmV0dXJuIFt0aGV0YVByaW1hLHRoZXRhUHJpbWFQcmltYV0KIAojIERlZmluaW1vcyB1bmFzIGNvbnN0YW50ZXMgYSBtYW5lcmEgZGUgZWplbXBsby4KIyBnOiBhY2VsZXJhY2nDs24gZGUgbGEgZ3JhdmVkYWQuCiMgTDogbG9uZ2l0dWQgZGVsIHDDqW5kdWxvCiMgdGhldGEwOiDDoW5ndWxvIGRlc2RlIGVsIHF1ZSBzZSBzdWVsdGEgZWwgcMOpbmR1bG8uCiMgdGhldGFQcmltYTA6IHZlbG9jaWRhZCBpbmljaWFsIGNvbiBsYSBxdWUgc2Ugc3VlbHRhLgogIApnID0gOS44CkwgPSAxLjAKdGhldGEwID0gc2NpcHkucGkvMiAjIMOBbmd1bG8gaW5pY2lhbDogOTAgZ3JhZG9zLgp0aGV0YVByaW1hMCA9IDAuMAogIApzb2x1Y2lvbmFkb3IgPSBzY2lweS5pbnRlZ3JhdGUub2RlKGR5ZHQpCiMgJ2RvcHJpNScgZXMgcGFyYSBpbmRpY2FybGUgYWwgc29sdWNpb25hZG9yIHF1ZSByZXN1ZWx2YQojIGxhIGVjdWFjacOzbiBwb3IgZWwgbcOpdG9kbyBkZSBSdW5nZS1LdXR0YSA0LDUKc29sdWNpb25hZG9yLnNldF9pbnRlZ3JhdG9yKCdkb3ByaTUnKQojIC5zZXRfZl9wYXJhbXMgZXMgbGEgZnVuY2nDs24gcGFyYSBwYXNhcmxlCiMgbG9zIGFyZ3VtZW50b3MgYWRpY2lvbmFsZXMgYSBsYSBlY3VhY2nDs24gZGlmZXJlbmNpYWwuCiMgUGFyYSBudWVzdHJvIGNhc28sIGxhIGxvbmdpdHVkIGRlbCBww6luZHVsbwpzb2x1Y2lvbmFkb3Iuc2V0X2ZfcGFyYW1zKEwpCiMgTGUgaW5kaWNhbW9zIGFsIHNvbHVjaW9uYWRvciBlbCB2YWxvciBpbmljaWFsIGRlbCBwcm9ibGVtYS4KIyB0aGV0YTA6IMOhbmd1bG8gZGVzZGUgZWwgcXVlIHNlIHN1ZWx0YSBlbCBww6luZHVsby4KIyB0aGV0YVByaW1hMDogdmVsb2NpZGFkIGluaWNpYWwgY29uIGxhIHF1ZSBzZSBzdWVsdGEuCnNvbHVjaW9uYWRvci5zZXRfaW5pdGlhbF92YWx1ZShbdGhldGEwLHRoZXRhUHJpbWEwXSAsIDApCiAKdGhldGFzID0gW10gIyBDcmVhciB2ZWN0b3IgdmFjw61vIGRlIG9yZGVuYWRhcy4KZHQgPSAwLjAxICMgUGVxdWXDsW8gaW50ZXJ2YWxvIGRlbHRhIGRlIHQuCndoaWxlIHNvbHVjaW9uYWRvci5zdWNjZXNzZnVsKCkgYW5kIHNvbHVjaW9uYWRvci50IDwgMTA6CiAgICBzb2x1Y2lvbmFkb3IuaW50ZWdyYXRlKHNvbHVjaW9uYWRvci50K2R0KQogICAgdGhldGFzLmFwcGVuZChzb2x1Y2lvbmFkb3IueVswXSoxODAvc2NpcHkucGkpICMgTGxlbmFyIGVsIHZlY3RvciBkZSBvcmRlbmFkYXMgY29uIHZhbG9yIGRlIHRoZXRhCiAKIyBDcmVhbW9zIHZlY3RvciBkZSBhYnNjaXNhcyBkZSB0aWVtcG8uCnQ9W10KZm9yIGkgaW4gcmFuZ2UobGVuKHRoZXRhcykpOgogICAgdC5hcHBlbmQoaSpkdCkKIApmaWcgPSBwbHQuZmlndXJlKCkgIyBDcmVhbW9zIHVuYSBudWV2byBncsOhZmljbyBjb24gTWF0cGxvdGxpYi4KcGx0LnBsb3QodCx0aGV0YXMpICMgR3JhZmljYW1vcyBlbCB2ZWN0b3IgZGUgYWJzY2lzYXMgdnMgZWwgdmVjdG9yIGRlIG9yZGVuYWRhcy4KcGx0LnhsYWJlbCgndGllbXBvIChzKScpICMgVMOtdHVsbyBlamUgeApwbHQueWxhYmVsKCdBbXBsaXR1ZCB0aGV0YSh0KScpICMgVMOtdHVsbyBlamUgeQpwbHQuc2hvdygpICMgTW9zdHJhbW9zIGVsIGdyw6FmaWNvIGVuIHBhbnRhbGxhLg==