You are currently browsing the monthly archive for octubre 2012.

Toda aplicación es habitual que tenga un fichero de configuración. Lo primero que hace el programa cuando empieza a funcionar es leer este fichero en el que se especifican una serie de parámetros que determinan una cierta configuración de funcionamiento.

Cuando somos nosotros los que estamos desarrollando una aplicación, la pregunta que se nos pasa por la cabeza es: ¿Qué formato debe tener un fichero de configuración? Ciertamente, dado que nosotros somos los desarrolladores, tenemos total libertad para determinar cual sera éste. Sin embargo, el hecho de seleccionar uno u otro nos permitirá aproximar nuestra aplicación a la manera estándar de funcionar el software, teniendo también la posibilidad de utilizar herramientas desarrolladas para estos estándares de facto.

Volviendo a nuestra libertad de elección, la primera manera que se nos ocurre es tener un fichero de texto con un formato determinado. Lo único que tenemos que hacer es desarrollar un código capaz de leer y escribir ficheros con este formato.

De unos años a esta parte, los formatos mas utilizados para los ficheros de configuración son los ficheros INI, XML, YAML y JSON.

INI: es el mas sencillo. Solo admite dos niveles y no esta preparado para binarios.

;config valSPH
...
numParticles=10000
kernel=cubic
odeSolver=rungeKutta23
...

XML: muy extendido. Soporta múltiples niveles, binarios… Hay muchas herramientas para parsearlo, por ejemplo Xerces de Apache, o también es fácil programarse uno basado en SAX.


<!-- config valSPH -->
...
<numParticles>10000</numParticles>
<kernel>cubic</cubic>
<odeSolver>rungeKutta23<odeSolver>
...

JSON: Dispone de MIME type, “application/json” con codificación y decodificación nativa en navegadores. Muy utilizado en projectos Ajax.

{

...
"numParticles": 10000,
"kernel": "cubic",
"odeSolver": "rungeKutta23",
...

}

JAML: nace como una generalización de JSON. Existe soporte en C++ con yaml-cpp

---# config valSPH
...
numParticles:  10000
kernel:        cubic
odeSolver:     rungeKutta23
...

Anuncios

Actualizamos la minianimación que teniamos incluyendo los algoritmos y las estructuras de datos necesarias para la NNPS.

Con el fin de poder comprobar la determinación de vecinos, el color de cada péndulo depende, en un instante dado, del número de vecinos que tiene a una distancia menor o igual a h.

En la siguiente animación se utiliza una LinkedList:

En el método SPH necesitamos conocer cual es su conjunto de vecinos de una partícula a. Este viene determinado por la smoothing length h, uno de los parámetros de las funciones kernel.

Supongamos que tenemos N partículas. Si N es lo suficientemente pequeño, lo único que tenemos que hacer es, para cada partícula, recorrer el resto de partículas calculando su distancia a nuestra partícula y tener en consideración solo aquellas que esten a distancia menor que h (all-pair search). La complejidad del algorítmo es O(N^2).

Sin embargo, cuando el número de partículas crece, y que lo ha de hacer pues en eso se basa la potencia del método, esta manera de trabajar es mejorable.

A continuación veremos algunas opciones, pero la idea general es que necesitamos tener una descomposición espacial de nuestro escenario que nos permita encontrar de manera eficiente el vecindario de cada una de nuestras partículas y que tendremos que mantener actualizada a medida que evoluciona el sistema.

Otra opción es la linked-list. Se trata de superponer una malla con un tamaño de celda 2h, de manera que, dada una partícula, solo tendremos que buscar las partículas que se encuentran en las celdas adyacentes. En C++ tenemos las clases list y vector. La primera por nombre y la segunda porque la clase vector puede crecer dinámicamente. Solo podremos utilizarlas si la h es constante para todas las partículas. La complejidad, con un número suficientemente pequeño de partículas por celda, es O(N).

Finalmente, podemos trabajar con árboles. Los árboles permiten complejidades genéricas en búsquedas de O(N \log N) y permiten trabajar con h diferentes para cada partícula. En astrofísica, por ejemplo, podríamos tener un arbol binario de búsqueda, pues tenemos que simular objetos autogravitantes y necesitamos resolver un problema de n-cuerpos, que parece ser que es lo que mas tiempo de CPU consume. Para hacer esto de manera eficiente necesitamos un BST. Se trata de utilizar este mismo para realizar las búsquedas. La clase map de C++ se implementa con un arbol binario balanceado (AVL).

Existen otras estructuras de datos que permiten hacer descomposiones espaciales que podria ser interesante analizar (UB-tree, Octrees, BST, kd-tree, R-tree).

En la siguiente simulación tenemos 1000 partículas. Cada partícula a representa un péndulo de masa m_a y posición (x_a(t),y_a(t),z_a(t)), con z_a(t) constante, sometido a las fuerzas:

  1. F_1 = -m_ag, la fuerza de la gravedad, con g constante.
  2. F_2, la fuerza elática de la suspensión, de sentido opuesta a la posición del péndulo y de magnitud k_1d, con k_1 > 0 y d la variación de la longitud del péndulo respecto de su longitud en reposo l.
  3. F_3, la fuerza de fricción, con sentido opuesta a la velocidad del péndulo y magnitud k_2v, con v = \sqrt{(x_a')^2+(y_a')^2} y k_2 > 0.

La ley de Newton F = ma queda, para cada uno de los péndulos, es decir, para cada partícula, como el sistema de dos ecuaciones diferenciales de segundo orden:

m_a \begin{bmatrix} x_a''(t) \\ y_a''(t) \end{bmatrix} = F_1 + F_2 + F_3

