Showing posts with label Programación. Show all posts
Showing posts with label Programación. Show all posts

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}]
Curvas de nivel a t fijo - radiación dipolo eléctrico (dipole radiation)

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]
Animación curvas de nivel a y fijo - radiación dipolo eléctrico (dipole radiation animation GIF)
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.

Mujeres en la ciencia

Me ha surgido inquietud por este tema del rol actual de la mujer en la ciencia (si se quiere también en cuestiones técnicas) por la siguiente fotografía

Imagen (recarga la página)

Se trata de la foto grupal de la EMMN 2013, acerca de la cual hablo en la entrada anterior. Aunque falta un grupo de 3 o 4 chicas por ahí, salta inmediatamente a la vista que la mayoría de asistentes somos hombres.

En la UAM-I por ejemplo, el número total de alumnos (sin importar género) en cada división de Ciencias Básicas e Ingeniería (CBI) es prácticamente la mitad del número de alumnos en la licenciatura en Administración. Además de eso, en CBI, a única excepción de la licenciatura en Química, todas las carreras tienen una mayor cantidad de hombres que mujeres, en mayor o menor grado.

Aún más, en la Olimpiada Internacional de Matemáticas (IMO) el mejor registro ha sido 59 mujeres contra 506 hombres en Alemania 2009. México ha llevado mujeres 4 veces, donde de 6 participantes, 3 veces una ha sido mujer y una vez dos han sido mujeres  (sigue acá el facebook de la Olimpiada Mexicana de Matemáticas). Aparentemente (no cuento con los datos precisos) situaciones similares se dan en la olimpiada de física, por ejemplo.

A pesar de estos datos de etapas tempranas, la situación histórica en general ha mejorado muchísimo y se han desechado bastantes prejuicios o nociones arcaicas que impedían una mayor participación, mejor desempeño y resultados de la mujer en cuestiones científicas o técnicas en general. Qué tanto se involucran los jóvenes en general en la ciencia, como tuve alguna vez oportunidad de comprobar en carne propia, es casi seguramente una situación puramente cultural, en la que al parecer, primero no se discrimina género y en segunda se enfoca a la mujer, que debe de satisfacer ciertos roles frente a los parámetros culturales establecidos, ya sea consciente o inconscientemente.

Imagen (recarga la página)

Como dije, no creo que sea cuestión ya de que se haga más o menos a la mujer en cuanto al rol que puede tener en la ciencia, si no más bien de patrones culturales bien arraigados. La imagen de arriba es de un número reciente de Nature llamado 'Women in Science', que en realidad me desagrada un poco por lo mismo que ya he dicho. Para quien desee ahondar más, seguramente abunda la información, y algo tendrán que decir los buenos amigos de las 'ciencias sociales' al respecto.

Acá comparto un artículo de dicho número, del cual prefiero reservar mi opinión, así como el lector tendrá la propia: Weird sisters?.
Como los hombres, las mujeres científicas tienen personalidades individuales e idiosincrasias, y ambos tienen debilidades así como capacidades extraordinarias –no porque sean mujeres, sino porque son seres humanos.
Finalmente un divertido cómic de Abstruse Goose en el que casualmente aparece una chica (el género no tiene relación alguna, quizá solo quería compartir el cómic: abstrusegoose.com/508) y resulta una buena motivación ;)

Imagen (recarga la página)

Caos Cuántico y la Escuela de Modelación y Métodos Numéricos 2013

Imagen (recarga la página)
Del 25 al 28 de junio del presente año se llevó acabo la Escuela de Modelación y Métodos Numéricos 2013 en el Centro de Investigación en Matemáticas (CIMAT) en Guanajuato, Guanajuato. La experiencia fue sumamente enriquecedora, de manera profesional y personal. La escuela trató los temas de Dinámica Molecular y Química Cuántica, Nanociencia y Nanotecnología. Me pareció curioso que un buen porcentaje de los asistentes al evento se dedican o a la biología o a la química; como físicos sólo asistimos dos compañeros y yo por parte de la misma institución, aunque hubo también varios matemáticos.

Acá comparto una presentación bastante breve que hicimos mis dos compañeros y yo en un minisimposio de mecánica cuántica para compartir un poco de lo tratado en la escuela y nuestra experiencia: Presentación EMMN-13

Imagen (recarga la página)

Ha sido el primer evento externo al que asisto, y vaya que fue divertido. Y para cerrar con broche de oro tanta química cuántica, un buen harlem shake molecular:

Bueno, y aprovechando la oportunidad para hablar en el minisimposio, decidí compartir un poco de información acerca de lo que es el llamado caos cuántico con esta presentación: ¿Qué es el caos cuántico?.
Imagen (recarga la página)
La plática se llevó acabo en unos 20 minutos, y lo que procuré fue sobre todo expresar mi interés por el tema y contagiar a algunos de mis compañeros. En este momento me encuentro con una gama de posibilidades enfrente sobre la dirección que puede tomar mi carrera, y ésta es una que me parece bastante atractiva, desde los sistemas dinámicos y la física no lineal, hasta el caos cuántico y aplicaciones biológicas. El tema es muy amplio y abunda información en la red, para quien desee ahondar en el tema, aunque quizá para ello antes se necesite afilar un poco el propio colmillo matemático.

Condicionales en Mathematica

Imagen (si ves este texto, recarga la página)

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,

Imagen (recarga la página)

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:

