Consideramos el sistema diferencial con condiciones iniciales siguiente: $$ \begin{cases} \dot{y_1} = y_2, \\ \dot{y_2} = -y_1, \\ y_1(0) = 1,\, y_2(0) = 0. \end{cases} $$
Para resolverlo numéricamente, utilizamos la función scipy.integrate.odeint que recibe como argumentos una función que representa la ecuación diferencial $\dfrac{\partial y}{\partial t}$ con $y = (y_1, y_2)$, el punto inicial $y(0) = (y_1(0), y_2(0))$ y un valor para el tiempo $t$.
A continuación, resolvemos el sistema diferencial para $t \in (0, 2 \pi)$ y graficamos la solución.
In [1]:
%pylab inline
In [2]:
import numpy as np
import matplotlib.pylab as plt
from scipy.integrate import odeint
def dydt(y, t):
"""Función que representa el sistema diferencial dy/dt = f(y, t)."""
return [y[1], -y[0]]
In [3]:
# Valores del tiempo.
t = np.linspace(0, 2*np.pi)
# Punto inicial.
y0 = [1, 0]
In [4]:
# Resolvemos el sistema diferencial para los valores de t.
y = odeint(dydt, y0, t)
# Graficamos la solución.
plt.plot(y[:,0], y[:,1])
plt.show()
Ahora consideramos el siguiente sistema diferencial $$ \begin{cases} \dot{y_1} = y_2, \\ \dot{y_2} = -y_1, \\ \dot{y_3} = 1, \\ y_1(0) = 1,\, y_2(0) = 0,\, y_3(0) = 0, \end{cases} $$
para $t \in (0, 4\pi)$.
De la misma manera que en la sección anterior, representamos el sistema diferencial como una función y resolvemos utilizando la función odeint para el punto inicial y valores de $t$.
In [5]:
def dydt(y, t):
return [-y[1], y[0], 1]
t = np.linspace(0, 4*np.pi)
y0 = [1, 0, 0]
y = odeint(dydt, y0, t)
Ahora, creamos una figura con 4 gráficas, 3 gráficas de las proyecciones en cada plano y la última la gráfica de la curva en 3D.
In [6]:
# Este import activa el soporte para graficación en 3D.
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
# Para esta gráfica especificamos la opcíon para graficar en 3D.
ax4 = fig.add_subplot(224, projection='3d')
# Graficamos la proyección en el primer plano.
ax1.set_title(u'Proyección X-Y')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.plot(y[:,0], y[:,1])
# Graficamos la proyección en el segundo plano.
ax2.set_title(u'Proyección Y-Z')
ax2.set_xlabel('y')
ax2.set_ylabel('z')
ax2.plot(y[:,1], y[:,2])
# Graficamos la proyección en el tercer plano.
ax3.set_title(u'Proyección X-Z')
ax3.set_xlabel('x')
ax3.set_ylabel('z')
ax3.plot(y[:,0], y[:,2])
# Graficamos la curva directamente en 3D.
ax4.set_title(u'Gráfica 3D')
ax4.set_xlabel('x')
ax4.set_ylabel('y')
ax4.set_zlabel('z')
ax4.plot(y[:,0], y[:,1], y[:,2])
# Quitamos la numeración para mejorar la presentación.
ax4.set_xticklabels([])
ax4.set_yticklabels([])
ax4.set_zticklabels([])
# Ajustamos el espaciado de las gráficas.
fig.tight_layout()
plt.show()
Para más información, ver los siguientes enlaces: