Una pelota de tenis mojada que rota

Hace poco me encontré esta imagen de una pelota de tennis mojada rotando sobre su propio eje (fuente):


La imagen me pareció interesante, más que por el sentido estético, porque como físico sé que las gotas deben salir por la tangente a cada punto de la superficie de la pelota y describir líneas rectas. Como sea, en la imagen los chorros de agua prácticamente parecen salen por la normal a cada punto y en realidad forman espirales.

¿Qué curva describen realmente las gotas de agua y cómo es compatible con la noción física usual?
Responderlo es relativamente sencillo utilizando únicamente mecánica clásica (matemáticamente no más que cálculo vectorial).

Para empezar podemos reducir la situación a una dos dimensional, y con la pelota siendo un disco (del cual sólo importa la superficie) y las gotas partículas puntuales. Por simplicidad, consideremos una pelota de radio unitario y que rota con una velocidad angular constante $\omega$, de modo que la posición $\vec{r}$ de una gota descrita en la superficie de la pelota, que idealizamos como circular, a cada tiempo $t\geq0$, está dada por
\begin{equation}\vec{r}(t)\equiv(\cos\omega{t},\sin\omega{t})\end{equation} y entonces, su velocidad tangencial es
\begin{equation}\vec{v}(t)\equiv\frac{d\vec{r}}{dt}=\omega(-\sin\omega{t},\cos\omega{t})\end{equation} Así entonces, una gota que a un momento $t_0\in[0,t]$ sale de la superficie de la pelota describe una recta dada por
\begin{align}\vec{\mathcal{D}}_{t_0}(t)&=\vec{r}(t_0)+\vec{v}(t_0)(t-t_0)\nonumber\\
&=(\cos\omega{t_0}-\omega(t-t_0)\sin\omega{t_0},\sin\omega{t_0}+\omega(t-t_0)\cos\omega{t_0})\end{align}

Rotación sentido antihorario
Vale, el problema realmente está resuelto, sólo queda considerar una familia de gotas e intentar visualizar y comprobar si podemos reproducir el fenómeno que se observa en la imagen. Consideremos primero rotaciones en sentido antihorario, i.e. $\omega>0$. Por simplicidad hagamos $\omega=1$.

El caso más simple es el de $t_0=0$,
\begin{align}\vec{\mathcal{D}}_0(t)&=(1,t)\end{align}

que es bastante aburrido ;-) Consideremos entonces una colección de gotas, donde la $n$-ésima gota está dada por
\begin{equation}\vec{r}_n(t)=\left(\cos(t+\theta_n),\sin(t+\theta_n)\right)\end{equation} cuya velocidad tangencial es
\begin{equation}\vec{v}_n(t)\equiv\frac{d\vec{r}_n}{dt}=\left(-\sin(t+\theta_n),\cos(t+\theta_n)\right)\end{equation} y de manera análoga, ahora
\begin{align}\vec{\mathcal{D}}_{t_0,\theta_n}(t)&=\vec{r}_n(t_0)+\vec{v}_n(t_0)(t-t_0)\end{align} y para $t_0=0$,
\begin{align}\vec{\mathcal{D}}_{0,\theta_n}(t)&=(\cos\theta_n-t\sin\theta_n,\sin\theta_n+t\cos\theta_n)\end{align}

Aquí he tomado los $\theta_n$'s igualmente espaciados por facilidad en Mathematica, pero en general las gotas podrían desprenderse en cualquier punto.

Viene entonces el caso que nos interesa, que es que se desprendan gotas casi continuamente de la pelota, es decir, considerar diversos valores $t_0$ en el intervalo $[0,t]$ muy ligeramente espaciados. Las ecuaciones para la gota en $\theta_n$ que se desprende al tiempo $t_0$ es
\begin{align}\vec{\mathcal{D}}_{t_0,\theta_n}(t)&=\left(\cos(t_0+\theta_n)-(t-t_0)\sin(t_0+\theta_n),\sin(t_0+\theta_n)+(t-t_0)\cos(t_0+\theta_n)\right)\label{dranticlock}\end{align} Graficando y visualizando en el tiempo varias de ellas se obtiene algo así


y prácticamente hemos reproducido el fenómeno de la fotografía. Vemos entonces que la pregunta que realmente queríamos responder es: ¿qué curva describen los 'chorros' de agua? Y eso es sencillo, ¡sólo hay que tomarle una foto a la ecuación (\ref{dranticlock})! Una forma de hacer esto es tomar un $t>0$ dado y graficar en $0\leq{t}_0\leq{t}$ para los diversos $\theta_n$, a modo de 'regresar la película' para cada $t_0$ ;-) El caso más sencillo es tomar $t=2\pi$; considera además por ahora $\theta_n=0$,
\begin{equation}\mathcal{C}_{\theta_n=0}(t_0)\equiv\vec{\mathcal{D}}_{t_0,0}(2\pi)=\left(\cos{t}_0-(2\pi-t_0)\sin{t}_0,\sin{t}_0+(2\pi-t_0)\cos{t}_0\right)\end{equation} que podemos escribir como
\begin{equation}\vec{\mathcal{C}}_{0}(t_0)=\left(\cos(2\pi-t_0)-(2\pi-t_0)\sin(t_0-2\pi),\sin(t_0-2\pi)+(2\pi-t_0)\cos(2\pi-t_0)\right)\end{equation} ya que $\cos(2\pi-x)=\cos{x}$ y $\sin(x-2\pi)=\sin{x}$, entonces haciendo $\tau\equiv{2\pi-t_0}$,
\begin{align}\vec{\mathcal{C}}_{0}(\tau)&=\left(\cos\tau-\tau\sin(-\tau),\sin(-\tau)+\tau\cos\tau\right)\nonumber\\
&=\left(\cos\tau+\tau\sin\tau,-\sin\tau+\tau\cos\tau\right)\end{align} que es la ecuación paramétrica de la curva llamada evolvente. Si consideramos el caso general para cualquier $\theta_n$, llegamos a
\begin{align}\vec{\mathcal{C}}_{\theta_n}(\tau)&=\left(\cos(\tau-\theta_n)+\tau\sin(\tau-\theta_n),-\sin(\tau-\theta_n)+\tau\cos(\tau-\theta_n)\right)\end{align}

y en efecto puede comprobarse que los chorros, i.e. la evolvente del círculo, es normal a cada punto de la superficie (e.g con $\vec{v}_n(0)\cdot\vec{\mathcal{C}}_{\theta_n}(0)=0$).

Rotación sentido horario
La única razón de haber tomado la rotación en sentido antihorario (además de la sospecha de que así se reproduciría la foto), es que no quería estar rastreando un signo negativo en las ecuaciones. Pero Mathematica no tiene ese problema, así que sólo hay que tomar $\omega<0$ y ya está.

Finalmente, cuando buscaba más imágenes e información sobre los chorros de agua, seguido me encontré con que se asume que los chorros describen o una espiral logarítmica o una espiral de Fibonacci, sobre todo -supongo- porque estas espirales son populares con otros fenómenos de la naturaleza. Como sea, acá estoy sobre-simplificando la cuestión y no descarto la posibilidad de que así sea (aunque no he encontrado ningún recurso que lo pruebe): el tratamiento de los chorros de agua de la pelota seguramente se puede complicar tanto como se quiera, permitiendo velocidades angulares grandes, incluyendo efectos como la interacción entre gotas, la fricción del aire, presencia de un campo gravitacional, etc...

Acá dejo el código que utilicé para hacer los GIF's en Mathematica:

Del Último Teorema de Fermat a los Espacios de Calabi-Yau

El Universo Elegante fue casi seguramente el primer libro de divulgación sobre Cuerdas que leí (quizá sólo luego de Hiperespacio de Michio Kaku) en la biblioteca de la UAM-A y seguramente uno de los que me inclinaron más hacia la ciencia que hacia la ingeniería. Como sea, recuerdo bien que en algún punto del libro se incluye una ilustración de una red de variedades Calabi-Yau, en la que se explica, están compactificadas las dimensiones espaciales extra.

Fuente: http://members.wolfram.com/jeffb/visualization/stringtheory

De manera fortuita (en realidad por el video al final ;-), hace poco encontré también esta charla de Andrew J. Hanson:


en la que se muestra también este vídeo.

Aunque mi expertise con los Calabi-Yau (el nombre debido a Eugenio Calabi y Shing-Tung Yau) es básicamente la misma que cuando leí el libro de Brian Greene, o sea, nula, al menos ya tengo acceso a las ideas básicas.

www.amazon.com/The-Shape-Inner-Space-Dimensions/dp/0465028373

Hasta donde sé, hay varias maneras de definir una variedad de Calabi-Yau; de cualquier modo, algunas de sus características son que es una variedad compleja, compacta, Kähler y de curvatura de Ricci nula.

Esto es,

• Un $d$-fold $\mc{M}$ de dimensión compleja $d$ (o $\dim_\mathbb{R}(\mc{M})=2d$) de coordenadas locales $z^i$, $z^\bar{\imath}$ con $i,\bar{\imath}=1,\ldots,d$
• Que puede cubrirse con una cantidad finita de parches coordenados
• Que admite una ${(1,1)}$-forma $\omega$ cerrada, i.e. tal que $d\omega=0$, que se relaciona con una métrica $g$ Hermitiana en $\mc{M}$, i.e. real y tal que $g_{ij}=g_{\bar{\imath}\bar{\jmath}}=0$, que localmente puede escribirse en términos de derivadas de una función real $K=K(z,\bar{z})$ llamada potencial de Kähler, $g_{i\bar{\jmath}}=\frac{\p{K}}{\p{z}^i\p{z}^\bar{\jmath}}$. En general esto y los dos puntos anteriores suelen decirse de una variedad Kähler.
• Cuyas componentes del tensor de Ricci $R_{i\bar{\jmath}}$ se anulan, $R_{i\bar{\jmath}}=0$, además, evidentemente (por el punto anterior) de $R_{ij}$, $R_{\bar{\imath}\bar{\jmath}}$

Éstas, según el autor, pueden tratarse como propiedades o como una definición (espero no hacer enojar a algún purista ;-) y es relevante mencionar que el último punto de hecho suele mencionarse más estrictamente como una consecuencia del llamado Teorema de Yau, proveniente de la llamada antes Conjetura de Calabi (y nacieron los Calabi-Yau ;-), que tiene que ver con que se anulen ciertos invariantes topológicos llamados primera clase de Chern, $c_1(\mc{M})$. El Teorema de Yau dice que dado $\mc{M}$ Kähler con $c_1(\mc{M})=0$ siempre existe una métrica tal que $R_{i\bar{\jmath}}=0$, por esto seguido sólo se dice que un $\mc{M}$ Kähler con $c_1(\mc{M})=0$ es un Calabi-Yau.

Como sea, mi intención es sólo dar una idea de lo que es una variedad de Calabi-Yau a quien -como yo- ya tiene estos conceptos a la mano. Para ahondar en el tema se puede consultar

• String Theory on Calabi-Yau Manifolds de Brian Greene: http://arxiv.org/abs/hep-th/9702155
• A mini-course on topological strings de Marcel Vonk: http://arxiv.org/abs/hep-th/0504147

En las notas de Vonk puede leerse que la última condición es interesante porque las variedades con tensor de Ricci nulo son soluciones de vacío a las ecuaciones de Einstein (con constante cosmológica nula), además de que se asegura que al compactificar 6 dimensiones extra en las variedades de Calabi-Yau, se preservan las supersimetrías de la teoría 4-dimensional, lo que (sea lo que sea) es relevante por razones técnicas y fenomenológicas. El punto es que esto lleva a que en cuerdas se empleen 3-folds Calabi-Yau (i.e. 6 dimensiones reales).

En las notas de Greene se muestra (2.10 y 2.11) que la hipersuperficie (quíntica o de quinto grado) en el espacio proyectivo complejo $\mathbb{C}P^4$ dada por
\begin{equation}z_1^5+z_2^5+z_3^5+z_4^5+z_5^5=0\end{equation} es un 3-fold de Calabi-Yau, mostrando únicamente que es Kähler con primera clase de Chern nula. Finalmente de esta forma es como se relaciona todo con el Último Teorema de Fermat y todo lo que se menciona en el vídeo de Hanson.

Andrew J. Hanson, Indiana University. www.cs.indiana.edu/~hanson

Reproducir los gráficos del corte 2D del Calabi-Yau de Hanson es más complicado de lo que parece. En el artículo

• A.J. Hanson. A construction for computer visualization of certain complex curves. Notices of the Amer. Math. Soc., 41(9):1156-1163, November/December 1994.

se menciona el procedimiento general. De cualquier modo uno siempre puede ocupar los Wolfram Demonstrations disponibles y destripar el código fuente si lo hay. Dos .cdf que pude hallar son

• Calabi-Yau Space: demonstrations.wolfram.com/CalabiYauSpace
• Solutions of Fermat's Equation: demonstrations.wolfram.com/SolutionsOfFermatsEquation

siendo el segundo el que tiene el código fuente disponible.

Finalmente, ya entrados en esto de los Calabi-Yau:

Geodésicas en AdS

Hace poco publiqué esta respuesta en Physics SE, misma que no fue elegida por el interrogador, que eligió en su lugar otra innecesariamente embrollada matemáticamente, como si el interrogador hubiese estado buscando cátedra en geometría diferencial cuando ignoraba siquiera lo que es un parámetro afín.

De cualquier modo, yo también me lié un poco con mi respuesta. En la parte de las geodésicas tipo luz básicamente copié el cálculo que hice en mi trabajo de servicio social, pero en el caso de las geodésicas tipo tiempo y tipo espacio me encontré con algunas dificultades (en realidad no pensé en la necesidad o plausibilidad de calcularlas explícitamente en el servicio social) y no pude ahondar más porque estaba en pleno periodo trimestral en la Universidad.

El objetivo de esta entrada entonces es básicamente reproducir analíticamente las geodésicas en el esquema del diagrama de Penrose de AdS que muestro acá:

y graficar con algún software numérico, e.g. Mathematica, para ello voy a emplear la métrica (con variables angulares fijas)
\begin{equation}ds^2=-(1+r^2)dt^2+(1+r^2)dr^2\end{equation} con $r=\tan\mc{R}\in\mathbb{R}^+$, y para la cubierta universal, $t\in\mathbb{R}$.

Para las geodésicas de tipo tiempo, tomando $\lambda$ como el tiempo propio, ya que (tomando signatura +---) $\dot{s}^2=-1$, se tiene que
\begin{equation}-(1+r^2)\dot{t}^2+(1+r^2)\dot{r}^2=-1\end{equation} Ahora bien, $\p_t$ (i.e. el vector de componentes $V^t=1$, $V^i=0$) es un vector de Killing global, entonces
\begin{equation}V_t\dot{t}=g_{tt}\dot{t}=-(1+r^2)\,\dot{t}\equiv-E\label{da}\end{equation} es una constante de movimiento (la energía, siendo que proviene de una simetría $t$-traslacional); donde el signo negativo es simplemente tal que $\dot{t}>0$. De aquí entonces
\begin{equation}\dot{r}^2+(1+r^2)-E^2=0\end{equation} es decir,
\begin{equation}\frac{dr}{\pm\sqrt{(E^2-1)-r^2}}=d\lambda\end{equation} cuya solución más general, con la condición (arbitraria) $r(0)=0$, es
\begin{equation}r(\lambda)=\pm\sqrt{E^2-1}\tan\lambda\,|\cos\lambda|\end{equation} que de cualquier modo, restringiendo $r(\lambda)>0$, se tiene que reducir inevitablemente el dominio de $\lambda$ a intervalos de $\pi/2$ respecto a $\lambda=0$, es decir,
\begin{equation}r(\lambda)=\begin{cases}\sqrt{E^2-1}\tan\lambda\,|\cos\lambda|,&\text{si}&\lambda\in[n\pi,\frac{(2n+1)\pi}{2}]\\[0.1in]-\sqrt{E^2-1}\tan\lambda\,|\cos\lambda|,&\text{si}&\lambda\in[\frac{(2n+1)\pi}{2},(n+1)\pi]\end{cases},\,n\in\mathbb{Z}\end{equation} esto por supuesto de algún modo está relacionado con que AdS no está extendido a-priori, o que está enrollado en la coordenada temporal. Tomando cualquier caso, de la ec. (\ref{da}) con $t(0)=0$, se sigue que
\begin{equation}t(\lambda)=\arctan(E\,\tan\lambda)\end{equation} y así entonces escribiendo $\lambda=\lambda(t)$, se sigue finalmente que (asumiendo $E>0$; si se asume lo contrario los intervalos se invierten correspondientemente),
\begin{equation}r(t)=\begin{cases}\sqrt{E^2-1}\frac{\tan{t}}{E\sqrt{1+\frac{\tan^2t}{E^2}}},&\text{si}&t\in[n\pi,\frac{(2n+1)\pi}{2}]\\[0.1in]-\sqrt{E^2-1}\frac{\tan{t}}{E\sqrt{1+\frac{\tan^2t}{E^2}}},&\text{si}&t\in[\frac{(2n+1)\pi}{2},(n+1)\pi]\end{cases},\,n\in\mathbb{Z}\end{equation} y simplemente hay que usar $\mc{R}=\arctan{r}$. De aquí uno puede procurar simplificar cosas e.g. tomar intervalos de $\pi$ respecto a $t=0$ haciendo $\tan{t}\to|\tan{t}|$ (esto equivale a lo que escribí en la respuesta de Physics SE para $r(\lambda)$), lo que puede justificarse siendo removible la discontinuidad entre cada $\pi/2$ respecto a 0; de algún modo esto recupera la $2\pi$-periodicidad de AdS en la cubierta universal. Esta cuestión de los intervalos vuelve un poco delicado el tratamiento de este problema de obtener analíticamente las expresiones para las geodésicas; el caso de las geodésicas de tipo espacio es enteramente análogo tomando $\dot{s}^2=1$ en el tiempo propio $\lambda$, pero nuevamente uno debe ser cuidadoso al considerar en dónde están definidas las soluciones físicamente relevantes.

Graficando las soluciones entonces para diversos valores de $E$, se obtiene (extendiendo a la cubierta universal)
como se esperaba, si se quiere rotar el gráfico, simplemente hay que graficar $t=t(\mc{R})$.

Finalmente una cuestión relevante es que, independientemente del signo que tenga, en el caso de geodésicas tipo tiempo, la cantidad $E$ debe satisfacer $|E|\geq1$. De esto no he encontrado nada al respecto en otros lados y no lo he podido discutir ampliamente con mi asesor, sin embargo no es descabellado que estas geodésicas tengan una restricción en la energía pues implican movimiento masivo y la cantidad 1 está asociada a un cambio de longitud de arco en el tiempo propio, de modo que para que una partícula se mueva en una de estas geodésicas es coherente que su energía deba exceder o igualar en proporción esta razón de cambio (en unidades correspondientes), si la iguala, simplemente no se mueve (i.e. se queda en el polo $\mc{R}=0$). También debe notarse que cuando $E\to\infty$, las geodésicas se tornan de tipo luz, como también se esperaba.

La geometría clásica del espacio de (Anti)-de Sitter

He pasado poco más de seis meses entrándole a la geometría (clásica) de los espacios (A)dS en mi proyecto de servicio social en la licenciatura. El producto son estas notas de casi 200 páginas que cubren los temas de Relatividad Especial, Relatividad General y los espacios de (Anti)-de Sitter. La razón de haberlo hecho así es que en realidad tuve que aprender todos estos temas a lo largo del servicio. Las notas me parecen bien digeridas para cualquier estudiante de licenciatura interesado en temas de gravitación y no creo que haya ningún problema para cualquiera que las quisiera leer. El trabajo presuntamente está libre de errores conceptuales mayores y a lo más podrían existir errores estéticos o de dedo, mismos que procuraré señalar aquí al encontrarlos.