Imagen (si ves este mensaje, recarga la página)

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
Imagen (si ves este mensaje, recarga la página)

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}]]

Factoriones

Hace tiempo que no entro a resolver problemas en Project Euler. Si te gusta resolver acertijos matemáticos, que seguido involucren habilidades de programación, es bastante recomendable. El punto es que recordé un problema que implícitamente pide encontrar factoriones.Los pide implícitamente pues sólo existen cuatro factoriones en base 10. Un factorion es un número que es igual a la suma de factoriales de sus dígitos, e.g. 145=1!+4!+5!

Para hallarlos simplemente hay que programar directamente el algoritmo para discriminar si un número es o no factorion. El verdadero problema es que el programa no pase días buscando factoriones. Para lograr esto, el procedimiento estándar es notar que k(9!) tendrá menos de k dígitos siempre que k>7, por tanto sólo hay que buscar factoriones hasta 7(9!)=2540160

Ahora bien, esto aún sigue sin ser muy eficiente. Esto es simplemente un puzzle de programación, como dije sólo existen cuatro factoriones, se conocen y punto.

El factorion más grande es 40585, así que probablemente exista una forma de recortar el límite superior a partir de un argumento lógico. Mi lenguaje ya predilecto es Python, y no sé si sea precisamente si es él el que hace lento el programa.

Acá comparto el código:
# 145 es un numero curioso, ya que 1! + 4! + 5! = 1 + 24 + 120 = 145.
# Encuentra la suma de todos los numeros iguales a la suma
# del factorial de sus digitos.
# Nota: como 1! = 1 y 2! = 2 no son sumas, no se incluyen.

import time
t0 = time.time()

def factorial(n):
if n > 1: return (n*factorial(n-1))
else: return 1

for i in range(1, 7*factorial(9)): # LIMITE SUPERIOR ÓPTIMO~1*(10^5)
n, s = str(i), 0
for j in range(len(n)):
s += factorial(int(n[j]))
if s == i: print(i)
print('Tiempo de ejecucion: ',time.time()-t0,'segundos.')

Tutorial LaTeX: ¿por dónde empezar?

Imagen (si ves este texto recarga la pag)$\LaTeX$ es prácticamente indispensable para escribir documentos científicos, pero aún más, es plausible que cualquier persona que genere documentos lo haga con LaTeX, desde un Curriculum Vitae, hasta un libro o publicaciones oficiales. Quizá lo más complicado -como en cualquier cosa- al querer comenzar a usar LaTeX para generar documentos -desde escolares hasta profesionales- es precisamente dar un paso adelante. Yo lo hice por mi mismo y debo aceptar que es algo embrollado al principio el cómo si quiera instalar LaTeX; uno busca "Latex" en google y no sabe ni por dónde.

LaTeX (pronunciado algo así como Lei-Tec o La-TeJ) es un sistema para composición de textos (basado en el lenguaje TeX), distinto por ejemplo a Microsoft Word, en que no es un sistema de WYSIWYG (what you see is what you get, lo que ves es lo que obtienes); con $\LaTeX$ compilas tu documento y de ahí generas una salida; LaTeX no es un procesador de texto común. Para esto no es necesario saber programación, hacerlo es muy sencillo, o bien hoy existen interfaces gráficas que lo pueden hacer por ti.

Algunas razones que considero para recomendar usar $\LaTeX$ son:
  • Calidad y Presentación: Un documento generado en LaTeX tiene un aspecto mucho más profesional y atractivo, independientemente de si el contenido es técnico o no. Si generas documentos que incluyen matemáticas, LaTeX es prácticamente una necesidad: un documento hecho en un procesador de textos es una pésima opción.
  • Costo: En general, LaTeX es software libre (incluyendo el no tener costo monetario, cosa que no siempre implica el software libre) y no son necesarias licencias extras para su uso y distribución.
  • Compatibilidad: Un documento .tex es simplemente un documento de texto que puedes abrir en cualquier editor de texto como el bloc de notas, sin importar el sistema operativo. La salida de un documento .tex por defecto es el DVI, que luego se convierte a al ya familiar pdf, que tiene la ventaja de no sufrir alteraciones de formato al cambiar de una plataforma a otra.
  • Contenido y Estilo: Con LaTeX simplemente escribes, separando el contenido del trabajo del estilo que le das; contrario a los procesadores de texto, en los cuales editas todo gráficamente.

Considera el trabajo sobre el Warp Drive de Miguel Alcubierre como ejemplo de una salida de LaTeX (si sabes sobre relatividad general, pues mucho mejor). En general los científicos utilizan LaTeX, y en arXiv se puede comprobar.

Lo primero que hay que hacer es descargarse MikTeX, aquí para Windows o acá (MacTeX) para Mac. En Linux yo utilizo TeX live. Estos son los motores (distribuciones) necesarios para generar documentos con LaTeX. Ya hecho esto, yo recomendaría instalar una interfaz gráfica, como TeXnicCenter, TeXworks o TeXmaker, después si así quieres podrás utilizar solamente un editor de texto y una ventana de comandos; yo me guiaré con TeXworks en adelante.

Hay muchos templates o plantillas para crear documentos más elaborados con un estilo predefinido. Por ahora describiré lo básico para generar un documento y así familiarizarse con LaTeX.