que es:

m_ax_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})x_a - k_2\frac{x'_a}{\sqrt{x_a'^2+y_a'^2}}

m_ay_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})y_a - k_2\frac{y'_a}{\sqrt{x_a'^2+y_a'^2}}

Podemos expresar este sistema como un sistemad de ecuaciones de primer orden introduciendo incongnitas adicionales: u_a = x_a' y v_a = y_a', con lo que nos queda:

\left \{  \begin{array}{rcl}  x_a' &=& u_a \\  y_a' &=& v_a \\  m_au_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)x_a - k_2\frac{u_a}{\sqrt{u_a^2+v_a^2}} \\  m_av_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)y_a - k_2\frac{v_a}{\sqrt{u_a^2+v_a^2}} - m_ag \\  \end{array}  \right .

En la siguiente simulación tomamos m_a aleatoria para cada partícula a, por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos l=1, k_1 = 10000 y k_2=1 en todos los péndulos. Cada uno de éstos se distribuye a lo largo de una espilar cilíndrica y las velocidades iniciales son, para todas ellas, u_{a_0}=0, v_{a_0} = -0.7 y v_{z_0} = 0. El intervalo de tiempo es desde t_i = 0 hasta t_f=10 que se divide en 400 subintervalos.

El resultado es el siguiente:

Para poder generar los call graph y los caller graph desde doxygen necesitamos tener instalado Graphviz.

Aquí tenemos la direccion desde donde nos lo podemos bajar. En el caso de mac, es un pkg que instala el fichero binario dot en el directorio /src/local/bin. Tendremos que añadir este path en la variable DOT_PATH en la pestaña Expert de doxygen.

A continuación, unos ejemplos del tipo de diagramas generados:

Además, en la segunda captura podemos observar que las clases pueden mostrarse con estilo UML.

A medida que codificamos y probamos nuestras clases las iremos documentando mediante doxygen.

La instalación y utilización básica de doxygen es extremadamente sencilla. El programa está disponible aquí.

A continuación tenemos un ejemplo del formato de los comentarios añadidos y el tipo de salida generado.

A continuación nuestra primera animación creada con VisIt a partir de ficheros silo generados a partir de llamadas a métodos de nuestras primeras clases.

Temporalmente, aunque es mejor utilizar nuestras partículas simplemente como nodos de interpolación a la hora de visualizar, las mostraremos como puntos.

A partir de este momento, podremos tener una referencia visual de nuestras partículas que nos ayudará mucho a la hora de evaluar como de bien estamos haciendo las cosas.

Sin mas dilación, nuestra minianimación:

No es mas que un conjunto de mil partículas desplazandose sobre un cilindro. La animación consta de 100 ficheros .silo y hemos generado el mpeg con el wizard de VisIt.

Ya dedicamos un post a VisIt. En el creamos una animación a partir de un conjunto de ficheros existentes.

¿Como generamos ficheros .silo desde nuestro código? En primer lugar, necesitamos tener las librerias correspondientes, que las podemos conseguir ejecutando:

./build_visit --console --no-visit --no-thirdparty  --thirdparty-path /usr/local --silo --hdf5 --szip

desde el terminal, y donde build_visit es un fichero que podemos conseguir aqui.

A continuación, necesitamos incluir la libreria <silo.h> en nuestro fichero y añadir:

PLATFORM=i386-apple-darwin10_gcc-4.2

SZIP_DIR=/usr/local/szip/2.1/$(PLATFORM)
SZIP_CPPFLAGS=-I$(SZIP_DIR)/include
SZIP_LDFLAGS=-L$(SZIP_DIR)/lib
SZIP_LIBS=-lsz

HDF5_DIR=/usr/local/hdf5/1.8.4/$(PLATFORM)
HDF5_CPPFLAGS=-I$(HDF5_DIR)/include $(SZIP_CPPFLAGS)
HDF5_LDFLAGS=-L$(HDF5_DIR)/lib $(SZIP_LDFLAGS)
HDF5_LIBS=-lhdf5 $(SZIP_LIBS) -lz

SILO_DIR=/usr/local/silo/4.6.2/$(PLATFORM)
SILO_CPPFLAGS=-I$(SILO_DIR)/include $(HDF5_CPPFLAGS)
SILO_LDFLAGS=-L$(SILO_DIR)/lib $(HDF5_LDFLAGS)
SILO_LIBS=-lsiloh5 $(HDF5_LIBS) -lm

LDFLAGS=$(LDFLAGS) $(SILO_LDFLAGS)
LIBS=$(SILO_LIBS)
CPPFLAGS=$(CPPFLAGS) $(SILO_CPPFLAGS)

a nuestro Makefile

Cuando eramos niños, nos pasabamos el curso esperando la llegada de las vacaciones. Sin embargo, a medida que pasaban, cada vez tenías mas ganas de volver a las clases.

En esto de programar pasa un poco lo mismo. Cuando estas metido en un proyecto en el que te pasas dias y dias escribiendo, probando código y peleandote con los ordenadores, lo que mas desea uno en el mundo es acabar con aquello. Sin embargo, una vez terminado y pasado un tiempo, después de disfrutar de un merecido descanso y del placer y de la potencia de ver aquello funcionando, nos entra el gusanillo de volver a la carga.

Existen muchos modelos para el ciclo de vida del software, es decir, diferentes maneras de enfrentarte al reto de como desarrollar sofware. Sin embargo, con el tiempo y despues de acumular muchas “horas de vuelo”, uno ya tiene mas o menos definida su pequeña estrategia. Son cosas de sentido común, pero que no va mal tenerlas presente.