Resumen
Los espacios de de Sitter y de Anti-de Sitter son básicamente espaciotiempos de curvatura constante con ciertas propiedades peculiares que resultan físicamente relevantes. Estos espacios pueden verse como una solución de las ecuaciones de campo de la Relatividad General cuando se considera una constante cosmológica ya sea positiva o negativa. Si bien el espacio ${AdS}$ juega un papel limitado en la Relatividad General, contrario al espacio ${dS}$, éste puede resultar de gran relevancia en otras teorías que incorporan gravedad como lo es teoría de cuerdas dentro del marco de la dualidad AdS/CFT. Aunque este trabajo no es exhaustivo, se busca introducir primeramente los elementos físicos y matemáticos básicos de la Relatividad Especial y General, dado que estos temas actualmente pertenecen únicamente a cursos optativos en la UAM-I. En general se cubren tres secciones básicas: 1) Relatividad Especial, 2) Relatividad General, y 3) Geometría y el espacio $(A)dS$.

En la sección 1 se introduce la notación tensorial y algunos conceptos básicos sobre tensores, luego se hace un resumen de los conceptos básicos de la Relatividad Especial y como complemento se tratan brevemente la electrodinámica relativista y las ecuaciones cuántico-relativistas. Finalmente se muestran diversos problemas resueltos, cuya formulación ha sido traducida de las diversas fuentes bibliográficas citadas.

En la sección 2 se presentan de manera más sólida los conceptos acerca de tensores y se introduce el tema de variedades diferenciables y curvatura a manera de formular la teoría de manera tensorial. Finalmente se obtienen las ecuaciones de campo por formulación Lagrangiana y asimismo se muestran algunos problemas resueltos, con el principal objetivo de complementar el contenido de la sección.

En la sección 3 se tratan finalmente los espacios clásicos de ${(A)dS_n}$ introduciendo diversos conceptos geométricos que aquí se dirigen principalmente a los campos de Killing y los espacios máximamente simétricos. A partir de aquí se estudian diversos aspectos de ${(A)dS_n}$ en varios sistemas coordenados y se ilustran los correspondientes diagramas de Penrose a modo de estudiar la estructura causal de cada espacio. Ambos espacios en los distintos sistemas coordenados se logran a través de diversas parametrizaciónes de hipersuperficies que luego se encajan dentro de un espacio coordenado de dimensión mayor. El tratamiento en ambos casos se hace en un número de dimensiones arbitrarias, de modo que el lector pueda simplemente considerar cualquier caso de interés. Finalmente se tratan los espacios ${(A)dS_n}$ en el contexto de la Relatividad General, considerando el papel que juegan desde un punto de vista cosmológico y se revisa la solución de $AdS$-Schwarzschild.

Enlace al documento en la Colección de Tesis Electrónicas de la UAM-I:
http://tesiuami.izt.uam.mx/uam/aspuam/presentatesis.php?recno=16576&docs=UAMI16576.pdf

Toros Invariantes KAM y Huecos Resonantes

Hace unos días presenté junto con una compañera este tema, al que sólo agregaría en el título, "en el Sistema Solar", ya que expusimos en específico el caso del cinturón de asteroides y de los anillos de Saturno. El tema es de lo más divertido; es como si uno adquiriese algo así como un elemento de entendimiento riquísimo acerca de la belleza del Sistema Solar, y en general, del universo.

Acá dejo pues, la presentación: Toros Invariantes KAM y Huecos Resonantes.
** Actualización: En mayo de 2014 presentamos este tema en el Seminario de Alumnos de la UAM-I; tiene casi el doble de diapositivas aunque éstas son más conceptuales y mucho más visuales: Toros Invariantes KAM y Huecos Resonantes en el Sistema Solar (Seminario)

Aunque para entender propiamente el tema se requiere conocer la teoría de Hamilton-Jacobi y por supuesto el Teorema KAM, ambas presentaciones son más bien cualitativas; a lo mucho se muestran las variables acción y de ahí las frecuencias para el problema de tres cuerpos restringido circular. En Mathematica ilustramos, con ayuda de Wolfram Demonstrations, los puntos de Lagrange y un mapeo desarrollado por Jan Frøyland para los anillos de Saturno.


Para los puntos de Lagrange, se considera el potencial efectivo para el problema de tres cuerpos restringido circular
\begin{equation}\Omega\equiv\frac{1}{2}(x^2+y^2)+\frac{1-\mu}{r_1}+\frac{\mu}{r_2}\end{equation} y de ahí se encuentran las raíces para la fuerza. En este documento que escribí hace un tiempo puedes encontrar más detalles: El problema de $n$-cuerpos en la entrada del mismo nombre.

(*** DEFINICION DEL POTENCIAL EFECTIVO ***)
V[x_,y_,M_]:=1/2 (x^2+y^2)+(1-M)/Sqrt[(x+M)^2+y^2]+M/Sqrt[(x+M-1)^2+y^2]
M=1/12;