Abre TeXworks y abre un nuevo documento con el ícono de hoja blanca. Escribe las líneas
\documentclass{article}
\begin{document}
Hola mundo
\end{document}
Los signos { } se utilizan, en general para agrupar contenido en LaTeX. A un lado del botón verde con el símbolo de Play, escoge PdfLatex (para generar un documento pdf) y da clic al botón verde (o Play). Se te pedirá guardar el documento .tex (éste será el código fuente, o documento fuente), escoge una carpeta (se incluirán varios documentos al momento de compilar) y guarda con el nombre que gustes, digamos "tutorial.tex". Cuando hay errores, TeXworks te los mostrará y no generará el pdf. En los otros programas probablemente te encuentres con las palabras "Build" o "Compile" para compilar y poder ver el documento pdf. El documento pdf (y otros documentos de compilación) se generará en la carpeta que creaste para el .tex.

Bueno, de aquí sólo queda ir descubriendo poco a poco el mundo de LaTeX. Es importante mencionar las líneas \documentclass[ ]{ }, que define la clase de documento (artículo, reporte, libro, etc...) y \usepackage[ ]{ } (que será necesaria más adelante), que define los llamados paquetes, que luego permiten ciertos comportamientos predefinidos para los documentos. Con la práctica incluso podrás crear luego los paquetes que se adapten a tus necesidades e incluirlos en un solo comando. Los corchetes [ ] se emplean para señalar características específicas, ya lo irás notando.

En fin, ésta es la parte difícil, ahora sólo es cuestión tuya. En Wikibooks hay un manual de LaTeX; muy bueno para cada pequeño detalle con que te vayas topando (en inglés está mucho más completo). En la cuestión de las matemáticas, es muy sencillo escribirlas. En la red abundan listas de símbolos matemáticos como ésta  y lo único que hay que hacer es escribirlos ya sea entre signos \$...\$ para ecuaciones dentro de párrafos o \$\$...\$\$ para ecuaciones fuera de párrafos. Hay opciones más elaboradas para escribir en este llamado mathmode y opciones más avanzadas para escribir matemáticas, como escribir matrices, arreglos, etc... Confío en que la habilidad del lector para escribir situaciones matemáticas más elaboradas llegará con la práctica misma.

El siguiente código puede servirte como plantilla para que inicies con algo más avanzado. Los comandos -como siempre- están hechos para sugerir su función, por lo que es muy sencillo deducir la mayoría de ellos. Intenta primero cambiando la clase de documento de report a article, o book. De ahí, juega con el código, sólo asegurate de descargar los paquetes incluidos (en Windows y supongo que en Mac también, MikTeX hará todo por ti).
\documentclass[11pt, a4paper]{report}
\usepackage[spanish,mexico]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amsmath} %Paquetes de la American Mathematical Society
\usepackage{amssymb}
\usepackage{hyperref} %Enlaces
\usepackage[nottoc]{tocbibind} %Bibliografía y nottoc evita que el índice se
% incluya a sí mismo

\newcommand*{\f}{\frac} % Cuando se quiera escribir una fracción,
%sólo será necesario escribir \f{}{} y no \frac{}{}

\begin{document}

\title{\textbf{El Teorema Integral de Cauchy}}
\author{Pedro Figueroa
\footnote{\href{mailto:pedrofigueroa@live.com.mx}
{\texttt{pedrofigueroa@live.com.mx}}}\\
Universidad Autónoma Metropolitana
\vspace{0.5in}
}

\date{\today}

\maketitle
\newpage
\setcounter{tocdepth}{1}
\pagenumbering{gobble} % Evita numeración de pgs
\tableofcontents % Índice

\newpage

\part{El Teorema}
%\section* hace que no se numere ni muestre en el índice la sección
\section{Teorema Integral de Cauchy}
El teorema integral de Cauchy se lee: Si $f(z)$ es analítica
en un dominio simplemente conexo $\mathcal{D}$,
para todo contorno $C$ en $\mathcal{D}$ se cumple
$$\oint_Cf(z)\,dz=0$$
Este teorema es de lo más conmovedor si además consideramos su demostración
para el caso $f^\prime(z)$ continua.

\vspace{0.25 in} %Nota el comportamiento de los espacios

Sabemos que \(f(z)=u(x,y)+iv(x,y)\). Si para la integral de línea se tiene
$$S_n=\sum_{m=1}^n(u+iv)(\Delta{x}_m+i\Delta{y}_m)$$
\begin{align*}
\lim_{n\to\infty}S_n
&amp;=\int_{\tilde{C}}f(z)\,dz\\
&amp;=\int_{\tilde{C}}\left(u\,dx-v\,dy\right)+
i\int_\tilde{C}\left(u\,dy+v\,dx\right)\end{align*}
% align, equation, multline son formas para
% escribir matemáticas fuera de párrafos,
% * evita que las ecuaciones se numeren
entonces análogamente, para un contorno \(C\)
$$\oint_Cf(z)\,dz=\oint_{C}\left(u\,dx-v\,dy\right)+
i\,\oint_{C}\left(u\,dy+v\,dx\right)$$
Como consideramos el caso en que ${f^\prime(z)}$ es continua,
$u$ y $v$ tienen derivadas parciales continuas en $\mathcal{D}$.

Lo conmovedor de esta demostración, es que es aplicable el teorema de Green,
que nos dice que si ${u,\,v}$ tienen derivadas parciales continuas en una región
abierta $R$ en $\mathcal{D}$, se cumple
$$\int_{\tilde{C}}(u\,dx-v\,dy)
=\iint\limits_R\left(-\f{\partial{v}}{\partial{x}}-
\f{\partial{u}}{\partial{y}}\right)\;dxdy$$
\ldots{}y como la gente sabe, para toda función analítica en $\mathcal{D}$
se deben cumplir las ecuaciones Cauchy- Riemann, a saber:
$$\f{\partial{u}}{\partial{x}}
=\f{\partial{v}}{\partial{y}}\hspace{0.5in}\f{\partial{u}}{\partial{y}}
=-\f{\partial{v}}{\partial{x}}$$
Y con ello se ha demostrado el teorema integral de Cauchy
para $f^\prime(z)$ continua.

\vspace{0.25 in}

Además, una de las consecuencias más importantes de este teorema es la llamada
fórmula integral de Cauchy, que dice que
$$f(z_0)=\f{1}{2\pi{i}}\,\oint_C\,\f{f(z)}{z-z_0}\,dz$$
con $C$ un contorno que encierra al punto $z_0$.

\section{Forma diferencial}
O bien, se tiene también la llamada forma diferencial de la fórmula,

$$f^{(n)}(z_0)=\f{n!}{2\pi{i}}\,\oint_C\,\f{f(z)}{(z-z_0)^{n+1}}\,dz,
\hspace{0.1in}\forall{n}\in\mathbb{N}$$

Édouard Goursat, un matemático francés, por el año 1900 hizo la demostración
del teorema integral de Cauchy sin la condición $f^\prime{(z)}$ continua,
que además contribuye, por ejemplo, al hecho de que la derivada de una función
analítica también sea analítica.

\vfill
\renewcommand{\bibname}{Referencias}
\begin{thebibliography}{2}
\bibitem{Ahlfors}
Ahlfors,
\emph{Complex Analysis}.
Mc Graw Hill,
3a Edition,
1979.
\end{thebibliography}

\end{document}

Comandos personalizados en LaTeX

Hace algunos días tuve que lidiar con demostraciones de identidades que involucran al operador ${\nabla}$ y por tanto derivadas parciales. Tal vez exista forma de demostrar dichas identidades sin hacer operaciones explícitamente; intenté generalizar a partir de identidades de simples vectores tomando en cuenta que $\nabla$ es un vector de operadores, sin embargo no lograban demostrar rigurosamente las identidades.

Una identidad a demostrar, por ejemplo, es:
$$\mathbf{\nabla\times\left(A\times{B}\right)=\left(B\cdot\nabla\right)A-B\left(\nabla\cdot{A}\right)-\left(A\cdot\nabla\right)B+A\left(\nabla\cdot{B}\right)}$$ con ${\mathbf{A,\,B}}$ dos funciones ${\mathbb{R}^3\rightarrow\mathbb{R}^3}$ bien comportadas cualesquiera. Bueno, ahora el problema no era demostrar las identidades, sino escribir en LaTeX tremendos engendros. Claro, pude haber elegido simplemente utilizar papel y lápiz como cualquier mortal, pero ya había comenzado antes y ahora era muy tarde, así que debía seguir con LaTeX.

Bueno, me han pasado el tip para crear comandos propios. Digamos que a cada paso tienes que escribir cosas como
\left(\frac{\partial{A_x}}{\partial{x}}-\frac{\partial{B_y}}{\partial{y}}\right)
o algo por el estilo. Bueno, te puedes ahorrar todo aquello haciendo
\newcommand*{\p}{\partial}
\newcommand*{\f}{\frac}
\newcommand*{\l}{\left}
\newcommand*{\r}{\right}
y así, aquello horrible se convierte en
\l(\f{\p{A_x}}{\p{x}}-\f{\p{B_y}}{\p{y}}\r)
Puedes obviamente incluir lo que tu desees en tus nuevos comandos, no sólo los utilizados en math mode. De manera similar puedes crear paquetes propios y demás, he aquí un enlace a wikibooks para el interesado.

De paso aprovecho para recomendar a aquellos que aún no se deciden a aprender LaTeX a probar TeXmaker; en lo personal uso TeXworks con MikTeX y ya me he acostumbrado, sin embargo en Linux descubrí el TeXmaker y es bastante cómodo el ambiente para cualquiera que se inicie con LaTeX, además de que en trabajos muy amplios ayuda mucho a diferenciar comandos o fórmulas del contenido en sí.

Química con LaTeX

Hace algún tiempo quise escribir un trabajo escolar de fisicoquímica en LaTeX, y me di cuenta que no era tanto inconveniente el abrir un editor y comenzar a escribir en LaTeX, sino el tener que escribir símbolos químicos, reacciones, o bien equilibrios.

Lo natural para, por ejemplo, un ion sulfato, sería escribir:
$\mathrm{SO_4^{-2}}$
con lo que se obtiene ${\textrm{SO}_4^{-2}}$, sin embargo la forma más práctica de escribir química usando LaTeX, tal vez sea utilizando el paquete mhchem, sobre todo para situaciones más elaboradas, por ejemplo:
que se obtiene con
$\ce{^{239}_{92}U}$
$\ce{A ->[\ce{H2O}] B}$
$\ce{SO4_{(aq)}^2- + Cu_{(aq)}^2+ -> CuSO4 v}$
$K=\frac{[\ce{CH3CO2-}][\ce{H3O+}]}{[\ce{CH3CO2H}]}$
$\ce{HbH+ + O2 <=> HbO2 + H+}$
Tal vez el único detalle de este paquete es el generar estructuras de Lewis, por que aunque es sencillo generar enlaces:
$\ce{A\bond{-}B\bond{=}C\bond{#}D}$
$\ce{A\bond{~}B\bond{~-}C\bond{~=}D}$
$\ce{A\bond{...}B\bond{...}C\bond{....}D}$
no sé de qué forma se puedan generar dichas estructuras; de cualquier modo se puede recurrir al paquete lewis y obtener símbolos del tipo
 
escribiendo
$\lewis{A}{1}{2}{3}{4}{5}{6}{7}{8}$
por ejemplo
$\lewis{O}{.}{.}{.}{}{}{}{}{.}\ce{\bond{=}C\bond{=}}\lewis{O}{}{}{}{.}{.}{.}{.}{}$

Sobre Ubuntu y la función atoi en C++

Al fin me he cambiado por completo a Ubuntu y todo ha ido bastante bien; tal vez lo más latoso ha sido hacer que un maldito multifuncional Epson trabaje correctamente. Lo que más me preocupaba era instalar software que no se encontraba en el centro de software o en extensión .deb (ahora sé que también el .rpm se convierte fácilmente a .deb con alien). A decir verdad era por que no podía instalar Mathematica -duh!- que quizá es el único software que no he podido remplazar por alguno libre. Ahora sólo me incomoda que al parecer no existe MiKTeX para Linux que no sea i386 y tengo que instalar paquetes de LaTeX manualmente. De cualquier modo recomiendo mucho cambiarse y abandonar Windows o Mac, pues aunque es seguro equivocarse o fallar (como por ejemplo borrar el maldito cargador de arranque, no poder ejecutar un maldito archivo o ni siquiera saber acceder a un maldito directorio), también lo es aprender bastante cómo funciona el maldito corazón de tu computadora.Y bueno, quisiera dejar de maldecir, pero acá me he encontrado otro maldito problema! Obviamente he conservado mis bonitos programas en C++, aunque desde que descubrí Python no suelo usar ni C++ ni Java. Bueno, el problema es que con Windows me había mal acostumbrado a Visual Studio y no tenía ni idea de cómo compilar/ejecutar desde la terminal; pues solucioné eso rápidamente y acomodé varios programas de Project Euler, hasta que me encontré con el código del problema 8. El desgraciado no daba el resultado, ¡sólo arrojaba un maldito 0! el problema: la función atoi( ).

En mi forma de resolver el problema, debo extraer un caracter de la enorme cadena y convertir a entero para poder operar con ella. Al parecer en Visual Studio podía extraer un sólo caracter de una cadena y convertirlo a entero con atoi, mientras en g++ (el compilador predeterminado en Ubuntu, o gcc) no. Por ejemplo
char a[]="123";
cout << atoi(&a[1]);
en g++ no imprimiría un maldito "2", mientras que en Visual Studio sí. Tal vez me equivoque pero ya estará el buen usuario de Visual Studio para desmentirme, pues si no es así, no sé cómo demonios di la respuesta del problema.

Pero bueno, "resolverlo" es sencillo e incluso mejora el maldito programa. En Ubuntu la salida de lo anterior sería "23", asi que usar atoi( ) me parece descartado. Lo mejor sería hacer algo como
char a[]="123";
cout << int(a[1])-'0';
lo que suele llamarse "type-casting" o "conversión de tipos" pues lo que hace es dar información sobre el número de caracter ASCII de a[1]; por ello es necesario restar el ASCII de '0'.

Regresando al tema de Ubuntu, es incluso un tema ya un tanto desgastado, pero aunque en general es un poco más complicado cambiarse de Windows a Linux sin guía alguna más que Google; si no requieres hacer computación más fuerte, Ubuntu satisface todo lo que Windows y tienes la posibilidad de aprender un poco más, mientras que si te gusta la computación, incluso aunque no seas desarrollador, resulta bastante fructífero el intentarlo, además de todos los beneficios que acarrea el software libre en la sociedad. Consulta la filosofía de Ubuntu.

Me gustaría saber si existe algún lugar donde conseguir ordenadores sin OS o con algún distro de Linux (no Mac) preinstalado. Es un poco frustrante querer conseguir un ordenador nuevo y verse obligado a pagar la licencia de Windows o Mac.

"How to switch from Windows to Linux"

La espiral de Ulam y project euler 58

Llevo dos días sin poder resolver el problema 51 de Project Euler, así que decidí intentar algún otro. Encontré el que muestro a continuación, que me pareció divertido y sencillo, el único detalle es que mi solución no es tan eficiente como pensé.

El motivo de esta entrada es desafiarte a resolverlo mejorando el análisis que yo hice (es decir, el tiempo que tardó mi programa si es que programarás en Python); no importa cómo lo resuelvas, lo importante es que la solución sea elegante, sencilla y rápida.

A continuación muestro el problema, mi resultado y tanto el análisis como el programa que hice para obtenerlo. Mi recomendación es que primero lo resuelvas a tu manera. Si te interesa resolver más problemas entra a project euler y lleva tu propio récord.

Comenzando con 1 y rotando en sentido contrario a las agujas del reloj de la siguiente manera, se forma una espiral cuadrada con lados de longitud 7.

37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49


Es interesante que los cuadrados impares se sitúan en la diagonal inferior derecha, pero aún más interesante es que 8 de los 13 números de las dos diagonales son primos, es decir, una proporción de 8/13 ≈ 62%.

Si enrollamos una nueva capa en esta espiral, se formará una espiral cuadrada con lados de longitud 9. Si este proceso continúa, ¿cuál es la longitud del lado de la espiral cuadrada para la cual la relación de primos en diagonales es inferior al 10%?

Mi respuesta:
Side length of the square spiral: 26241
Time elapsed: 26.3150000572 seconds

Análisis y código fuente:
Comenzando desde 1, el número de la $i$-ésima esquina superior derecha de la espiral cuadrada, $k_i$ está dado por
$$k_i=k_{i-1}+8(i-1)+2\hspace{0.25in}\text{con}\hspace{0.25in}k_0=1$$ Para las 3 esquinas restantes (que etiqueto con 1, 2 y 3), en cada vuelta
\begin{align*}k_{(i,1)}&=k_i+2i\\
k_{(i,2)}&=k_i+4i\\
k_{(i,3)}&=k_i+6i\end{align*} Así entonces, generamos
\begin{align*}k_1&=k_0+2=1+2=3\\
k_{(1,1)}&=k_1+2(1)=3+2=5\\
k_{(1,2)}&=k_1+4(1)=3+4=7\\
k_{(1,3)}&=k_1+6(1)=3+6=9\end{align*} Falta revisar si los números son primos y contar. En el problema se menciona que en la diagonal inferior izquierda todos son cuadrados impares, así que éstos se pueden ignorar.

Mi código en Python se ve así:
from __future__ import division
import time

t0 = time.time()

def isprime(x, i=2):
while i <= (x/i):
if x%i == 0: return False
i += 1
return True

prime, total = 0, 1
k, n = 3, 2
while prime/total >= 0.1 or prime == 0:
for i in range(3):
if isprime(k+i*n) == True:
prime += 1
k = k + 4*n + 2
total += 4
n += 2

print 'Side length of the square spiral:',n-1
print 'Time elapsed:',time.time()-t0,'seconds'

Apenas me he enterado que la espiral utilizada se llama Espiral de Ulam, para el que esté interesado en conocer más. La espiral en general, por ejemplo, muestra la siguiente distribución de números primos:

Números de Lychrel, Algoritmo196 y Project Euler 55

Uno de los problemas abiertos de Teoría de Números es encontrar números de Lychrel. Estos números hipotéticos desafían el llamado algoritmo196, -una secuencia que utiliza el viejo truco de sumar y agregar para obtener números palíndromos o capicúas. El algoritmo196 toma cualquier número natural, invierte sus dígitos y le suma el original; a nuestro interés, si la suma no es capicúa repite el proceso a partir de dicha suma, por ejemplo

338 + 833 = 1171
1171 + 1711 = 2882

Los números de Lychrel son entonces, los que no producen capicúas al aplicar el algoritmo196. De nuevo, esto es un problema abierto, y precisamente 196 es el primer número de Lychrel que se ha conjeturado. Consulta La conjetura del 196 en Gaussianos. En seguida muestro una parte del problema 55 de Project Euler y mi solución, que es vulgar y silvestre lógica potenciada con el lenguaje Python, en realidad no hice un análisis más allá.

[…] Aunque nadie lo ha demostrado aún, se cree que algunos números, como el 196, nunca producen un capicúa. Un número que nunca forma un capicúa mediante el proceso de sumarse a sí mismo invertido se denomina un número de Lychrel. Debido a la naturaleza teórica de estos números, y para el propósito de este problema, vamos a suponer que un número es Lychrel hasta que se demuestre lo contrario. Además, se te da la información de que cada número inferior a diez mil, (i) será un capicúa en menos de 50 iteraciones, o, (ii) nadie, con toda la potencia de cálculo que existe, ha logrado por ahora asociarlo a un capicúa.

De hecho, 10677 es el primer número que requiere más de cincuenta iteraciones antes de producir el capicúa: 4668731596684224866951378664 (53 iteraciones, 28 dígitos).

Sorprendentemente, existen números capicúa que son a su vez números de Lychrel; el primer ejemplo es 4994.

¿Cuántos números de Lychrel hay inferiores a diez mil?
Solución en Python:
import time

t0 = time.time()
L, M = 0, 10000
for i in range(1,M):
s = str(i+int(str(i)[::-1]))
lychrel = True
j = 1
while j &amp;lt;= 50:
if s == s[::-1]:
lychrel = False
break
j += 1
s = str(int(s)+int(s[::-1]))
if lychrel == True:
L += 1
print 'Hay',L,'números Lychrel inferiores a',M
print 'Tiempo de cómputo:',time.time()-t0,'segundos'
Hay 249 números Lychrel inferiores a 10000
Tiempo de cómputo: 0.345999956131 segundos
Como he dicho, la solución no son en realidad números de Lychrel, y asimismo el encontrarlos o negarlos quizás dependa más de una astuta demostración que de potentísimos ordenadores.

La millonésima permutación lexicográfica

De nuevo Project Euler, tal vez con el programa más difícil que me he encontrado hasta ahora. Desde que empecé a programar, la recursión (recursividad o recurrencia) me ha parecido de las formas más difíciles y sin embargo más elegantes de programar (programación funcional). Para hacer permutaciones me parece inevitable utilizar recursión. Cito el problema 24 de Project Euler:
Una permutación es una disposición ordenada de objetos. Por ejemplo, 3124 es una  posible permutación de los dígitos 1, 2, 3, 4. Si todas las permutaciones se listan numérica o alfabéticamente, a esta disposición le llamamos orden lexicográfico. Las permutaciones lexicográficas de 0, 1, 2 son: 012   021   102   120   201   210. ¿Cuál es la millonésima permutación lexicográfica de los dígitos 0, 1, 2, 3, 4, 5, 6, 7, 8, 9?
El problema parece sencillo; después de programar los problemas 19 y 28 de forma casi trivial, me topé con este problema que simplemente no podía lograr. Prácticamente siempre empiezo un programa de la forma más sencilla, lo pruebo así y lo llevo después a la forma complicada. El ejemplo que da Project Euler puede resultar engañoso, y como yo, cualquiera puede caer en un pozo sin fondo.

Sabemos que el número de permutaciones posibles de un conjunto de $n$ objetos es igual a $n!=n\cdot(n-1)\cdot3\cdot2\cdot1$. El problema de comenzar con casos sencillos, es que, por ejemplo para $n=3$, un simple bucle anidado de 2 niveles es suficiente para repetir 3 veces un bucle que se repite 2 veces, obteniendo las $3!=6$ permutaciones buscadas, en cambio, si tomamos el caso $n=4$, necesitaríamos un bucle de 3 niveles que hiciera 4 veces un bucle que repite 3 veces otro bucle que se repite 2 veces, obteniendo $4!=24$ permutaciones... y así sucesivamente para cualquier $n$ usando este razonamiento, lo cual no es práctico ni mucho menos, recordando que el problema que queremos resolver nos pide $10!$ permutaciones.

Bueno, pues ya que sabemos que se obtendrán $10!$ permutaciones, debemos recordar cómo utilizar la técnica de recursividad. Regularmente cuando se introduce este concepto se utiliza el ejemplo del factorial. En Python se ve algo así:

def f(n):
if n > 1 :
return (n*f(n-1))
else:
return 1
Bien pues se necesita usar este tipo de razonamiento pero buscando formar permutaciones. Una buena forma de lograrlo es intentarlo primero con un conjunto de cuatro objetos, lo más sencillo a empezar con 0123. Un razonamiento que funciona es intercambiar un par externo e ir avanzando hacia el siguiente término, como se observa el patrón que a continuación muestro. Las primeras permutaciones se verían así:
0123
0132
0213
0231
0312
0321
1023
1032
1203
1230
1302
1320
2013
$\vdots$
Notamos fácilmente que para la primer posición de izquierda a derecha, el 0 se mantiene fijo y los términos adyacentes permutan $6=3!$ veces, y de igual modo, si nos fijamos en la posición siguiente, en todas las permutaciones que comienzan con 0, cada objeto tiene $2=2!$ permutaciones, y esto fácilmente lo podemos predecir para cualquier combinación siguiente. Si observáramos una posición aún anterior, veríamos que 0, 1, 2, y 3 forman exactamente $4!$ permutaciones, que es lo que buscamos. También es fácil notar que no necesitamos ordenar permutaciones para formar un orden lexicográfico si comenzamos con una cadena ordenada, i.e. 0123456789.

Ahora el trabajo arduo (al menos para mí lo fue) es describir esto en un lenguaje estricto y obtener lo que se busca. Recomiendo primero encontrar el patrón que describen las permutaciones y luego intentar codificarlo. Casi siempre que resuelvo un nuevo problema de Project Euler parece que ha sido el más retador, sin embargo éste se ha llevado definitivamente el título, además de que la recursividad como dije, no me parece algo tan sencillo.

Python (y supongo que otros lenguajes más) tienen funciones predefinidas que dan permutaciones, sin embargo ¿qué de divertido habría en sólo buscar la millonésima de éstas, o si acaso tener que ordenar las permutaciones?.

Este tipo de problemas es digno de discusión con profesores o colegas/compañeros si es que te estancas; en mi caso pasé 2 días y 2 noches pensando y lidiando con la recursividad, aunque el patrón sea muy sencillo de descifrar.

Mi código es un tanto feo aunque no es tan largo y tuve algunos problemas con las posiciones del arreglo que contiene cada permutación. Dentro del programa cree un archivo en el cual escribí las permutaciones y al final la millonésima permutación lexicográfica, lo que como vimos arriba sería innecesario considerando que $10!>1,000,000$ pudiéndonos detener entonces antes de que el programa se ejecute totalmente.

Aquí dejo las permutaciones del conjunto 012345678 para que te puedas guiar con más datos acerca del patrón formado y el algoritmo que debes codificar.

HAPPY CODING!

:)

Números Amigos

Imagen (si ves este texto, recarga la pág)

A todos nos resulta familiar el conjunto de los números naturales o el conjunto de los números enteros, aunque a algunos niños cueste tanto trabajo dominar cuando se introducen los enteros negativos, en fin, los conjuntos de números más comunes, que conocemos como $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{Q}$ o $\mathbb{C}$. Pero hay más formas de clasificar a los números, como por ejemplo en números hambrientos, números vampiros, números narcisistas, etc… que cumplen ciertas características que los hacen peculiares (Gaussianos tiene una entrada muy buena sobre estos tipos de números).

Una de estas tantas clasificaciones peculiares son los números amigos. Una pareja de números amigos son dos números enteros positivos a, b tales que a es la suma de divisores propios de b y b es la suma de divisores propios de a, con a diferente de b (cuando a=b se les llama números perfectos o amigos de sí mismos).

La verdad es que de nuevo, me encontré con ellos en un problema de Project Euler, en el cual se dice:
Sea d(n) la suma de divisores propios de n. Si d(a)=b y d(b)=a, donde a≠b, entonces (a,b) es una pareja amigable y tanto a como b se llaman números amigos.
Los números amigos al parecer se conocen desde la época de Pitágoras y desde entonces han sido objeto de estudio de los matemáticos, de los cuales Euler participó en la generalización de una regla para obtenerlos aunque al parecer no funciona para cualquier número y da al menos dos parejas que no son números amigos.

Aún más sutil es el caso de las parejas regulares de números amigos. Sean (a,b) una pareja amigable con a<b y sea w el mayor común divisor de esta pareja tal que a=Aw y b=Bw. Si A y B son primos relativos y libres de cuadrados, producto de i y j factores primos, respectivamente, entonces la pareja (a,b) es regular y se dice que es de tipo (i,j). Si esto no se cumple simplemente se dice que son irregulares.

El ejemplo más sencillo a exponer sería el de la primer pareja de números amigos: (220,284). Usando el enunciado del problema de Project Euler: podemos notar que los divisores propios de 220 son 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 y 110; así entonces, d(220)=284. Los divisores propios de 284 son 1, 2, 4, 71 y 142; entonces d(284)=220. Además de esto, notamos que el mayor común divisor de ambos es w=4, entonces nota que A=55, B=71, los cuales son primos relativos y libres de cuadrados, producto de 2 y 1 factores primos, A=11$\times$5 y B=71, entonces (a,b)=(${4\times11\times5}$, ${4\times71}$), esto es, la pareja (220,284), es una pareja amigable regular de tipo (2,1). Esta pareja es la más antigua conocida y no hay pareja alguna anterior a ella.

Una forma usada para encontrar números amigos, números perfectos y números sociales es la sucesión alícuota en la que cada término es la suma de los divisores propios del término anterior. Apenas con el advenimiento de las computadoras se ha potencializado la búsqueda de números amigos y ya se conocen millones de ellos.

Bueno pues comparto el enunciado completo del problema de Project Euler y un código en Java que soluciona el problema. No hice nada extraordinario y el programa funciona bien para números amigables debajo de 100,000 tardando a lo mucho 1 minuto.

****¿Puedes realizar un programa que diga cuándo un par amigable es regular y de qué tipo?

Problema 21 (Project Euler):
Sea d(n) la suma de divisores propios de n (números estrictamente menores que n que dividen exactamente a n). Si d(a)=b y d(b)=a, donde a≠b, entonces a y b es una pareja amigable y tanto a como b se llaman números amigos. Por ejemplo, los divisores propios de 220 son 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 y 110; así entonces, d(220)=284. Los divisores propios de 284 son 1, 2, 4, 71 y 142; entonces d(284)=220. Evalúa la suma de todos los números amigos menores a 10000.

public class euler21
{
public static void main(String[] asrg)
{
final int MAX=100000;
int a[],i,j,x=0;
long s1=1,s2=1,t=0;
a=new int[MAX];
for(i=6;i<MAX;i++,s1=1,s2=1)
{
for(j=2;j<=(i/2);j++)
{ if(i%j==0) {s1+=j;} }
for(j=2;j<=(s1/2);j++)
{ if(s1%j==0) {s2+=j;} }
if(s2==i && s2!=s1) {a[x]+=i;x++;}
}
System.out.println("Los numeros amigos debajo de "+MAX+" son:");
for(i=0;a[i]!=0;i++)
{
t+=a[i];
System.out.println(a[i]);
}
System.out.println("\nSe encontraron en total "+i+"numeros amigos");
System.out.println("Y la suma de numeros amigos debajo de "+MAX+" es: "+t);
}
}

The first triangle number to have over a thousand divisors

So this is the twelfth problem of Project Euler (I mean with a 1000 instead of a 500), and though the reasoning to get it right seems rather easy, it becomes puzzling when you make a little-hard-way program, run it, and watch it go blank for several looooooooong hours. So here’s the original description of the problem (register in Project Euler & keep you own record!):

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be

1 + 2 + 3 + 4 + 5 + 6 + 7 = 28

The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:
1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

So… the hardest way to do it is by checking divisors up to the nth triangle number. We obviously know that doesn’t work and we’re so clever that we notice right away that a fine optimization would be to check up to n/2. Ok! So go and try doing it this way! You’ll only get valuable time wasted and you'll wait for hours for the naughty number to magically appear (I mean, it’s not that this was my case ¬¬). So, how to do it? I suggest looking for a pattern in the sequence of divisors. A hint is that all triangle numbers seem to have an even number of divisors (I haven’t proved this sentence and I don’t know if someone already has).

So, work it out and when you’re ready, check out the algorithm below to see how I did it (it’s a really short one! and it works fine with MAX=1000).
long long int t=28, c=0, i, j;
for(i=8; c<MAX; i++)
{
c=0;
t+=i;
for(j=2; j<=t/j; j++)
{
if(t%j==0) {c++;}
}
c=2*(c+1);
}