Para empezar, es muy importante tener una visión global del sistema, de manera que, nuestro primer diseño, es interesante que lo abarque todo. Para ello van bien, por ejemplo, realizar el diagrama de casos de uso, el diagrama de clases, y algunos diagramas de comunicación y de secuencia. Además, existe software que sera capaz de generarnos, a partir de ellos, la arquitectura de la aplicación, es decir, el conjunto de todos los ficheros de clases que necesitaremos y la interrelación entre ellos.

Sin embargo, en el momento de empezar a implementar el software, a escribir el programa, es mejor seguir una estrategia en espiral.

cuya idea de fondo es sencilla: no pretender programarlo todo de golpe sino hacerlo de manera incremental. Para ello, empezamos seleccionando un pequeño conjunto del global de los objetivos de nuestro software, seleccionaremos las clases que nos permitiran asumirlos, las implementaremos y las probaremos. Llegados a este punto habremos dado una vuelta en la espiral ya que estaremos en condiciones de volver a empezar la misma rutina con otro pequeño conjunto de objetivos pero con la diferencia de que nuestro programa habrá crecido con respecto a la vuelta anterior y lo que allí nos propusimos estará funcionando.

Pues nada, teniendo ésto presente y sin mas preámbulos, nos ponemos el mono de programador y, por fin, ¡a programar!

En el artículo A modified SPH approach for fluids with large density differences de F. Ott y E. Schnetter se explica una nueva modificación del método SPH que permite la interacción entre fluidos con densidades muy diferentes.

En la aproximación del SPH estandar tenemos:

\rho_i = \sum_j m_j W_{ij}

con W_{ij}:=W(\boldsymbol{x}_i - \boldsymbol{x}_j), o también:

\frac{d}{dt}\rho_i = \sum_j m_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

donde \boldsymbol{v}_{ij}:=(\boldsymbol{v}_i - \boldsymbol{v}_j) y \nabla W_{ij}:=(\nabla W)(\boldsymbol{x}_i - \boldsymbol{x}_j), que tiene la ventaja de que los valores iniciales para \rho_i pueden escogerse libremente permitiendo la modelización de discontinuidades.

Ellos proponen utilizar,

\frac{d}{dt}\rho_i = m_i \sum_j \boldsymbol{v}_{ij} \cdot \nabla W_{ij}

deducida suavizando la densidad de partículas n_i = \frac{1}{V_i} en lugar de la densidad de masa \rho_i, donde V_i = \frac{m_i}{\rho_i} y \frac{d}{dt}m_i = 0.

justificandolo con tests positivos en tres problemas: ecuación de advección, onda sonora y tubo de choque.

En el review que Rosswog hace sobre el método SPH, en especial en sus aplicaciones a la astrofísica, hay un apartado dedicado al tratamiento de los shocks. El tratamiento de los strong shocks es uno de los puntos débiles del método SPH (penetración de partículas parcialmente resuelto mediante XSPH). En este apartado comenta que, básicamente, exiten dos maneras de tratarlos:

  1. Utilizar soluciones (analíticas) de un problema de Riemann entre dos partículas adyacentes.
  2. Ampliar la discontinuidad hasta una longitud resoluble numéricamente, para lo que se añaden terminos disipativos artificiales extra al flujo.

Comenta que la mayoria de las implementaciones de SPH utilizan esta segunda aproximación, consiguiendo esta amplicación o mediante el método de discretización utilizado, utilizando métodos que modifican las ecuaciones originales de una manera que nos puede venir bien, tipo el esquema Lax; o mediante la adición de terminos ad-hoc, algo bastante habitual, de tipo presión.

Pero también existen métodos que implementan la primera, y cita dos artículos:

  1. El artículo Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver de Inutsuka y Shu-Ichiro.
  2. El artículo Implementations and tests of Godunov-type particle hydrodynamics de Cha y Whitworth.

En el primero comenta como consiguen, tras algunos intentos fallidos,  introducir Riemann Solver en las ecuaciones de manera similar a como se consigue en los métodos basados en malla, manteniendo la formulación conservativa de las mismas. Para empezar, comentan que mientras el error en la mayoria de los métodos basados en malla proviene de la truncación de la serie de de Taylor, en SPH, si consideramos:

\langle f \rangle(\boldsymbol{x}) := \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

como la convolución de la función f(x) con el kernel. Mediante Taylor, de f alrededor de x=x_i, podemos ver que:

\langle f\rangle(\boldsymbol{x}_i) - f(\boldsymbol{x}_i) = \frac{h^2_{ef}}{4} \nabla^2 f(\boldsymbol{x}) + O(h^4_{ef}).

donde

h^2_{ef}:=2 \int x^2 W(\boldsymbol{x},h)d\boldsymbol{x},

por lo que el error proviene de la convolución con la función simétrica W.

También se tiene:

