Indice | Previo: Mecanica.Cinematica.SolucionEdM | Siguiente: Mecanica.Dinamica.Cantidades
5.2.6. Integración numérica de la e.d.m.¶
En aquellos sistemas dinámicos en los que resolver la e.d.m. o encontrar todas las constantes de movimiento (resolver por cuadratura el sistema), sea imposible matemáticamente o simplemente muy difícil, es posible buscar una solución aproximada usando los métodos numéricos que estudiamos al final de la Seccion Ecuaciones diferenciales.
Como vimos allí, para hacerlo, es necesario primero escribir la e.d.m. como el conjunto de 6 ecuaciones diferenciales de primer orden con la forma general (Ec. ecuaciones_reducidas):
\begin{equation} \label{eq:edm_ecuaciones_reducidas} \{\dot{Y_i} = f_i(t,\{Y_k\})\}_{6} \end{equation}Donde $Y_i$ ($i=0,1,2,\ldots,5$) son las denominadas funciones auxiliares que reemplazan aquí a las cantidades claves del sistema dinámico (las componentes de $\vec r$ y $\vec v$). En física, llamaremos a las Ecs. (edm_ecuaciones_reducidas), ecuaciones de movimiento reducidas del sistema o e.d.m.r.
Si comparamos la forma general de las e.d.m.r. en la Ec. (edm_ecuaciones_reducidas) con las ecuaciones de primer orden (edm_primer_orden_escalar), podemos hacer la siguiente identificación para las variables auxiliares $Y_i$:
\begin{equation} \label{eq:variables_auxiliares_edm} \begin{array}{ccc} Y_0 = x, & Y_1 = y, & Y_2 = z\\ Y_3 = v_x, & Y_4 = v_y, & Y_5 = v_z\\ \end{array} \end{equation}Por otro lado, las funciones $f_i$ serán:
$$ \begin{array}{ccc} f_0(t,\{Y_k\}) = v_x = Y_3, & f_1(t,\{Y_k\}) = v_y = Y_4, & f_2(t,\{Y_k\}) = v_z = Y_5\\ f_3(t,\{Y_k\}) = a_x, & f_4(t,\{Y_k\}) = a_y, & f_5(t,\{Y_k\}) = a_z\\ \end{array} $$Con esta identificación, una forma compacta de escribir las e.d.m.r., muy útil a la hora de preparar algoritmos, es:
\begin{equation} \label{eq:edm_reducidas_particula} \dot Y_i = \left\{ \begin{array}{ccc} Y_{3+i} & {\rm ,} & 0 \leq i < 3 \\ a_{i-3} & {\rm ,} & 3\leq i < 6 \end{array} \right. \end{equation}Donde hemos introducido la notación $\{a_0,a_1,a_2\}\equiv\{a_x,a_y,a_z\}$.
Integración numérica de las e.d.m. del ejemplo 1¶
El sistema dinámico del ejemplo 1 introducido en la Sección Integración por cuadraturas se caracteriza por tener una aceleración del tipo $\vec a:(j_0t,0,0)$. En términos de la parametrización de la e.d.m.r. en la Ec. (edm_reducidas_particula), esto significa que $\{a_0,a_1,a_2\}=(j_0t,0,0)$.
El sistema de ecuaciones diferenciales que describe el sistema puede implementarse en Python usando la rutina:
def edm_ejemplo1(Y,t,j0=1):
dYdt=[0,0,0,0,0,0]
dYdt[:3]=Y[3:]
dYdt[3:]=[j0*t,0,0]
return dYdt
Aunque al final de la Seccion Ecuaciones diferenciales nos habíamos familiarizado con este tipo de rutinas, el diseño de está en particular merece algunos comentarios:
Como sabemos, el propósito de esta rutina es calcular la lista de los valores de las funciones $f_i$ (lado derecho de las Ecs. edm_ecuaciones_reducidas) que son iguales a las primeras derivada en el tiempo de las variables auxiliares $\dot{Y}_i$. En la rutina, para hacer más explícito su significado, hemos decidido llamar a esta lista
dYdten lugar def. Los nombres de las variables no afectan la funcionalidad de las rutinas, pero pueden hacerla más legible y modificable.Para asignar los valores de la lista
dYdthemos aprovechado el poder dePythonpara sacar trozos, porciones o tajadas (slices en inglés) de listas y arreglos. Así el trozodYdt[:3]corresponde a las primeras tres componentesdYdt[0],dYdt[1],dYdt[2](nótese que, por empezar en 0, este trozo no incluye la componentedYdt[3].) Por otra parte el trozoY[3:]de esta lista, corresponde a las componentesY[3],Y[4],Y[5]. Por tanto, la igualdaddYdt[:3]=dYdt[3:]equivale a escribir explícitamentedYdt[0]=Y[3],dYdt[1]=Y[4],dYdt[2]=Y[5]que justamente implementa la parte $\dot{Y}_i=Y_{i+3}$ de la Ec. (edm_reducidas_particula)
Una vez escrita la rutina, la solución se obtiene siguiendo los algoritmos introducidos al final de la Seccion Ecuaciones diferenciales:
#Valor del tiron
j0=0.5
#Condiciones iniciales
from numpy import array
Yos=array([1,0,0,-3,0,0])
#Tiempos para obtener la solución
from numpy import linspace
ts=linspace(0.0,5.0,20)
#Solución con odeint
from scipy.integrate import odeint
Ys=odeint(edm_ejemplo1,Yos,ts,args=(j0,))
Podemos, finalmente, visualizar la solución a la e.d.m.r. haciendo un gráfico de la coordenada $x$ (columna Ys[:,0] de la matriz de solución) como función del tiempo ts:
%matplotlib inline
import matplotlib.pyplot as plt
plt.figure();
plt.plot(ts,Ys[:,0],'ko');
plt.xlabel("$t$");
plt.ylabel("$x$");
plt.show();
¿Cómo saber si la solución obtenida con odeint y mostrada en la Figura (code:ejemplo1_grafico_solucion) es la correcta?
Existen dos manera de comprobarlo. La primera es verificar que las posiciones y velocidades obtenidas satisfagan las constantes de movimiento que escribimos en las Ecs. (constantes_ejemplo1).
Así por ejemplo, podemos verificar que el valor de $I_x$ y $S_x$ sean efectivamente constantes:
#Extraemos los valores de x y dxdt de la solución
xs=Ys[:,0]
xdots=Ys[:,3]
#Fórmula de la constante C_Ix
C_Ixs=xdots-0.5*j0*ts**2
#Fórmula de la constante C_Sx
C_Sxs=xs-(1./6)*j0*ts**3-C_Ixs*ts
Comprobamos así que las fórmulas $C_{Ix}=\dot{x}-j_0t^2/2$ y $C_{Sx}=x-j_0t^3/6-I_x t$, tienen el mismo valor para todos los tiempos en los que integramos la e.d.m.r., es decir, son, por definición, constantes de movimiento. La solución numérica, por tanto,satisface nuestras expectativas matemáticas.
La segunda manera de verificar que nuestra solución numérica coincida con la analítica es compararla con la solución explícita escrita en las Ecs. (solucion_ejemplo1):
$$ x(t)=\frac{1}{6}j_0 t^3+I_x t+S_x $$Aquí, los valores de $I_x$ y $S_x$ pueden obtenerse de las condiciones iniciales.
Una comparación gráfica entre ambas soluciones se consigue con este algoritmo:
#Valor de las constantes de movimiento
Ix=Yos[3]-0.5*j0*ts[0]**2
Sx=Yos[0]-(1./6)*j0*ts[0]**3+Ix*ts[0]
#Lista más completa de valores del tiempo
from numpy import linspace
tas=linspace(ts[0],ts[-1],100)
#Solución analítica
xs=(1./6)*j0*tas**3+Ix*tas+Sx
#Gráfico
plt.figure();
plt.plot(ts,Ys[:,0],'ko',label="Solución numérica");
plt.plot(tas,xs,'r-',label="Solución analítica");
plt.xlabel("$t$")
plt.ylabel("$x$")
plt.legend();
plt.show();
La coincidencia entre la solución analítica y la solución numérica mostrada en la Figura (code:compara_numerica_analitica) es casi perfecta.
Integración numérica de las e.d.m. del ejemplo 2¶
Usando las mismas herramientas y algoritmos análogos a los usados antes, podemos ahora resolver el ejemplo 2 de la Sección Integracion de la e.d.m..
De nuevo, las ecuaciones reducidas del sistema serán, como en el ejemplo 1, las mismas de la Ec. (edm_reducidas_particula), pero ahora $\{a_i\}=\{-\omega Y_0,0,0\}$ (nótese que hemos reemplazado $x$ por la variable auxiliar $Y_0$ de acuerdo a las reglas en Ec. variables_auxiliares_edm).
La rutina que implementa las edm en este caso será:
def edm_ejemplo2(Y,t,omega):
dYdt=[0,0,0,0,0,0]
dYdt[:3]=Y[3:]
dYdt[3:]=[-omega*Y[0],0,0]
return dYdt
La solución al sistema, una comprobación de que la constante $C_{Ix}$ en la Ec. (Ix_ejemplo2) es en realidad una constante, y una gráfica de la posición como función del tiempo, se muestra en el siguiente algoritmo:
#Propiedades del sistema
omega=2.5
#Condiciones iniciales
from numpy import array
Yos=array([1,0,0,0,0,0])
#Tiempos
from numpy import linspace
ts=linspace(0,10.0,100)
#Solución
Ys=odeint(edm_ejemplo2,Yos,ts,args=(omega,))
#Constante de movimiento
xs=Ys[:,0]
xdots=Ys[:,3]
C_Ixs=0.5*xdots**2+0.5*omega*xs**2
#Gráfico
fig=plt.figure();
plt.plot(ts,Ys[:,0]);
plt.xlabel("$t$");
plt.ylabel("$x$");
plt.show();
Podemos comprobar al examinar la Figura (code:solucion_numerica_ejemplo2) nuestra intuición inicial de que la dinámica del sistema correspondía a la de un movimiento oscilatorio.
Indice | Previo: Mecanica.Cinematica.SolucionEdM | Siguiente: Mecanica.Dinamica.Cantidades