(*** ENCUENTRA PUNTOS DE LAGRANGE Y GRAFICA CON EL MAPA DE CONTORNO DE V ***)
Module[{sols, deriv, L1, L2, L3, L4, L5},
deriv = D[V[x, 0,M], x];
L1 = FindRoot[deriv == 0, {x, 0}];
L2 = FindRoot[deriv == 0, {x, 1.5}];
L3 = FindRoot[deriv == 0, {x, -1}];
L4 = FindRoot[{D[V[x,y,M] ,x]== 0, D[V[x,y,M],y] == 0}, {{x, 0}, {y, 1}}];
L5 = FindRoot[{D[V[x,y,M],x] == 0, D[V[x,y,M],y] == 0}, {{x, 0}, {y, -1}}];
sols = {{x /. L1, 0,V[x /.L1, 0,M]},
{x /. L2, 0, V[x /.L2, 0,M]},
{x /. L3, 0, V[x /.L3, 0,M]},
{x /. L4, y /.L4,V[x /.L4, y /.L4,M]},
{x /.L5, y /.L5, V[x /.L5, y /.L5,M]}};

ContourPlot[
V[x, y,M], {x, -2, 2}, {y, -2, 2}, MaxRecursion -> Automatic, 
Frame -> False, PerformanceGoal->"Quality", ImageSize->500,
ContourShading->None, ContourStyle->Black, Contours->20,
Epilog -> {PointSize[.025], Red, Point[Most[#] & /@ sols]}]]

Para el modelo de Frøyland, se describe la distribución de partículas entre los radios ${r_0\approx185.5\times10^3\,\mathrm{km}}$ (distancia del centro de Saturno a Mimas, el satélite perturbativo) y ${r_s\approx60.5\times10^3\,\mathrm{km}}$ (radio de la superficie de Saturno). De la tercera ley de Kepler, se sabe que una partícula de masa $m$ en órbita circular de radio $r$ alrededor de Saturno tendrá un periodo $\tau$ tal que ${\tau^2=\frac{4\pi^2}{Gm}r^3}$. Como radio de referencia se considera ${r_0}$, de modo que cada vez que Mimas completa una órbita alrededor de Saturno, su posición angular cambia en ${2\pi}$, mientras que para una partícula en ${r_0<{r}_n<{r}}$ tras la $n$-ésima órbita, con periodo ${\tau_n}$, su posición angular $\theta_n$ variará ligeramente de la de Mimas. Frøyland entonces describe el ángulo ${\theta_{n+1}}$ módulo ${2\pi}$ para la ${(n+1)}$-órbita vía \begin{equation}\theta_{n+1}=\theta_n+2\pi\left(\frac{\tau_0}{\tau_n}\right)=\theta_n+2\pi\left(\frac{r_0}{r_n}\right)^{3/2}\end{equation} Ésta es la mitad del mapeo, pues la posición radial también debería cambiar debido al efecto perturbativo de Mimas. Frøyland considera la solución de la ecuación ${\ddot{r}=\mc{F}_n}$ con $\mc{F}_n$ fuerza radial por unidad de masa (segunda ley de Newton) por diferencias finitas, promediando la aceleración en un periodo completo de Mimas (${\Delta{t}=\tau_0}$)
\begin{equation}\mc{F}_n=\frac{r_{n+1}-2r_n+r_{n-1}}{\tau_0^2}\end{equation} también, por método de Euler, ${\mc{F}_n(r_n,\theta_n)=f_n(r_n,\theta_n)/\tau_0^2}$, de donde propone
\begin{equation}f_n=-A\frac{\cos\theta_n}{(r_0-r_n)^2}\end{equation} con $A$ una constante positiva. Presuntamente la forma de ${f_n}$ es mucho más complicada si se consideran los efectos de otras lunas además de Mimas, por ello se permite un rango bastante amplio para $A$. De aquí se sigue finalmente que
\begin{equation}r_{n+1}=2r_n-r_{n-1}-A\frac{\cos\theta_n}{(r_0-r_n)^2}\end{equation} Nos hemos enterado primero de este trabajo en Orden y caos en sistemas complejos, que desafortunadamente no obtiene el mapa de Frøyland, sino que únicamente lo presenta. No he logrado encontrar el trabajo original de Frøyland, ni trabajos posteriores al respecto, e. g. de Gould y Tobochnik, sin embargo la referencia en donde encontré básicamente el razonamiento (un tanto de handwaving, como se dice en inglés) que aquí muestro es It's a nonlinear world de Richard Enns, que me parece suficiente para saber cómo funciona la cosa, a menos que uno decida aventarse un clavado más profundo en el tema.

El mapeo de Frøyland puede visualizarse como sigue en Mathematica.

(*** NUMERO DE ITERACIONES Y VALOR ARBITRARIO DE LA CONSTANTE A ***)
nit=4000;
A=180;

(* MAPEO EN COORDENADAS POLARES. j->Condiciones iniciales *)
Quiet[With[{Puntos=Map[#1[[1]] {Cos[#1[[2]]],Sin[#1[[2]]]}&,
Table[
RecurrenceTable[
{r[n+1]==2 r[n]-r[n-1]-(A Cos[T[n]])/(r[n]-185.7)^2,
T[n+1]==T[n]+2Pi (185.7/r[n])^(3/2),
r[0]==66+3 j, r[1]==r[0], T[0]==0}, {r,T}, {n,1,nit}],
{j,40}],
{2}]},
(* DIBUJA LOS PUNTOS DEL MAPEO *)
ListPlot[Puntos, AspectRatio->1, Axes->False,
PlotStyle->{PointSize[.001]}, PlotRange->180{{-1,1},{-1,1}},
ImageSize->500 ,
(* PLANETA *)
Epilog->{RGBColor[0,0,0], Disk[{0,0},60]}]]]

Interferencia en Rendijas Múltiples

Considérese una fuente de luz monocromática que pasa a través de tres rendijas paralelas separadas entre sí por una distancia $d$. Por simplicidad además, considérese que las ondas tienen la misma amplitud $\mathcal{E}$, la misma longitud de onda $\lambda$, y así la misma frecuencia angular $\omega$ y una diferencia de fase constante ${\phi=\frac{d\sin\theta}{\lambda}}$ donde ${\theta\ll1}$ es el ángulo entre la normal a las rendijas y el vector al punto de incidencia.

Se tiene entonces que las tres ondas que emergen de las rendijas son de la forma
\begin{align}\mathbf{E}_1&=\boldsymbol{\mathcal{E}}\sin\omega{t}\\
\mathbf{E}_2&=\boldsymbol{\mathcal{E}}\sin(\omega{t}+\phi)\\
\mathbf{E}_3&=\boldsymbol{\mathcal{E}}\sin(\omega{t}+2\phi)\end{align} Para conocer la intensidad, se necesita calcular el promedio temporal
\begin{equation}\left\langle\mathbf{E}^2\right\rangle=\frac{1}{\tau}\int_{t}^{t+\tau}|\mathbf{E}(t^\prime)|^2\,dt^\prime\end{equation} donde ${\mathbf{E}=\sum_i\mathbf{E}_i}$ y el periodo $\tau$, en general para funciones armónicas es ${2\pi/\omega}$.

Para lograrlo entonces, hay que calcular la norma de la resultante del campo eléctrico al cuadrado. Se tiene que
\begin{equation}\mathbf{E}=\boldsymbol{\mathcal{E}}\left[\sin\omega{t}+\sin(\omega{t}+\phi)+\sin(\omega{t}+2\phi)\right]\end{equation} y empleando la identidad
\begin{equation}\sin{A}+\sin{B}=2\cos\left(\frac{A-B}{2}\right)\sin\left(\frac{A+B}{2}\right)\end{equation} puede verse que es conveniente realizar la suma
\begin{equation}\mathbf{E}_1+\mathbf{E}_3=2\boldsymbol{\mathcal{E}}\cos\left(\phi\right)\sin(\omega{t}+\phi)\end{equation} de modo que
\begin{align}\mathbf{E}&=\boldsymbol{\mathcal{E}}\left[\sin(\omega{t}+\phi)+2\cos(\phi)\sin(\omega{t}+\phi)\right]\nonumber\\&=\boldsymbol{\mathcal{E}}(1+2\cos\phi)\sin(\omega{t}+\phi)\end{align} donde es inmediato realizar el promedio,
\begin{align}I&\propto\left\langle\mathbf{E}^2\right\rangle=\mathcal{E}^2\frac{(1+2\cos\phi)^2\omega}{2\pi}\int_0^{2\pi/\omega}\sin^2(\omega{t}+\phi)\,dt\nonumber\\&=\mathcal{E}^2\frac{(1+2\cos\phi)^2\omega}{4\pi}\left(\frac{2\pi}{\omega}\right)\nonumber\\&=\mathcal{E}^2\frac{(1+2\cos\phi)^2}{2}\end{align} i.e. simplemente
\begin{equation}\left\langle\sin^2(\omega{t}+\phi)\right\rangle=\frac{1}{2}\label{ast}\end{equation} De aquí entonces puede verse que la intensidad máxima ${I_m}$ ocurre cuando ${\cos\phi=1}$, i.e.
\begin{equation}I_m\propto\frac{9}{2}\mathcal{E}^2\end{equation} entonces puede escribirse
\begin{equation}\frac{I}{I_m}=\frac{(1+2\cos\phi)^2}{9}\end{equation} es decir, explícitamente
\begin{equation}\frac{I}{I_m}=\frac{1}{9}\left[1+2\cos\left(\frac{d\sin\theta}{\lambda}\right)\right]^2\end{equation} con esto entonces uno puede graficar el comportamiento de la intensidad respecto a la fase $\phi$, e.g. en Mathematica:

Plot[(1 + 2 Cos[2Pi F])^2/9, {F, -Pi, Pi},
 Axes -> None, Frame -> True, FrameLabel -> {"2Pi F", "I/Im"},
 FrameTicks -> {{{0, 1/2, 1}, None}, {{-Pi, -Pi/2, 0, Pi/2, Pi}, None}},
 GridLines -> Automatic, GridLinesStyle -> Directive[Gray, Dashed]]

La generalización natural de la triple rendija es considerar el caso de n rendijas. Considérense entonces ahora, de manera análoga, las n ondas
\begin{equation}\mathbf{E}_\alpha=\boldsymbol{\mathcal{E}}\sin\left(\omega{t}+\alpha\phi\right),\hspace{0.25in}\alpha=0,1,\ldots,n-1\end{equation} en cuyo caso entonces, debe calcularse el promedio ${\langle\mathbf{E}^2\rangle}$ con
\begin{equation}\mathbf{E}=\boldsymbol{\mathcal{E}}\sum_{\alpha=0}^{n-1}\sin(\omega{t}+\alpha\phi)\end{equation} entonces recordando de la serie geométrica que $\displaystyle{\sum_{a=0}^{b-1}r^a=\frac{1-r^b}{1-r}}$, se tiene,
\begin{align}\sum_{\alpha=0}^{n-1}\sin(\omega{t}+\alpha\phi)&=\sum_{\alpha=0}^{n-1}\left(\sin{\omega{t}}\cos\alpha\phi+\sin\alpha\phi\cos{\omega{t}}\right)\nonumber\\[0.1in]
&=\sin{\omega{t}}\,\Re\left[\sum_{\alpha=0}^{n-1}\mathrm{e}^{i\alpha\phi}\right]+\cos{\omega{t}}\,\Im\left[\sum_{\alpha=0}^{n-1}\mathrm{e}^{i\alpha\phi}\right]\nonumber\\[0.1in]
&=\sin{\omega{t}}\,\Re\left[\frac{1-\mathrm{e}^{in\phi}}{1-\mathrm{e}^{i\phi}}\right]+\cos{\omega{t}}\,\Im\left[\frac{1-\mathrm{e}^{in\phi}}{1-\mathrm{e}^{i\phi}}\right]\nonumber\\[0.1in]
&=\sin{\omega{t}}\,\Re\left[\frac{\mathrm{e}^{in\phi/2}\left(\mathrm{e}^{-in\phi/2}-\mathrm{e}^{in\phi/2}\right)}{\mathrm{e}^{i\phi/2}(\mathrm{e}^{-i\phi/2}-\mathrm{e}^{i\phi/2})}\right]+\cos{\omega{t}}\,\Im\left[\frac{\mathrm{e}^{in\phi/2}\left(\mathrm{e}^{-in\phi/2}-\mathrm{e}^{in\phi/2}\right)}{\mathrm{e}^{i\phi/2}(\mathrm{e}^{-i\phi/2}-\mathrm{e}^{i\phi/2})}\right]\nonumber\\[0.1in]
&=\sin{\omega{t}}\,\Re\left[\mathrm{e}^{i(n-1)\phi/2}\frac{\sin\left(n\frac{\phi}{2}\right)}{\sin\left(\frac{\phi}{2}\right)}\right]+\cos{\omega{t}}\,\Im\left[\mathrm{e}^{i(n-1)\phi/2}\frac{\sin\left(n\frac{\phi}{2}\right)}{\sin\left(\frac{\phi}{2}\right)}\right]\nonumber\\[0.1in]
&=\sin{\omega{t}}\frac{\sin\left(\frac{n\phi}{2}\right)\cos\left[\frac{(n-1)\phi}{2}\right]}{\sin\left(\frac{\phi}{2}\right)}+\cos{\omega{t}}\frac{\sin\left(\frac{n\phi}{2}\right)\sin\left[\frac{(n-1)\phi}{2}\right]}{\sin\left(\frac{\phi}{2}\right)}\nonumber\\[0.1in]
&=\frac{\sin\left(\frac{n\phi}{2}\right)}{\sin\left(\frac{\phi}{2}\right)}\left\{\sin{\omega{t}}\cos\left[\frac{(n-1)\phi}{2}\right]+\cos{\omega{t}}\sin\left[\frac{(n-1)\phi}{2}\right]\right\}\nonumber\\[0.1in]
&=\csc\left(\frac{\phi}{2}\right)\sin\left(\frac{n\phi}{2}\right)\sin\left[\omega{t}+\frac{n-1}{2}\phi\right]\end{align} y por tanto, empleando el resultado (\ref{ast}), se tiene simplemente que
\begin{equation}I\propto\left\langle\mathbf{E}^2\right\rangle=\frac{\mathcal{E}^2}{2}\csc^2\left(\frac{\phi}{2}\right)\sin^2\left(\frac{n\phi}{2}\right)\end{equation} donde la amplitud máxima ahora depende del valor de n. Uno puede investigar fácilmente cómo se relaciona la amplitud máxima graficando la función ${\csc{x}\sin{nx}}$,

donde se hace evidente que la amplitud máxima será
\begin{equation}I_m\propto\frac{\mathcal{E}^2}{2}n^2\end{equation} lo que concuerda con lo hallado para 3 rendijas. Así entonces
\begin{equation}\frac{I}{I_0}=\frac{1}{n^2}\csc^2\left(\frac{\phi}{2}\right)\sin^2\left(\frac{n\phi}{2}\right)\end{equation} y nuevamente uno puede visualizar los resultados para un distinto número de rendijas, e.g. con Mathematica,

Trazando las curvas de nivel de los campos de radiación

Apenas va terminando mi malísimo curso de Radiación y Óptica, y sin embargo he rescatado algunas cosas al trabajar por mi cuenta. Lo primero fue graficar el mapa de contorno de la componente de los campos eléctrico $\B{E}$ y/o magnético $\B{B}$ para un dipolo eléctrico oscilante. Lograrlo de hecho fue sencillísimo, ya que uno llega a las expresiones del tipo
\begin{equation}\B{E}=-\alpha\frac{\sin\theta}{r}\cos\omega\tau\,\boldsymbol{\hat{\theta}},\hspace{0.75in}\B{B}=-\frac{\alpha}{c}\frac{\sin\theta}{r}\cos\omega\tau\,\boldsymbol{\hat{\varphi}}\end{equation} con $\alpha$ constante, $\theta$ el ángulo polar, $\varphi$ el ángulo azimutal y ${\tau=t-r/c}$ el tiempo de retardo, entonces uno puede simplemente graficar curvas de nivel para la correspondiente componente angular manteniendo algún parámetro dado fijo, e.g. $t$, y haciendo ${\alpha=1}$ por simplicidad,

Table[ContourPlot[-(Sin[ArcCos[z/Norm[{x, y, z}]]]/Norm[{x, y, z}])
 Cos[t - Norm[{x, y, z}]] /. {t -> Pi}, {x, -25, 25}, {z, -25, 25},
  MaxRecursion -> 5, ContourShading -> None, FrameLabel -> {x, z}, ContourStyle -> Black, 
   PlotLabel -> "y=" ~~ ToString[y]], {y, 0, 9, 1}]

o bien, para $y$ fijo y variando $t$,

GIF = Table[
   Manipulate[
    ContourPlot[-(Sin[ArcCos[z/Norm[{x, y, z}]]]/Norm[{x, y, z}]) Cos[
        t - Norm[{x, y, z}]] /. {y -> 0}, {x, -10, 10}, {z, -10, 10}, 
     Contours -> 10, ContourShading -> None, ContourStyle -> Black, 
     Exclusions -> {x == 0}, FrameTicks -> None, MaxRecursion -> 4],
     {t, k, Pi, ControlType -> None}], {k, 0, 9Pi/10, Pi/10}];

Export["Campo.gif", GIF, "DisplayDurations" -> 0.2]

Sin embargo, uno bien pudo haber decidido pasar el campo (ya sea $\B{E}$ o $\B{B}$) a coordenadas cartesianas y graficar en dichas coordenadas alguna componente arbitraria, esto es, por ejemplo para el campo $\B{E}$, sabiendo que
\begin{equation}\boldsymbol{\hat{\theta}}=\begin{pmatrix}\cos\varphi\cos\theta\\\sin\varphi\cos\theta\\-\sin\theta\end{pmatrix}\end{equation} además, en los anteriores gráficos he hecho explícitamente ${\theta=\arccos\frac{z}{r}}$, por lo que de manera análoga ahora se sustituye explícitamente ${\varphi}$ en términos de ${x,y}$ (véase atan2); y uno puede llevarse una sorpresa al querer graficar las 3 componentes del campo como se hizo con la componente angular, por ejemplo con la componente $x$, uno obtiene

aunque en cierto modo era de esperarse, pues se trata de la componente $x$ de un campo que sólo cambia en la dirección polar graficada en el plano ${\{x,z\}}$; de hecho la componente cartesiana más parecida a la componente polar es la componente $z$, como también es de esperar, por ello resulta difícil interpretar cualitativamente el gráfico en este modo, mientras que es muy sencillo hacerlo con el gráfico de la componente polar.

Esto lo digo porque después de estudiar el caso del dipolo, estudié el caso de una carga en movimiento circular uniforme (clásico), y la visualización a primeras no arrojó nada cualitativamente bueno, presuntamente por el detalle de las coordenadas que menciono aquí. Uno puede ir y encontrar los potenciales explícitos de Liénard–Wiechert para el problema, que son de la forma
\begin{align}V(\B{r},\tau)&=\frac{1}{\sqrt{r^2+\rho^2-2\rho\left(x\cos\tau+
y\sin\tau\right)}-\rho\left(y\cos\tau-x\sin\tau\right)}\\
\B{A}(\B{r},\tau)&=\rho\,V(\B{r},\tau)\,\boldsymbol{\hat{\varphi}}\end{align} donde simplemente tomé todas las constantes como la unidad, excepto el radio de la órbita clásica ${|\boldsymbol{\rho}|=\rho}$, que puse en el plano ${\{x,y\}}$ y donde ${\tau=t-\frac{|\B{r}-\boldsymbol{\rho}|}{c}}$ con ${\B{r}=(x,y,z)}$. De aquí uno entonces puede obtener los campos vía
\begin{equation}\B{E}=-\nabla{V}-\p_t\B{A},\hspace{0.75in}\B{B}=\nabla\times\B{A}\end{equation} tomando antes, por supuesto, el tiempo de retardo $\tau$ explícitamente en función de $t$, y haciendo todo de una buena vez en coordenadas esféricas, i.e. con ${x=r\cos\varphi\sin\theta}$, ${y=r\sin\varphi\sin\theta}$ y el gradiente y el rotacional en esféricas. Finalmente al graficar, debe regresarse a las variables cartesianas, aunque ya se podrá graficar cualquier componente esférica. De manera análoga uno puede partir directamente de los campos de radiación sin pasar por los potenciales, si uno cuenta con las expresiones explícitas.

El cómputo es notablemente caro con un ordenador promedio, por lo que hay que tener algo de paciencia. Lo ideal sería generar una imagen .gif para cada componente para poder interpretar claramente cada componente, como hice con el dipolo, pero para eso haría falta bastante tiempo o un ordenador más rápido. De cualquier modo no pienso quitarle al lector la diversión de hacerlo por su cuenta, por lo que solo comparto una de las salidas para la componente angular del campo $\B{E}$, la que exhorto a verificar, pues meter la pata puede ser bastante fácil, además aparentemente no gané mucho al pasarme a las componentes del campo en coordenadas esféricas, pues las curvas de nivel son muy parecidas a las de las componentes cartesianas. En color rojo marco la posición de la partícula, que describe una órbita circular de radio 3.

Péndulo con soporte en una parábola oscilante

Este es un problema bastante divertido, cuya solución es análoga a la del péndulo doble. El punto de suspensión de masa $M$ de un péndulo simple de longitud $\ell$ y masa $m$ está restringido a moverse sobre una parábola oscilante dada por \begin{equation}y=\alpha{x}^2+\sin\omega{t}\end{equation} en el plano vertical.

Lo que se quiere es

1. Obtener el Lagrangiano y el Hamiltoniano del sistema.
2. Obtener las ecuaciones de movimiento de Lagrange y de Hamilton.
3. Resolver las ecuaciones de movimiento y visualizar las soluciones.

Para la descripción Lagrangiana se tiene

• Posición del soporte: $\vec{R}=\begin{cases}X=X(t)\\Y=\alpha{X}^2+\sin\omega{t}\end{cases}$
• Posición de la masa del péndulo: $\vec{r}=\begin{cases}x=X+\ell\sin\theta\\y=Y-\ell\cos\theta\end{cases}$
• Elección de coordenadas generalizadas: $\{q_1,q_2\}=\{X(t),\theta(t)\}$
• Energía Cinética: $T=\frac{1}{2}m\left(\dot{x}^2+\dot{y}^2\right)+\frac{1}{2}M\left(\dot{X}^2+\dot{Y}^2\right)=T\left(X,\dot{X},\theta,\dot{\theta},t\right)$
• Energía Potencial: $V=g\left(my+MY\right)=V\left(X,\theta,t\right)$
• Lagrangiana: $\mathcal{L}\equiv{T-V}=\mathcal{L}\left(X,\dot{X},\theta,\dot{\theta},t\right)$

Para pasar a la descripción Hamiltoniana

• Las fuerzas en el sistema pueden derivarse de $V$
• $\mathcal{L}=\mathcal{L}\left(\vec{q},\dot{\vec{q}},t\right)\,\Longrightarrow$ No hay coordenadas cíclicas $\Longrightarrow$ No se conserva cantidad alguna
• $\frac{\partial\vec{R}}{\partial{t}}\neq\vec{0},\,\frac{\partial\vec{r}}{\partial{t}}\neq\vec{0}\;\Longrightarrow\,\mathcal{H}\neq{T+V}=E$
• $\frac{\partial{E}}{\partial{t}}\neq{0}\,\Longrightarrow\,E=E(t)$
• Se puede obtener el Hamiltoniano directamente de la definición por transformada de Legendre:
  \begin{align}\mathcal{H}&\equiv\sum_ip_i\dot{q}_i-\mathcal{L}\nonumber\\
  &=p_{_X}\dot{X}+p_{_\theta}\dot{\theta}-\mathcal{L}\left(X,\dot{X},\theta,\dot{\theta},t\right)\end{align} donde
  \begin{equation}p_{_j}=\frac{\partial\mathcal{L}}{\partial\dot{q}_{_j}}\;\Longrightarrow\;\dot{q}_{_j}=\dot{q}_{_j}(q_{_j},p_{_j},t)\;\Longrightarrow\;\mathcal{H}=\mathcal{H}\left(X,p_{_X},\theta,p_{_\theta},t\right)\end{equation}
• O se puede obtener de la forma
  \begin{equation}\mathcal{H}=\frac{1}{2}\left(\vec{p}-\vec{b}\right)^\mathrm{T}\mathbb{M}^{-1}\left(\vec{p}-\vec{b}\right)
  -\mathcal{L}_0\end{equation}

Para lograr el punto 3., preferí utilizar Mathematica

Mi archivo en Mathematica se ve como sigue (da clic derecho + Ver Imagen para ver el tamaño completo):

Instalando Mathematica 9 en Ubuntu

Ya he confesado antes que Mathematica ha sido el único software que no he podido reemplazar por alguna versión libre. Apenas he reinstalado mi sistema y he recordado que el proceso de instalación puede ser un tanto desconcertante para quien se aventura por vez primera a utilizar un sistema basado en Linux y que no necesariamente es un geek aficionado. El proceso se reduce básicamente a ejecutar un archivo file.sh, lo que no tiene mayor problema, sólo se dan permisos para ejecutar

chmod +x 'directorio original/file.sh'

y entonces se ejecuta el archivo,

./'directorio original/file.sh'

El único detalle que puede hacer fallar la instalación, es que primero hay que copiar o trasladar el archivo a la carpeta /opt, en donde se almacenan paquetes externos,

sudo mkdir '/opt/Mathematica 9'
sudo cp -r 'directorio original/file.sh' '/opt/Mathematica 9/file.sh'

Y listo, sólo hay que ejecutar el archivo y concluir la instalación. Esto no sucede en general con los archivos .sh, que hasta donde sé, pueden tener formas y finalidades muy diversas, en general uno está acostumbrado a instalar paquetes .deb, pero bueno, esto podría ser útil con cualquier caso parecido al de Mathematica.

Condicionales en Mathematica

Lo que el gráfico muestra es una función
$$f(x)=\left\{\begin{array}{ll}16\,\mathrm{e}^{-x^2/3},&x<0\\[0.1in]\lfloor{x}\rfloor,&0\leq{x}<8\\[0.1in]25-x,&8\leq{x}<14\\[0.1in]5\sin{x},&x\geq14\end{array}\right.$$ lograda no con Piecewise, sino con un vulgar y silvestre If (anidado):

Plot[If[x > 0, If[x < 8, Floor[x], If[x < 14, -x + 22, 5 Sin[x]]],
16 Exp[-1/3 x^2]], {x, -5, 7 Pi}, PlotStyle -> {Red, Thick},
Exclusions -> {x == 0, x == 8, x == 14}, AxesLabel -> {x, f[x]}]

Esto surge por una pregunta que me encontré en Mathematica SE (por la que decidí por fin registrarme). Se pregunta cómo evaluar sumatorios con condiciones, por ejemplo
$$\sum_{\substack{i=-\infty\\i\neq0}}^{\infty}\,\sum_{\substack{j=2\\j|i\\j\neq{i}}}^{n}f(i,j)$$ puede lograrse con una instrucción del tipo

Sum[If[Divisible[i, j] && i != 0 && i != j, f[i, j], 0],
{i, -Infinity, Infinity}, {j, 2, n}]

Lo que hace If[cond, V, F] es regresar V si cond es cierta o F si cond es falsa. Ésta, y en general los condicionales, son funciones bien conocidas por quienes gustan programar, y en Mathematica a veces puede no ser obvio el cómo implementarlas. Acá se muestran los condicionales comunes en Mathematica, y como se dice, algunos resultan más económicos que un If anidado. Curiosamente en SE la mejor respuesta no fue la mía, sino una que proponía utilizar Boole. Y pues en efecto, Boole[A] regresa 1 si A es cierta o 0 si A es falsa. La suma anterior entonces se escribiría como

Sum[f[i, j] Boole[i != 0] Boole[Divisible[i, j]] Boole[i != j],
{i, -Infinity, Infinity}, {j, 2, n}]

La salida, por ejemplo, para ${-10\leq{i}\leq10,\,2\leq{j}\leq10}$ de esta suma, es

f[-10, 2] + f[-10, 5] + f[-10, 10] + f[-9, 3] + f[-9, 9] + f[-8, 2] +
f[-8, 4] + f[-8, 8] + f[-7, 7] + f[-6, 2] + f[-6, 3] + f[-6, 6] +
f[-5, 5] + f[-4, 2] + f[-4, 4] + f[-3, 3] + f[-2, 2] + f[4, 2] +
f[6, 2] + f[6, 3] + f[8, 2] + f[8, 4] + f[9, 3] + f[10, 2] + f[10, 5]

Tal vez es sólo cuestión de gustos, pero en general ésta me parece la forma más sencilla de usar los condicionales en Mathematica, simplemente incorporarlos al argumento de la función en cuestión.

Partícula en un cilindro

Jocosamente, la situación de la partícula en una caja es clásica en mecánica cuántica. Comparto el caso de una caja cilíndrica, sobre todo por la parte de la visualización de los resultados en Mathematica, que no tuve oportunidad de hacer con más cuidado dentro de mi primer curso de cuántica.

Se tiene una partícula confinada en una caja cilíndrica de radio $\mathcal{R}$ y altura $h$, esto es, en coordenadas cilíndricas ${(r,\theta,z)}$, la partícula está sujeta al potencial
$$V=\left\{\begin{array}{l}0\,\,\text{si}\,\,0\leq{z}\leq{h},\,0\leq{r}\leq\mathcal{R}\\\infty\,\,\text{de otro modo}\end{array}\right.$$ entonces dentro del cilindro la ecuación estacionaria de Schrödinger es
$$-\frac{\hbar^2}{2m}\nabla_{r\theta{z}}^2\psi=E\psi$$ con $\nabla_{r\theta{z}}^2$ el Laplaciano en coordenadas cilíndricas, y cuya solución puede hallarse por el método de variables separables, i.e. es de la forma ${\psi(r,\theta,z)=R(r)\Theta(\theta)Z(z)}$. Encontrar explícitamente estas funciones es precisamente el problema a resolver en los cursos básicos de cuántica, así que no lo mostraré aquí. Para normalizar la función de onda te puede ser útil esta entrada. Finalmente se llega a que los eigenvalores, i.e. los niveles de energía, están dados por
$$E=\frac{\hbar^2}{2m}\left[\left(\frac{\gamma\,\pi}{h}\right)^2+\lambda_k^2\right]$$ donde ${\gamma=1,2,\ldots}$ y $\lambda_k$ satisface ${J_\mu(\lambda_k\mathcal{R})=0}$ en el k-ésimo cero para ${\mu=0,\pm1,\pm2,\ldots}$ con ${J_\mu}$ funciones Bessel de primera especie y orden $\mu$ (todo esto sólo se hace evidente resolviendo el ejercicio uno mismo), mientras que las eigenfunciones,
$$\psi_{\gamma\mu{k}}(r,\theta,z)=N_{\gamma\mu{k}}J_\mu(\lambda_k\,r)\mathrm{e}^{i\mu\theta}\sin\left(\frac{\gamma\pi}{h}\,z\right)$$ donde ${N_{\gamma\mu{k}}=\pm\left(\frac{2}{h\,\pi}\right)^{1/2}\left[\mathcal{R}\,J_\mu^\prime(\lambda_k\mathcal{R})\right]^{-1}}$ es la constante de normalización.

Para tener concretamente los niveles de energía únicamente se requiere el k-ésimo cero de la función ${J_\mu}$, mismo que será necesario para visualizar la función de onda (o en general la distribución de probabilidad), para la cual se puede graficar e.g. la parte real de diversas superficies de nivel en el eje $z$ en las mismas coordenadas cilíndricas. Una forma de hacer esto en Mathematica con valores arbitrarios para el radio, la altura y los números cuánticos, es la siguiente:

BesselJZero[n,k] encuentra el k-ésimo cero de la función Bessel de primera especie y orden n, BesselJ[n,x] da la función Bessel de primera especie y orden n en x. Utilizo además ParametricPlot3D[] para graficar en coordenadas cilíndricas (con la parametrización correspondiente). Puedes también encerrar todo el comando de los gráficos en la función Timing[] para saber cuánto tarda la evaluación. Ahora bien, utilizo Table[], y al menos en mi ordenador la evaluación es extremadamente lenta (o el código es ineficiente, probablemente esto sea por manejar los números cuánticos como parámetros, siendo sincero por ahora ignoro si hay una mejor manera de hacerlo), de aprox 8 minutos, y no de muy buena calidad (si probara dando un valor más grande a MaxRecursion[] probablemente se pasaría todo el día evaluando), si en el tuyo se ejecuta más rápidamente puedes probar con Manipulate[] para obtener una mejor visualización. La salida del código es la siguiente

viendo de cerca uno de los casos, con MaxRecursion->5,

Con un poco más de potencia computacional podría hacerse un gráfico de un cilindro y la (parte real y/o imaginaria) función de onda y se tendría una muestra muy mona de la solución del problema, casi como este programita de Wolfram Demonstrations, pero con superficies de nivel en lugar de curvas de nivel; además de poderse emplear Manipulate[] para visualizar cómo cambia la función de onda "en tiempo real".

De aquí no debe haber ningún problema para visualizar la densidad de probabilidad en un espacio análogo. Quizá de lo más ilustrativo de estas visualizaciones es el comportamiento respecto a los números cuánticos; conforme decidí jugar con esto, por ejemplo, me di cuenta que en mi curso de cuántica se consideró ${\lambda_k}$ por sí mismo como un número cuántico, cuando de manera precisa, k (que denota el k-ésimo cero de las funciones Bessel) es el número cuántico. Prueba dando distintos valores a los números cuánticos y observa cómo afectan en la solución.

Excluir asíntotas en Mathematica

Seguido es necesario graficar funciones con discontinuidades esenciales, y en Mathematica eso seguido significa ver las líneas verticales en los puntos de discontinuidad, por ejemplo, al introducir el comando

Plot[Tan[x],{x,-2Pi,2Pi}, PlotStyle -> Thick]

uno obtiene

Bueno, pues uno se deshace de estas líneas, o en general de cualquier punto, con la función Exclusions[ ] de la cual puedes leer detalles aquí o en la misma sección de ayuda de Mathematica. El uso del comando es bastante sencillo, por ejemplo, con la línea

Plot[Tan[x],{x,-2Pi,2Pi}, PlotStyle -> Thick, Exclusions -> {Cos[x] == 0}]

obtenemos ya el gráfico deseado excluyendo aquellos puntos en los que cos(x)=0, i.e. en los que la función no está definida

De manera análoga para funciones con saltos como Floor[ ] o HeavisideTheta[ ] en las que Mathematica no muestra las líneas verticales, puede escribirse, por ejemplo

Plot[Floor[x], {x, 0, 10}, PlotStyle -> Thick, Exclusions -> None]

para obtenerlas de vuelta. El comando puede ser de utilidad en diversas situaciones, yo en particular he notado su utilidad al visualizar las soluciones que llevan a los estados de energía permitidos para una función de onda cuántica unidimensional, sin las estorbosas asíntotas de unas funciones cotangente y tangente (por eso el nombre de la entrada, particularmente para excluir asíntotas), y la diferencia es sustancial, uno pasa de visualizar algo como esto

a algo como esto

pero bueno, en Mathematica las posibilidades parecen ser inagotables, por lo que muy probablemente te llegue a ser útil esta función en diversas situaciones.

El péndulo doble por formalismo de Lagrange

El péndulo doble es un problema que puede simularse fácilmente utilizando el formalismo de Lagrange o de Hamilton. En realidad la generalización a partir de aquí se sigue de forma muy sencilla, de modo que pueden formarse sistemas de muchos péndulos más. Además por supuesto, siempre puede hacerse más específico el problema para alguna situación dada.

Acá muestro el caso más sencillo en formalismo de Lagrange para el péndulo doble en un campo gravitacional constante, en el cual ambas masas son iguales y las longitudes de los trozos de cuerda que las unen con sus orígenes son iguales:


Se tienen dos grados de libertad, ya que cada partícula se ha restringido al plano ${\{x,y\}}$ y a mantener una longitud máxima $\ell$ con su respectivo origen. Así pues, a partir de la figura, se eligen por conveniencia, las coordenadas generalizadas generalizadas, ${\{\theta_1,\;\theta_2\}}$.

Sean $(x_1,y_1)$, $(x_2,y_2)$ las posiciones de las masas ${m_1}$ y ${m_2}$, respectivamente (se sabe que ${m_1=m_2=m}$, por simplicidad, utilizo la notación sólo para distinguir una de otra), entonces
\begin{equation}\begin{array}{ll}x_1=\ell\sin\theta_1&\hspace{0.5in}x_2=\ell\sin\theta_2+x_1\\y_1=-\ell\cos\theta_1&\hspace{0.5in}y_2=-\ell\cos\theta_2+y_1\end{array}\end{equation} y de este modo, se tienen la energía cinética $T$ y la energía potencial $V$ del sistema,
\begin{align}T&\equiv\frac{1}{2}\sum_im_i\left(\dot{x}_i^2+\dot{y}_i^2\right)\nonumber\\&=\frac{1}{2}m\ell^2\left[2\dot{\theta}_1^2+\dot{\theta}_2^2+2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2)\right]\\[0.25in]V&\equiv\sum_im_igy_i\nonumber\\&=-mg\ell(2\cos\theta_1+\cos\theta_2)\end{align} y por tanto, la Lagrangiana del sistema es
\begin{align}\mathcal{L}&\equiv{T}-V\nonumber\\&=\frac{1}{2}m\ell^2\left[2\dot{\theta}_1^2+\dot{\theta}_2^2+2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2)\right]+mg\ell(2\cos\theta_1+\cos\theta_2)\end{align} de donde, ya que $\displaystyle{\frac{\partial\mathcal{L}}{\partial{t}}=0}$, se conserva el Hamiltoniano $\mathcal{H}$ del sistema, esto es
\begin{align}\mathcal{H}&\equiv\sum_i\frac{\partial\mathcal{L}}{\partial\dot\theta_i}\dot\theta_i-\mathcal{L}\nonumber\\[0.1in]&=m\ell^2\dot{\theta}_1\left[2\dot{\theta}_1+\dot{\theta}_2\cos(\theta_1-\theta_2)\right]+ml^2\dot{\theta}_2\left[\dot{\theta}_2+\dot{\theta}_1\cos(\theta_1-\theta_2)\right]-T+V\nonumber\\[0.1in]&=m\ell^2\left[2\dot{\theta}_1^2+\dot{\theta}_2^2+2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2)\right]-T+V\nonumber\\[0.1in]&=T+V\end{align} es decir, la energía total del sistema, lo que se pudo haber deducido por simple inspección, ya que
\begin{equation}\frac{\partial(x_i,y_i)}{\partial{t}}=0,\hspace{0.25in}\frac{\partial{V}}{\partial\dot{\theta}_1}=\frac{\partial{V}}{\partial\dot{\theta}_2}=0\end{equation} además ésta es, aparentemente, la única constante que puede obtenerse, ya que ambas, $\displaystyle{\frac{\partial\mathcal{L}}{\partial\theta_1}},\;\displaystyle{\frac{\partial\mathcal{L}}{\partial\theta_2}}\neq{0}$, por esta razón entonces se recurre a obtener directamente las ecuaciones de movimiento por ecuaciones de Euler-Lagrange,
\begin{align}\frac{\partial\mathcal{L}}{\partial\theta_i}-\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\theta}_i}=0\end{align} de donde se siguen las ecuaciones de movimiento
\begin{align}\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)+2\frac{g}{\ell}\sin\theta_1&=2\ddot{\theta}_1+\ddot{\theta}_2\cos(\theta_1-\theta_2)+\dot{\theta}_2\sin(\theta_1-\theta_2)(\dot{\theta}_1-\dot{\theta}_2)\nonumber\\&\\[0.1in]\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)+\frac{g}{\ell}\sin\theta_2&=\ddot{\theta}_2+\ddot{\theta}_1\cos(\theta_1-\theta_2)+\dot{\theta}_1\sin(\theta_1-\theta_2)(\dot{\theta}_1-\dot{\theta}_2)\nonumber\\&\end{align} cuya solución puede encontrarse numéricamente con algún software como Mathematica.