\langle \frac{\partial f}{\partial x} \rangle (\boldsymbol{x}) = \int \frac{\partial f(\boldsymbol{x'})}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= - \int f(\boldsymbol{x'}) \frac{\partial}{\partial x'} W(\boldsymbol{x'}-\boldsymbol{x},h)d\boldsymbol{x'}

= \frac{\partial}{\partial x} \int f(x')W(x'-x,h)dx'

De manera que \nabla f = \nabla \langle f \rangle. Si definimos:

\rho(\boldsymbol{x}):=\sum_j m_j W(\boldsymbol{x}-\boldsymbol{x}_j,h),

entonces:

\sum_j \frac{m_j}{\rho(\boldsymbol{x})}W(\boldsymbol{x}-\boldsymbol{x}_j,h) = 1,

\sum_j m_j \nabla \big [ \frac{W(\boldsymbol{x}-\boldsymbol{x}_j,h)}{\rho(\boldsymbol{x})} \big ]= 0

que son la clave del método, pues ahora se puede hacer:

f_i:= \langle f \rangle (\boldsymbol{x}_i) = \int f(\boldsymbol{x'})W(\boldsymbol{x'}-\boldsymbol{x}_i,h)d\boldsymbol{x'}

= \int \sum_j m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'

\sum_j \Delta f_{i,j}

donde:

\Delta f_{i,j}:=\int m_j \frac{f(\boldsymbol{x}')}{\rho(\boldsymbol{x}')}W(\boldsymbol{x}'-\boldsymbol{x}_i,h)W(\boldsymbol{x}'-\boldsymbol{x}_j,h)d\boldsymbol{x}'.

En el método SPH estandar nos encontramos, por una parte, con:

\sum_j \frac{m_j}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}-\boldsymbol{x}_j,h) \approx 1

y por otra, con:

\Delta f_{i,j} \approx m_j \frac{f(\boldsymbol{x}_j)}{\rho(\boldsymbol{x}_j)}W(\boldsymbol{x}_i-\boldsymbol{x}_j,h)

y que, evidentemente, son aproximaciones mas pobres a las anteriormente obtenidas.

A continuación derivan las ecuaciones para la posición:

\ddot{\boldsymbol{x}}:=\int \frac{d\boldsymbol{v}(\boldsymbol{x})}{dt}W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= - \int \frac{1}{\rho(\boldsymbol{x})} \nabla P(x) W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

y para la energía

\dot{u}_i := \int \frac{du(\boldsymbol{x})}{dt} W(\boldsymbol{x}-\boldsymbol{x}_i,h)d\boldsymbol{x}

= \sum_j m_j \int \frac{P}{\rho^2}[\boldsymbol{v}-\dot{\boldsymbol{x}}_i] \cdot [\nabla W(\boldsymbol{x}-\boldsymbol{x}_i,h)W(\boldsymbol{x}-\boldsymbol{x}_j,h) - W(\boldsymbol{x}-\boldsymbol{x}_i,h) \nabla W(\boldsymbol{x}-\boldsymbol{x}_j,h) ] d\boldsymbol{x},

que se modifica a

$latex$

para tener en cuenta la disipación

En el segundo artículo define e implementa el método Godunov-type Particle Hydrodynamics (GPH). Básicamente GPH = SPH + Godunov upwind schema.

Hablemos de generalizaciones, cosa esencial en matemáticas:

  • Desde el punto de vista del análisis funcional, ¿qué pasa si nuestro espacio tiene asociada una métrica no definida positiva?
  • En el caso finito, ¿qué pasa cuando en lugar de trabajar con \mathbb{R}^n o \mathbb{C}^n tengo variedades diferenciables mas generales?
  • En analisis funcional tengo espacios de funciones donde estas cumplen una propiedad asociada con una medida y especificada normalmente mediante una integral. ¿Puedo tener espacios de variedades cumpliendo propiedades de este tipo? (pues tiene sentido hablar de integración en variedades orientadas)
  • Al hablar de variedades de Riemann, asociamos una métrica a una variedad para poder hablar de distancias, areas, angulos, etc. En realidad, asociamos un producto escalar, un tensor 2 veces covariante, del que deriva una norma a partir de la que podemos especificar una distancia.

Teorema de Laurent.

Definición (serie de Laurent): Sea z_0 \in \mathbb{C}. Una serie de Laurent en z_0 es formalmente una expresión de la forma:

\sum_{n \in \mathbb{Z}} a_n (z-z_0)^n con a_n \in \mathbb{C}.

Llamaremos parte analítica a \sum_{n=0}^\infty a_n (z-z_0)^n y parte principal a \sum_{n=1}^\infty a_n (z-z_0)^{-n}. Diremos que la serie es convegente en z si lo son su parte analítica y su parte principal.

Definición (anillo): Llamaremos anillo centrado en z_0 de radio menor r y radio mayor R a:

A(z_0;r,R):= \{ z \in \mathbb{C}: r < |z-z_0| < R\}.

En el caso que r=0 entonces A(z_0;0,R) = D(z_0,R)-\{z_0\} = D'(z_o,R) tenemos un disco perforado. Si r=0 y R=\infty entonces A(z_o;0,\infty) = \mathbb{C}-\{z_0\}. Finalmente, si r \neq 0 y R = \infty entonces A(z_0,r,\infty = \mathbb{C}-\overline{D(z_o,r)}.

Teorema (de Laurent): Si f \in \mathcal{H}(A(z_0;r,R)) entonces existe a_n \in \mathbb{C} tal que:

f(z)= \sum_{n \in \mathbb{Z}} a_n(z-z_0)^n para todo $z \in A(z_0;r,R)$.

demostración:

Clasificación de singularidades aisladas

Teorema de los residuos.

Definición (función meromorfa): Sea f:U-A \longrightarrow \mathbb{C} una función analítica con U \subset \mathbb{C} un abierto y A el conjunto de singularidades aisladas de f. Diremos que f es una función meromorfa en U.

Teorema (de los residuos)

Integrales tipo I, II, III, IV, V i VI

Postulado 1: La representación del estado de un sistema físico.

P1. La máxima información posible sobre un sistema físico en un instante dado t es su estado cuántico, que se representa por un vector (ket) de norma 1 y de fase arbitraria en un espacio de Hilbert separable.

