Péndulo esférico

anima.sh

#!/usr/bin/env bash
#
# Crea un vídeo de una serie de imágenes .png
# Juan Luis Cano Rodríguez <juanlu001@gmail.com>

ffmpeg -i fig%03d.png -b:v 512k animation.ogv

pendulo_esferico.py

# coding: utf-8
# 
# Péndulo esférico
# Juan Luis Cano Rodríguez <juanlu001@gmail.com>

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from numpy import sin, cos, pi

# du / dt = f(u, t)
# Ecuaciones de Hamilton para theta, psi, p_theta, p_psi
def f(u, t):
    return np.array([
        u[2],
        u[3] / (sin(u[0]) ** 2),
        u[3] ** 2 * cos(u[0]) / (sin(u[0]) ** 3) + 10 * sin(u[0]),
        0
    ])

# Vector de condiciones iniciales
u0 = np.array([pi / 3, 0, 0, 0.3])

N = 1000
t = np.linspace(0, 10, N + 1)

sol = odeint(f, u0, t)

theta = sol[:, 0]
psi = sol[:, 1]

# Representación gráfica
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.autoscale(False)
ax.set_xlim(-1.0, 1.0)
ax.set_ylim(-1.0, 1.0)
ax.set_zbound(-1.0, 1.0)

# Se producen N / step imágenes para montar el vídeo
step = 10
for i in xrange(N / step + 1):
    ax.plot(sin(theta[:i * step]) * cos(psi[:i * step]), sin(theta[:i * step]) * sin(psi[:i * step]), cos(theta[:i * step]), color='#204a87')
    plt.savefig('figs/fig{:03}.png'.format(i))
    ax.lines.pop()