Acá muestro comandos de solución en Mathematica para los valores arbitrarios
$$m=1,\,g=9.81,\;\ell=10,\;\theta_1(0)=\frac{\pi}{3},\;\theta_2(0)=0,\;\dot{\theta}_1(0)=\dot{\theta}_2(0)=-1,\;0\leq{t}\leq{100}$$

SetAttributes[{l, m, g}, Constant];

(*Define las posiciones de ambas masas*)
x1[t_] := l Sin[T1[t]]
y1[t_] := -l Cos[T1[t]]
x2[t_] := l Sin[T2[t]] + x1[t]
y2[t_] := -l Cos[T2[t]] + y1[t]

(*Define la lagrangiana del sistema*)
T = 1/2 m (x1'[t]^2 + x2'[t]^2 + y1'[t]^2 + y2'[t]^2);
V = m g (y1[t] + y2[t]);
L = T - V;

Needs["VariationalMethods`"]
(*Antes compruébese que las ecuaciones de movimiento son correctas*)

g = 9.81; l = 10; tmax = 100; m=1;

(*Resuelve numéricamente para valores arbitrarios*)
Sol = NDSolve[{EulerEquations[L, T1[t], t],
EulerEquations[L, T2[t], t], T1[0] == Pi/3,
T2[0] == 0, T1'[0] == T2'[0] == -1},
{T1[t], T2[t]}, {t, tmax}];

(*Grafica en el espacio [x,y]*)
GraphicsRow[{ParametricPlot[{x1[t], y1[t]} /. {Sol}, {t, 0, tmax},
AxesLabel -> {"x1", "y1"}],
ParametricPlot[{x2[t], y2[t]} /. {Sol}, {t, 0, tmax},
AxesLabel -> {"x2", "y2"}]}, ImageSize -> 1000]

(*Grafica en el espacio [T1,t], [T2,t]*)
GraphicsRow[{Plot[T1[t] /. {Sol}, {t, 0, tmax},
AxesLabel -> {t, T1}],
Plot[T2[t] /. {Sol}, {t, 0, tmax},
AxesLabel -> {t, T2}]}, ImageSize -> 1000]

(*Grafica el espacio de configuraciones*)
ParametricPlot[{T2[t], T1[t]} /. {Sol}, {t, 0, tmax},
AxesLabel -> {T2, T1}, ImageSize -> 1000]

(*Observa la evolución del sistema en el plano [x,y]*)
Manipulate[
ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}} /. {Sol}, {t, 0, 0 + a},
AxesLabel -> {"x", "y"}], {{a, 0.1, "Animación"}, 0, \[Infinity],
ControlType -> Trigger, PerformanceGoal -> "Quality"}]

Una de las salidas generadas en este caso, en el plano ${\{x,y\}}$ se ve así, para cada masa

Paréntesis de Poisson y el vector de Laplace-Runge-Lenz

Aprovechando que he tenido que hablar sobre el paréntesis de Poisson en un curso de Mecánica Clásica, acá comparto el cálculo del paréntesis de Poisson del vector de Laplace-Runge-Lenz (LRL) con el hamiltoniano, que -como se espera- resulta ser nulo, demostrando que el vector LRL se conserva. El camino que tomé es un tortuoso, pero funciona y ayuda a familiarizarse con las propiedades del paréntesis de Poisson. Por lo extenso de algunas ecuaciones, recomiendo leer la entrada desde un ordenador.

Sabiendo que el hamiltoniano es $\displaystyle{\mathcal{H}=\displaystyle{\frac{\mathbf{p}^2}{2\mu}-\frac{k}{r}}}$, se quiere verificar que el vector LRL, $\displaystyle{\mathbf{A}=\mathbf{p}\times\mathbf{L}-\mu{k}\displaystyle{\frac{\mathbf{r}}{r}}}$ es una constante de movimiento, calculando explícitamente el paréntesis ${\{\mathbf{A},\mathcal{H}\}}$.

Descomponiendo el vector LRL,
\begin{align}A_i&=\epsilon_{ijk}p_jL_k-\mu{k}\frac{r_i}{r}\nonumber\\[0.1in]&=\epsilon_{ijk}\epsilon_{kmn}p_jr_mp_n-\mu{k}\frac{r_i}{r}\nonumber\\[0.1in]&=\epsilon_{kij}\epsilon_{kmn}p_jr_mp_n-\mu{k}\frac{r_i}{r}\nonumber\\[0.1in]&=\left(\delta_{im}\delta_{jn}-\delta_{in}\delta_{jm}\right)p_jr_mp_n-\mu{k}\frac{r_i}{r}\nonumber\\[0.1in]&=\mathbf{p}^2r_i-(\mathbf{r}\cdot\mathbf{p})p_i-\mu{k}\frac{r_i}{r}\nonumber\\[0.1in]&=\left(\mathbf{p}^2-\mu{k}\frac{1}{r}\right)r_i-(\mathbf{r}\cdot\mathbf{p})p_i\end{align}
Así entonces, calculando el paréntesis de Poisson,
\begin{align}\{A_i,\mathcal{H}\}&=\left\{\left(\mathbf{p}^2-\mu{k}\frac{1}{r}\right)r_i-(\mathbf{r}\cdot\mathbf{p})p_i,\frac{\mathbf{p}^2}{2\mu}-\frac{k}{r}\right\}\nonumber\\[0.1in]&=\frac{1}{2\mu}\left(\left\{\mathbf{p}^2r_i,\mathbf{p}^2\right\}-\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\mathbf{p}^2\right\}\right)+k\left(\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\frac{1}{r}\right\}-\left\{\mathbf{p}^2r_i,\frac{1}{r}\right\}-\frac{1}{2}\left\{\frac{r_i}{r},\mathbf{p}^2\right\}\right)+\mu{k}^2\left(\left\{\frac{r_i}{r},\frac{1}{r}\right\}\right)\end{align} así, para el primer sumando,
\begin{align}\left\{\mathbf{p}^2r_i,\mathbf{p}^2\right\}&=\mathbf{p}^2\left\{r_i,\mathbf{p}^2\right\}+r_i\left\{\mathbf{p}^2,\mathbf{p}^2\right\}\nonumber\\[0.1in]&=\mathbf{p}^2\left\{r_i,\mathbf{p}^2\right\}\nonumber\\[0.1in]&=\mathbf{p}^2\,\sum_{\alpha}\left(\frac{\partial{r_i}}{\partial{r_\alpha}}\frac{\partial{\mathbf{p}^2}}{\partial{p_\alpha}}-\frac{\partial{r_i}}{\partial{p_\alpha}}\frac{\partial{\mathbf{p}^2}}{\partial{r_\alpha}}\right)\nonumber\\[0.1in]&=2p_i\mathbf{p}^2\end{align} y también
\begin{align}\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\mathbf{p}^2\right\}&=\left(\mathbf{r}\cdot\mathbf{p}\right)\left\{p_i,\mathbf{p}^2\right\}+p_i\left\{(\mathbf{r}\cdot\mathbf{p}),\mathbf{p}^2\right\}\nonumber\\[0.1in]&=p_i\left\{(\mathbf{r}\cdot\mathbf{p}),\mathbf{p}^2\right\}\nonumber\\[0.1in]&=p_i\,\sum_\alpha\left(\frac{\partial}{\partial{r_\alpha}}(\mathbf{r}\cdot\mathbf{p})\frac{\partial}{\partial{p_\alpha}}\mathbf{p}^2-\frac{\partial}{\partial{p_\alpha}}(\mathbf{r}\cdot\mathbf{p})\frac{\partial}{\partial{r_\alpha}}\mathbf{p}^2\right)\nonumber\\[0.1in]&=2p_i\mathbf{p}^2\end{align}
Entonces podemos ir reduciendo el vector LRL a
\begin{equation}\{A_i,\mathcal{H}\}=k\left(\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\frac{1}{r}\right\}-\left\{\mathbf{p}^2r_i,\frac{1}{r}\right\}-\frac{1}{2}\left\{\frac{r_i}{r},\mathbf{p}^2\right\}\right)+\mu{k}^2\left(\left\{\frac{r_i}{r},\frac{1}{r}\right\}\right)\end{equation}
Para el siguiente sumando,
\begin{align}\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\frac{1}{r}\right\}&=(\mathbf{r}\cdot\mathbf{p})\left\{p_i,\frac{1}{r}\right\}+p_i\left\{(\mathbf{r}\cdot\mathbf{p}),\frac{1}{r}\right\}\nonumber\\[0.1in]&=(\mathbf{r}\cdot\mathbf{p})\,\sum_\alpha\left(\frac{\partial{p_i}}{\partial{r_\alpha}}\frac{\partial{r^{-1}}}{\partial{p_\alpha}}-\frac{\partial{p_i}}{\partial{p_\alpha}}\frac{\partial{r^{-1}}}{\partial{r_\alpha}}\right)+p_i\,\sum_\alpha\left(\frac{\partial}{\partial{r_\alpha}}(\mathbf{r}\cdot\mathbf{p})\frac{\partial{r^{-1}}}{\partial{p_\alpha}}-\frac{\partial}{\partial{p_\alpha}}(\mathbf{r}\cdot\mathbf{p})\frac{\partial{r^{-1}}}{\partial{r_\alpha}}\right)\nonumber\\[0.1in]&=(\mathbf{r}\cdot\mathbf{p})\frac{r_i}{r^3}+\frac{p_i}{r}\nonumber\\[0.1in]&=\left\{\mathbf{p}^2r_i,\frac{1}{r}\right\}\nonumber\\[0.1in]&=\mathbf{p}^2\left\{r_i,\frac{1}{r}\right\}+r_i\left\{\mathbf{p}^2,\frac{1}{r}\right\}\nonumber\\[0.1in]&=\mathbf{p}^2\,\sum_\alpha\left(\frac{\partial{r_i}}{\partial{r_\alpha}}\frac{\partial{r^{-1}}}{\partial{p_\alpha}}-\frac{\partial{r_i}}{\partial{p_\alpha}}\frac{\partial{r^{-1}}}{\partial{r_\alpha}}\right)+r_i\,\sum_\alpha\left(\frac{\partial}{\partial{r_\alpha}}\mathbf{p}^2\frac{\partial{r^{-1}}}{\partial{p_\alpha}}-\frac{\partial}{\partial{p_\alpha}}\mathbf{p}^2\frac{\partial{r^{-1}}}{\partial{r_\alpha}}\right)\nonumber\\[0.1in]&=2\frac{r_i}{r^3}(\mathbf{r}\cdot\mathbf{p})\end{align} y también
\begin{align}\left\{\frac{r_i}{r},\mathbf{p}^2\right\}&=r_i\left\{\frac{1}{r},\mathbf{p}^2\right\}+\frac{1}{r}\left\{r_i,\mathbf{p}^2\right\}\nonumber\\[0.1in]&=-2r_i\frac{1}{r^3}(\mathbf{r}\cdot\mathbf{p})+\frac{2}{r}p_i\end{align} entonces se tiene
\begin{align}k&\left(\left\{(\mathbf{r}\cdot\mathbf{p})p_i,\frac{1}{r}\right\}-\left\{\mathbf{p}^2r_i,\frac{1}{r}\right\}-\frac{1}{2}\left\{\frac{r_i}{r},\mathbf{p}^2\right\}\right)=k\left[\frac{p_i}{r}+(\mathbf{r}\cdot\mathbf{p})\frac{r_i}{r^3}-2(\mathbf{r}\cdot\mathbf{p})\frac{r_i}{r^3}-\frac{1}{2}\left(\frac{2}{r}p_i-2(\mathbf{r}\cdot\mathbf{p})\frac{r_i}{r^3}\right)\right]=0\end{align} por lo que la derivada temporal del vector LRL se reduce a
\begin{equation}\{A_i,\mathcal{H}\}=\mu{k}^2\left\{\frac{r_i}{r},\frac{1}{r}\right\}\end{equation} que evidentemente es nulo, por tanto se concluye que
\begin{equation}\mathbf{\dot{A}}=\left\{\mathbf{A},\mathcal{H}\right\}=\mathbf{0}\end{equation} y en efecto $\mathbf{A}$ es constante de movimiento.

En Mathematica es sencillo programar el paréntesis de Poisson para verificar lo anterior. Las siguientes líneas hacen esa tarea para cada componente.

(*Definición paréntesis de Poisson*)
PB[u_, v_, q_Symbol, p_Symbol] := D[u, q] D[v, p] - D[v, q] D[u, p]

(*Definiciones*)
p = {p1, p2, p3}; r = {r1, r2, r3};
P = FullSimplify[Norm[p], p1 > 0 && p2 > 0 && p3 > 0]
R = FullSimplify[Norm[r], r1 > 0 && r2 > 0 && r3 > 0]
A := Cross[p, Cross[r, p]] - r/R
H = P^2/2 - 1/R

(*El vector LRL es ortogonal al momento angular*)
A . Cross[r, p] == 0 // Simplify

(*El vector LRL se conserva*)
Simplify[Sum[{PB[A[[1]], H, r[[j]], p[[j]]], PB[A[[2]], H, r[[j]],
p[[j]]], PB[A[[3]], H, r[[j]], p[[j]]]}, {j, 1, 3}]]

Strong Mathematical Induction and Fibonacci Numbers

I had never seen the second principle of induction before, also called strong mathematical induction or complete induction, but the name suggests really all that it is: a wider variant of mathematical induction. When we use the principle of finite induction (mathematical induction), we know that a proposition ${P(n),\,\forall{n}\in\mathbb{N}}$, is true whenever there is some ${n=k}$ such that ${P(k)}$ and ${P(k+1)}$ holds and the proposition is true for the base case (usually 1). To see this is rather easy and intuitive, since it is just a generalization for all natural numbers of some particular proposition. But sometimes the assumption that ${P(k)}$ holds for ${n=k}$ is not enough and you’ve got to assume that the proposition holds for all ${n\geq{k}}$; that is, for all the elements of a subset ${A=\{n_0\in\mathbb{N}\,|\,n\geq{n_0,n_0+1,\ldots,k-1,k}\}}$ where ${n_0}$ is the base case. It often just entails proving the base case for each element.The following example, for Fibonacci numbers, illustrates the proof of a proposition by means of strong mathematical induction. The original sentence can be found in “Elements of the theory of numbers” by Joseph & Thomas Dence.

The Fibonacci numbers denoted ${F_n}$ are defined recursively by ${F_1=F_2=1}$, ${F_n=F_{n-1}+F_{n-2}}$ for ${n>2}$. Show that for ${n\in\mathbb{N}}$, ${F_{n+1}\leq[(1+\sqrt{5})/2]^n}$.

It can be easily verified for the base case ${n=1}$ to begin the proof.

Let us assume that for ${k\leq{n}}$, the proposition ${P(k):\,F_{k+1}\leq\phi^k}$, where $\phi\equiv(1+\sqrt{5})/2$ (the golden ratio), is true; then for ${k+1}$, *
$$F_{k+2}=F_{k+1}+F_{k}\leq\phi^{k-1}(\phi+1)$$ now notice that
$$\phi+1=\frac{3+\sqrt{5}}{2}=\frac{6+2\sqrt{5}}{4}=\frac{1+2\sqrt{5}+5}{4}=\phi^2$$ so the previous inequality becomes evident as
$$F_{k+2}\leq\phi^{k+1}$$ and so ${P(k+1)}$ is also true. Hence we say that by induction ${P(n)}$ is true ${\forall{n}\in\mathbb{N}}$.

This illustrates the situation in Mathematica for continuous functions

with the command

<<PlotLegends`

Plot[{Fibonacci[x+1], GoldenRatio^x}, {x,1,20},
Filling -> {1->{2}}, PlotStyle -> Thick,
PlotLegend -> {"F_n+1", "phi^n"},
LegendPosition -> {-0.75,0}]

Notice the exponential growth of the Fibonacci numbers. Now dare you to prove that ${F_{n+1}\geq\phi^{n-1},\,\forall{n}\geq{1}}$.

* This assumption is the key of the strong mathematical induction. In this problem the cases ${n=1,\ldots,k-1,k}$ are considered. Notice that for the base case ${n=1}$, each proposition holds.