Superposición en un producto tensorial de espacios de Hilbert

Postulado 2: La representación de las magnitudes medibles.

P2. Cualquier magnitud del sistema que en pricipio se pueda medir tiene asociada un operador lineal autoadjunto definido sobre el espacio vectorial de los estados. La totalidad de los autovalores recibe el nombre de espectro y los autovectores definen una base del espacio de Hilbert.

Conjunto completo de observables compatibles

Postulado 3: El resultado de la medida.

P3. El resultado de una medida de un observable A sobre un estado |\psi \rangle solo puede ser uno de sus autovalores. La probabilidad de que al medir A, de descomposicón espectral A=\sum_i \lambda_i \Pi_i,\, (\lambda_i \neq \lambda_j), sobre un estado |\psi \rangle obtengamos el resultado \lambda_i es:

\mathcal{P}_{| \psi \rangle} (A:\lambda_i) = ||\Pi_i | \psi \rangle ||^2 = \langle \psi | \Pi_i | \psi \rangle

donde A:\lambda_i indica que la medida del observable A nos da el valor \lambda_i

Dualidad onda-particula

Incertidumbre de Heisenberg

Postulado 4: La proyección, reducción o colapso.

P4. Cualquier estado | \phi \rangle sobre el que se realiza una medida del observable A que es filtrante respecto del autovalor \lambda_i, pasa a estar inmediatamente despues de la medida en el estado

\frac{\Pi_i | \phi \rangle}{||\Pi_i | \phi \rangle||}

con probabilidad ||\Pi_i | \phi \rangle||^2, o se destruye durante el proceso de medida.

La no localidad de la mecánica cuántica

Postulado 5: La evolución temporal.

P5. Entre medidas, el sistema evoluciona según

i\hbar \frac{\partial}{\partial t} | \varphi(t) \rangle = H | \varphi (t) \rangle

donde H es el operador hamiltoniano del sistema.

Postulado 6: los operadores de posición y de momento.

P6. Los operadores de posición, en coordenadas cartesianas, y de momento satisfacen las reglas de conmutación:

[X_i,X_j]=0, \, [P_i,P_j]=0, \, [X_i,P_j]=i \hbar \delta_{ij} I.

La función de onda en el espacio de posiciones

La función de onda en el espacio de momentos

Dimensiones y teorema de Ehrenfest

Los valores de los operadores X_i y P_i evolucionan según las leyes clásicas pero con fuerzas medias.

Ayer asistí al primer Friday’s miniWorkshop del IVICFA (L’institut Valencià d’Investigació Coorporativa de Física Avançada) que trataba sobre supercomputación y computación GRID. Forma parte de un ciclo de seminarios sobre “Fronteras de la Física” y contaba ayer, entre otros atractivos, con la apertura de los mismo por parte de Ian Bird, Project Leader of the Worldwide LHC Computing GRID (WLCG), con “Petabyte scale computing for the LHC”, o el seminario “Modeling complex solar magnetodynamical phenomena using supercomputing and visualization techniques” de Fernando Moreno Insertis, del Instituto de Astrofísica de Canarias (IAC) y PI of the European Solaire Network.

Aunque ya hicimos algunos comentarios en un post anterior, aprovechamos la ocasión para volver a hablar sobre estos temas y así aclararnos las ideas.

Un computador computa, por lo que un supercomputador supercomputa ;-). Como el elemento fundamental para computar en la arquitectura Von Neumann son las CPUs, en un supercomputador dispondremos de multiples de éstas. El otro elemento fundamental en la arquitectura Von Neumann, que es la que las convierte en máquinas de proposito general, es la memoria. En función de como las CPUs comparten esta memoria, nos encontramos con dos modelos de supercomputación: memoria compartida y memoria distribuida. Las dos interfaces que nos permiten programar en estos ambientes son OpenMP y MPI respectivamente.

Tenemos que aclarar que, cuando hablamos de supercomputación, nos estamos refiriendo a que estas entidades, múltiples procesadores y memorias, existen físicamente, pues podemos encontrarnos todos estos elementos de manera virtual en ordenadores mas sencillos con sistemas operativos multiproceso.

La última tendencia en arquitecturas paralelas son las arquitecturas vectoriales, arquitecturas SIMD: una instrucción múltiples datos, debido a su enorme desarrollo en la evolución de las tarjetas gráficas, que es donde aparecen de manera natural al tener que aplicar la misma operación a múltiples píxeles. CUDA es un estandar de facto en la extensión de estas operaciones a cualquier tipo de datos no necesariamente gráficos.

Por tanto, en el mundo de los supercomputadores actuales es fácil encontrarnos con arquitecutras paralelas heterogeneas en donde conviven simultanemente la memoria compartida con la memoria distribuida de múltiples unidades de proceso que soportan un juego de instrucciones vectorial. La paralelización de estos códigos se realiza ad-hoc, pero, al igual que sucedió con el triunfo de las arquitecturas RISC sobre las CISC, la última palabra es posible que la tengan los compiladores, aunque es un hecho objetivo que la paralelización automática es instrínsecamente muy compleja.

Un compilador, aunque se suele asociar a los traductores de lenguajes de alto nivel a lenguaje máquina, en realidad son traductores entre dos lenguajes, sean del tipo que sean, por lo que es totalmente lógico pensar en compiladores de programas secuenciales en un lenguaje determinado a programas paralelos heterogeneos en el mismo lenguaje. Mientras se programó en lenguaje máquina, las arquitecturas CISC dominaron el mercado. Con la aparición de los lenguajes de alto nivel y de compiladores potentes para los mismos que solo utilizaban un subconjunto muy reducido del amplísimo conjunto de instrucciones disponibles, entraron en escena las arquitecturas RISC y desplazaron a las primeros.

Por otra parte, ¿cúal es la propuesta del GRID computing? Pués su aparición es, al igual que por ejemplo los protocolos TCP/IP de Internet, una solución de facto a un problema existente y es, por una parte, el hecho de que existan en un momento determinado multitud de recursos repartidos a lo largo del planeta, y por otra, a la posibilidad de compartirlos de manera transparente por multitud de usuarios igualmente dispersos.

Por ejemplo, la realidad es que un grupo de la Universidad de las Palmas de Gran Canaria tiene un supercomputador fruto de sus necesidades en un determinado momento y otro grupo de la University of Tasmania tiene otro con los mismos recursos. Si en un momento dado se dan cuenta de que podrían compartir sus recursos, lo cual a priori siempre es positivo desde el punto de vista de la teoría de juegos, la tecnología grid ofrece una capa de abstracción, la middleware grid, gracias a la cual pasamos a tener un único sistema con el total de los recursos.

Pensando trivialmente, solo por aclarar ideas, y sabiendo que es un caso totalmente irreal, si ambas máquinas solo se utilizasen durante las ocho horas laborables locales, supongamos de 8h a 18h, obviamente los dos grupos salen muy beneficiados, pués con una diferencia horaria de 10 horas no habría conflictos de acceso y los dos pasarían a disponer del doble de recursos de los que tenían.

Por tanto, los supercomputadores existen y existirán, pues son la mejor solución a nivel local, pero la tecnología GRID, sin entrar en consideraciones sobre Cloud, es la mejor solución a nivel global, ya que permite la interconexión transparente de estos óptimos locales que pueden llegar a ser muy heterogéneos.

El teorema de los residuos es una parte fundamental de la variable compleja y es una generalización del teorema integral de Cauchy y la fórmula integral de Cauchy.

Fórmula integral de Cauchy

Sean \Omega \subset \mathbb{C} un abierto simplemente conexo (I(\gamma,z) = 0 si z \notin \Omega), f \in \mathcal{H}(\Omega) y \gamma un ciclo tal que \gamma^* \subset \Omega. Entonces, para z \in \Omega-\gamma^*:

\frac{1}{2 \pi i} \int_\gamma \frac{f(u)}{u-z}du = f(z) I(\gamma,z).

demostración:

Teorema integral de Cauchy.

Sean \Omega \subset \mathbb{C} un abierto simplemente conexo (I(\gamma,z) = 0 si z \notin \Omega), f \in \mathcal{H}(\Omega) y \gamma un ciclo tal que \gamma^* \subset \Omega. Entonces, para z \in \Omega-\gamma^*:

\int_\gamma f(u)du = 0

demostración:

Fórmula integral de Cauchy para las derivadas.

Sean \Omega \subset \mathbb{C} un abierto simplemente conexo (I(\gamma,z) = 0 si z \notin \Omega), f \in \mathcal{H}(\Omega) y \gamma un ciclo tal que \gamma^* \subset \Omega. Entonces, para z \in \Omega-\gamma^*:

\frac{n!}{2 \pi i} \int_\gamma \frac{f(u)}{(u-z)^{n+1}}du = f^{n)}(z) I(\gamma,z)

demostración:

Para formular matemáticamente la mecánica cuántica necesitamos hablar de espacios de Hilbert y de operadores lineales. La rama de las matemáticas que trata estos temas es el análisis funcional.

Procederemos a dar cada una de las definiciones en el momento que las necesitemos, de manera que, como lo que nos interesan son los espacios de Hilbert y los operadores lineales, definiremos estos en primer lugar. Sin embargo, como estas definiciones se basan en otras definiciones, que por abstracción adquieren identidad propia, iremos introduciendo las mismas a medida que nos vayan apareciendo.

Espacios de Hilbert

Definición (Espacio de Hilbert): Un espacio de Hilbert es un espacio prehilbertiano completo.

A esto nos referiamos, pués de momento, nos hemos quedado tal y como estabamos, pues no sabemos que significa ser prehilbertiano y tampoco completo. Estos nuevos conceptos aparecen porque al estudiar los espacios de Hilbert y abstraer ciertas partes, estas adquieren interes per se.

Definición (Espacio prehilbertiano): Un espacio prehilbertiano es un espacio vectorial H sobre un cuerpo \mathbb{K} = \mathbb{R}, \mathbb{C} en el que tenemos definido un producto escalar.

Definición (Producto escalar): Un producte escalar en un \mathbb{K}-espacio vectorial H es una aplicación:

\langle \cdot , \cdot \rangle: H \times H \longrightarrow \mathbb{K}

tal que es:

  1. lineal en la primera componente: \langle x + y , z \rangle = \langle x, z \rangle + \langle y, z \rangle y $$
  2. simétrica en \mathbb{R}\langle x,y \rangle = \langle y,x \rangle, o hermítica en \mathbb{C}, \langle x,y \rangle = \overline{ \langle y,x \rangle}.
  3. definida positiva: \langle x,x \rangle \geq 0 y \langle x,x \rangle = 0 \Leftrightarrow x=0

Ejemplos

  • (\mathbb{R}^n,\langle \cdot, \cdot \rangle) con \langle x,y \rangle = \sum_{i=1}^n x_i y_i es un espacio prehilbertiano.
  • (\mathbb{C}^n, \langle \cdot,\cdot \rangle) con \langle z,w \rangle = \sum_{i=1}^n z_i \overline{w_i} es un espacio prehilbertiano.
  • Sea A un conjunto igual a \{ 1, \ldots , n\}, \mathbb{N}, \mathbb{Z}. El espacio \ell^2(A) es un espacio prehilbertiano con el producto escalar \langle x,y \rangle := \sum_{\alpha \in A} x_\alpha \overline{y_\alpha} con x = \{ x_\alpha \}_{\alpha \in A} y y = \{ y_\alpha \}_{\alpha \in A}. En el caso que A = \mathbb{N} o \mathbb{Z} escribiremos \ell^2(A) = \ell^2. Si A = \{1,\ldots,n\} entonces \ell^2(A) \equiv \mathbb{K}^n que tambien se denota por \ell^2(n).

Vamos a demostrar este último ejemplo. Probaremos, en primer lugar, que \ell^2 es un espacio vectorial, y posteriormente, que \langle \cdot,\cdot \rangle es un producto interior.

Definición (Notación de Dirac): Una notación alternativa para el producto escalar introducida por Dirac y ampliamente utilizada en mecánica cuántica por su versatilidad es la notación bra-ket. Podemos escribir el producto escalara en notación matricial:

\langle \cdot | \cdot \rangle :

Espacios de Banach

Para empezar, recordaremos que un espacio normado es un par (E,|| \cdot ||_E) formado por un espacio vectorial E sobre un cuerpo K, que puede ser \mathbb{R} o \mathbb{C}, y una norma ||\cdot||_E sobre E, que es una función:

||\cdot||_E : E \longrightarrow [0,+\infty[

que satisface las propiedades de:

  1. separación: ||x||_E = 0 \Leftrightarrow x = 0.
  2. homogeneidad: ||\lambda x|| = |\lambda|||x||_E.
  3. desigualdad triangular:

Uno de los grandes logros de la variable compleja es su efectividad a la hora de permitir demostrar resultados muy alejadas de su aparente alcance. Un claro ejemplo de esto es su aplicación a la demostración del teorema fundamental del algebra. A continuación presentamos una demostración del mismo junto con unos resultados previos que necesitaremos para la misma.

Desigualdad de Cauchy

Sea f una función analítica, es decir, f = \sum_{n=0}^\infty a_n(z-z_0)^n con z \in D(z_0,R). Sea 0 < r < R y

M(r) := \max (|f(z)|:|z-z_0|=r).

Entonces

|a_n| \leq \frac{M(r)}{r^n}

demostración:

a_n = \frac{f^{n)}(z_0)}{n!}

ya que f'(z) = \sum_{n=1}^\infty n\,a_n (z-z_0)^{n-1}, por lo que f'(z_0) = 1 \cdot a_1. De la misma manera, f''(z) = \sum_{n=2}^\infty (n-1)n\,a_n (z-z_0)^{n-2}, por lo que f''(z_0) = 2 \cdot 1 \cdot a_2, etc.

Si aplicamos ahora la fórmula integral de Cauchy para las derivadas:

f^{n)}(z) = \frac{n!}{2 \pi i} \int_{C(z_0,R)} \frac{f(u)}{(u-z)^{n+1}} du,\, \forall z \in D(z_0,R)

tenemos:

a_n = \frac{1}{2 \pi i} \int_{C(z_0,r)} \frac{f(z)}{(z-z_0)^{n+1}} dz

por lo que si aplicamos la propiedad de que |int_\gamma f(z) dz| \leq L(\gamma) \, \max \{ |f(z)| : z \in \gamma([a,b]) \}, nos queda:

|a_n| \leq \big | \frac{1}{2 \pi i} \big | (2 \pi r) \max_{|z-z_0|=r} \big | \frac{f(z)}{(z-z_0)^{n+1}} \big |.

Como |\frac{1}{2 \pi i}| = \frac{1}{2 \pi} y

\max_{|z-z_0|=r} \big | \frac{f(z)}{(z-z_0)^{n+1}} \big | = \max_{|z-z_0|=r} \frac{|f(z)|}{r^{n+1}}

entonces:

|a_n| \leq \frac{r}{r^{n+1}} \max_{|z-z_0|=r} |f(z)| = \frac{M(r)}{r^n} \,\Box

Teorema de Liouville

Sea f \in \mathcal{H}(\mathbb{C}), f no constante. Entonces f no está acotada.

demostración:

Vamos a demostrar el reciproco, es decir, que si f\in \mathcal{H}(\mathbb{C}) y f está acotada, entonces f es necesariamente constante.

Para empezar, como f es holomorfa, entonces es analítica y

f(z) = \sum_{n=0}^\infty a_n z^n, \, \forall z \in \mathbb{C}.

Fijamos 0 < r < +\infty. Si le aplicamos ahora la desigualdad de Cauchy, nos queda

|a_n| \leq \frac{M(r)}{r^n}.

que, como f está acotada, es decir, \exists C > 0: |f(z)| \leq C, \, \forall z \in \mathbb{C}, que no depende de r, nos queda:

|a_n| \leq \frac{C}{r^n}.

Finalmente, si n \geq 1 entonces |a_n| \leq \lim_{r \rightarrow \infty} \frac{C}{r^n} = 0, por lo que f(z) = a_0 que implica que f es constante \forall z \in \mathbb{C} \, \Box

Lo que nos está diciendo el teorema de Liouville, por ejemplo, es que la función sin(z) no está acotada (a diferencia de lo que ocurre en variable real), pues solo lo estan las funciones constantes.

Teorema fundamental del Álgebra

Sea p(z) = a_0 + a_1 z + a_2 z^2 + \ldots + a_n z^n un polinomio no constante. Entonces la ecuación p(z)=0 tiene solución.

demostración:

Vamos a proceder por reducción al absurdo. Supondremos que p(z) \neq 0, \, \forall z \in \mathbb{C} y consideraremos la función f(z) = \frac{1}{p(z)}. Es fácil ver que f(z) \in \mathcal{H}(\mathbb{C}) y que no es constante.

Veremos que f(z) está acotada, lo que supondra una contradicción con Liouville. Efectivamente, para p(z) = a_0 + a_1 z + a_2 z^2 + \ldots + a_n z^n, como n \geq 1 y a_n \neq 0, tiene sentido escribir p(z) como:

p(z) = z^n (a_n + \frac{a_{n-1}}{z} + \ldots + \frac{a_0}{z^n})

con lo que, si tomamos valores absolutos y utilizando que |a+b| \geq |a|-|b|, nos queda:

|p(z)| \geq |z|^n (|a_n| - \frac{|a_{n-1}|}{|z|} - \ldots - \frac{|a_0|}{|z|^n})

\,\Box

Proposición:

Sea p(z) un polinomio de grado n \geq 1. Entonces:

p(z) = a(z-z_1)^{k_1}(z-z_2)^{k_2}\ldots (z-z_m)^{k_m} con \sum_{i=1}^m k_i = n.

demostración: (por inducción sobre n)

a) Si n=1 entonces p(z)=az+b con a \neq 0. Sea z_1 tal que p(z_1) =0. Entonces:

p(z) = az -az_1 + az_1 + b = a(z-z_1), pues az_1 + b = 0.

b) Suponemos que para grado \leq n es cierto y sea p(z) un polinomio de grado n+1.

Fórmula integral de Cauchy para las derivadas.

Sea f(z) \in \mathcal{H}(A) y \overline{D(z_0,R)} \subset A. Entonces:

f^{n)}(z) = \frac{n!}{2 \pi i} \int_{C(z_0,R)} \frac{u}{u-z}du, \, \forall z \in D(z_0,R).

demostración:

Sabemos que f(\xi) = \frac{1}{2 \pi i} \int_{C(z_0,R)} \frac{f(u)}{u-\xi}du siempre que |\xi - z_0|<R.

Fijamos ahora z \in D(z_0,R) y escogemos 0 < r < R: D(z,r) \subset D(z_0,R). Tenemos un teorema que nos asegura que si |\xi - z| < r entonces:

f(\xi) = \frac{1}{2 \pi i} \sum_{n=0}^\infty (\int_{C(z_0,R)} \frac{f(u)}{(u-z)^{n+1}} du) (\xi - z)^n,

por lo que:

\frac{f^{n)}(z)}{n!} = \frac{1}{2 \pi i} \int_{C(z_0,R)} \frac{u}{u-z}du

Aplicaciones:

1. Calcular

\int_{C(0,\frac{1}{2})} \frac{e^z}{z(1-z)^2}dz.

En este caso, la única singularidad que queda dentro de la circunferencia es el 0. Tenemos que aislarla, por lo que escribimos la integral de la siguiente forma:

\int_{C(0,\frac{1}{2})} \frac{\frac{e^z}{(1-z)^2}}{z-0} dz = \int_{C(0,\frac{1}{2})} \frac{f(z)}{(z-0)}

con f(z):=\frac{e^z}{(1-z)^2}. Ahora tenemos que f(z) es holomorfa en todos los puntos excepto en el 1, que no está dentro de la circunferencia que consideramos en la integral i, por tanto, podemos aplicar el teorema integral de Cauchy para las derivadas y obtenemos

\int_{C(0,\frac{1}{2})} \frac{f(z)}{z} = 2 \pi i f(0) = 2 \pi i

ya que f(0) = \frac{e^0}{(1-0)^2} = 1.

2. Calcular

\int_{C(1,\frac{1}{2})} \frac{e^z}{z(1-z)^2}dz.

Ahora, la singularidad que queda dentro de la circunferencia es el 1 y tenemos que aislarlo para poder aplicar el teorema integral para derivadas. Escribimos

\int_{C(1,\frac{1}{2})} \frac{\frac{e^z}{z}}{(z-1)^2} = 2 \pi i f'(1)

aplicando el teorema. De esta manera, como f'(z) = \frac{ze^z - e^z}{z^2} y f'(1)=0, donde f(z):=\frac{e^z}{z}, nos queda que la integral vale 0.

3. Calcular

I:=\int_{C(0,2)} \frac{e^z}{z(1-z)^2}dz.

Ahora las dos singularidades están dentro de la circunferencia considerada. Lo que haremos es decomponer \frac{1}{z(z-1)^2} en fracciones simples:

\frac{1}{z(z-1)^2} = \frac{A}{z} + \frac{B}{z-1} + \frac{C}{(z-1)^2}.

Esto nos lleva a un sistema cuya solución es A=1=C, B=-1, por lo que podemos escribir:

I = \int_{C(0,2)} \frac{e^z}{z} dz - \int_{C(0,2)} \frac{e^z}{z-1} dz + \int_{C(0,2)} \frac{e^z}{(z-1)^2}

de manera que podemos aplicar la fórmula integral de Cauchy para las derivadas, con f(z):=e^z que es holomorfa en todo el plano, a cada integral:

I = 2 \pi i (f(0) - f(1) + f'(1)) = 2 \pi i (1 - e + e) = 2 \pi i.

octubre 2012
L M X J V S D
« Sep   Nov »
1234567
891011121314
15161718192021
22232425262728
293031