You are currently browsing the category archive for the ‘Uncategorized’ category.

Cuando estudiaba en la universidad había toda una secuencia de asignaturas (Introducción a los Computadores, Estructura de Computadores I, Estructura de Computadores II, Arquitectura de Computadores, Multiprocesadores, Arquitecturas Vectoriales…) en las que nos enseñaban a “diseñar” un ordenador desde cero, llegando a ser tan complejo como podamos imaginar. Empezamos con circuitos combinacionales (puertas lógicas, comparadores, (semi)sumadores, (de)multiplexores, (de)codificadores, ROMs, PAL, PLA, PLDs, etc.) que se caracterizaban porque su salida solo dependía de la entrada en ese momento y que nos permitían construir, entre otras cosas, Unidades de Proceso, de las que forman parte las  Unidades Aritmético Lógicas o las Unidades de Coma Flotante (es en estas donde se ejecutan las instrucciones del lenguaje máquina de nuestro ordenador);  seguíamos con circuitos secuenciales (biestables, contadores síncronos y asíncronos, registros serie/paralelo, RAMs, etc.) cuya salida dependía además del estado en el que se encontraba un sistema cuando le llegaba la entrada (una máquina de estados) y que nos permitía saber hacer, entre otras cosas, la Unidad de Control y las memorias; y de aquí a diseño de procesadores (cableados o microprogramados), buses, memorias cache, lenguajes máquina, procesadores aritméticos, segmentación de procesadores…

Cuando diseñas lenguajes máquina, estos deben contener conjuntos de instrucciones para manejar diferentes tipos de datos, y uno de los tipos indispensables en infinidad de campos son los reales. Por una parte, esos reales van a tener que operarse en procesadores aritméticos, y por otra, los compiladores deben traducir instrucciones de lenguajes de mas alto nivel que contengan reales a instrucciones del LM. Es por toda esta interrelación entre diferentes niveles de abstracción que se hizo necesaria la especificación de un estándar (algo que todo el mundo va a cumplir si quiere continuar en el juego): el estándar IEEE 754.

Está claro que en unas aplicaciones, por ejemplo a nivel de átomos, vamos a necesitar números reales muy pequeños, mientras que en otras, si trabajamos con estrella o con galaxias, lo que necesitaremos son reales extremadamente grandes. Es por este motivo que se adopta una representación en coma flotante (notación científica donde tenemos una mantisa, una base y un exponente) frente a una representación de punto fijo mucho mas intuitiva al principio.

Si solo utilizáramos naturales, con el binario nos apañaríamos, pero en el cuanto que aparecen los enteros, necesitamos representar el signo mas y el signo menos también con ceros y unos. Ésto nos lleva a tener que reservar un bit para el signo, puesto que el número total de bits siempre va a estar determinado de antemano, y a los conceptos de representación explícita y valor representado.

Para concretar ideas, supongamos el caso sencillo de que tenemos tres bits. Los números binarios representados, los valores explícitos del sistema de representación, son del 0 al 7. Sin embargo, si estamos reservando el primer bit para el signo, donde un 0 significará + y un 1 significará -, los valores que representan esos valores explícitos son, en realidad, del +0 al +3 los cuatro primeros binarios (hasta aquí todo normal), pero ahora el valor explícito 4 representa al número entero -0, el valor explícito 5 representa al -1, el 6 se corresponde con el -2 y el 7 con el -3. Esta representación de los enteros se llama de signo y magnitud. El mayor problema que tiene, además de la duplicidad del 0, es que el algoritmo de suma es diferente del de los naturales (si hemos diseñado el sumador completo de 3 bits para naturales y ahora queremos utilizarlo para enteros, este no funciona: (+1) + (-1) tienen valores explícitos 001 y 101, que son los que vamos a meter al sumador, cuya suma da 110. Pero resulta que este valor explicito representa al -2 y no al 0…)

Existen otras representaciones de enteros, en donde lo que va a cambiar es la relación entre el valor explicito y el valor representado, que intentan resolver los diferentes problemas que nos encontramos: duplicidades del cero, algoritmos aritméticos o lógicos diferentes, rangos asimétricos, etc. Las mas conocidas son el complemento a 1 (donde el algoritmo de suma si es igual para naturales y para enteros), el complemento a 2 (donde conseguimos un solo cero a costa de perder la simetría de los rangos) y las representaciones en exceso, bias en inglés (donde los negativos tendrán un 0 delante y los positivos un 1, lo cual nos permitirá comparar enteros en comparadores binarios diseñados para naturales).

Y por fin, con todo esto, ya podemos entender las especificaciones de doble precisión del estándar IEEE 754. El estándar IEEE 754 define como representar números reales en simple y doble precisión. Como siempre se ha dicho, en un ordenador todo se almacena como ceros y unos, y el caso de los reales no podía ser menos. En simple precisión dedicaremos 32 bits (del bit 0 al bit 31) al almacenamiento del real, mientras que en doble precisión lo que tendremos serán 64 bits (bits 0, 1, …, 63). Ni uno mas. Ni uno menos… Y por cuestiones de espacio y no duplicidad, nos centraremos a partir de ahora en la doble precisión.

Para empezar, el bit 63 nos especificará el signo del número real que estamos representando.

A continuación, se reservan 11 bits para el exponente: del bit 62 al bit 52. Este exponente está en representación por exceso. Se consiguen buenas propiedades cuando este exceso es 2^{n-1}-1, que es lo que precisamente hace el estándar, de manera que 2^{11-1}-1 = 1023. Para obtener el valor explícito del valor que queremos representar, se le suma a éste el bias. Restaremos el bias para obtener el valor representado a partir del valor explicito. Con 11 bits, el rango de valores explícitos es del 0 al 2047. Por tanto, al restar el 1023 obtenemos el rango de valores representados que va de -1023 al 1024 (observar que, además, los exponentes negativos son mas pequeño que los positivos, algo bueno para utilizar comparadores binarios para comparar exponentes a la hora de hacer operaciones con números en coma flotante). Sin embargo, como el exponente mas pequeño (00000000000 que representa el -1023) y el mas grande (11111111111 que representa 1024) se reservan para casos especiales (representación del 0,  números denormalizados, mas y menos infinito y los Not a Number (NaN)), el rango de valores validos para los exponentes en el caso de números normalizados es de -1022 a 1023.

Finalmente, los últimos 52 bits, del bit 51 al 0, representan la mantisa. Recordemos que, al disponer de un exponente que varia al desplazar la coma (de ahí lo de coma flotante), la representación de un número, de entrada, no es única. Para conseguir esta unicidad, imponemos lo que se llama normalización, que quiere decir que la coma va a estar a la izquierda o la derecha de la primera cifra significativa. Esto en binario quiere decir que, dado cualquier número, siempre vamos a modificar el exponente de manera que la mantisa quede representada como 1. lo que sea o 0.1 lo que sea. Aunque la normalización a nivel general se puede hacer tanto por la izquierda como por la derecha, solo una de estas es la utilizada por el estándar, y es 1.xxx, con el primer uno a la izquierda de la coma. Además, como siempre vamos a tener este bit a 1, no lo vamos a almacenar de manera explicita y vamos a tener lo que se llama representación con bit oculto (aunque no hemos mencionado nada al respecto, lo mismo pasa con la base, que como ya sabemos que es dos, también queda almacenada de manera implícita). Con todo esto, dada una secuencia de 52 bits

Llegados a este punto, ¿para que sirven los exponentes que nos hemos reservado? Pues para tratar casos especiales. El primero es la representación del 0. Si todo número representado se tiene que interpretar como 1.xxx, esta claro que nunca vamos a conseguir, si no hacemos nada al respecto, una representación explicita del cero. Lo que hacemos es que, cuando todos los bits del exponente están a 0 (exponente -1023) y toda la mantisa también está a 0, entonces este número se interpreta como el 0. Para el mismo exponente pero mantisas diferentes de 0, tenemos los números denormalizados (se van a interpretar con bit oculto a 0, es decir, su valor sera 0.xxx con xxx siendo el valor de la mantisa en lugar de interpretarse como 1.xxx). El matlab, por ejemplo, soporta algo de esto, pues realmin nos da el mínimo numero representable en coma flotante, 2.225073858507201e-308, pero si hacemos 0.0000001/(10^308), obtenemos 9.999999984816838e-316, que es mas pequeño. La explicación es que el primero solo es el mas pequeño de los normalizados, mientras que el segundo es un denormalizado. El mas pequeño posible de los denormalizados nos lo da el comando matlab con eps(0.0). Lo mismo ocurre con el exponente mayor, que es 1024: se utiliza para representar + \infty o -\infty cuando toda la mantisa esta a cero y en función del signo; y para los NaN cuando la mantisa es diferente de 0.

Con todo esto, sea s el bit de signo, e los bits de exponente y m los de mantisa. El valor real es -(1)^s (2^{e-1023})(1.m) si e \in (0,2047). Si e=0 y m \neq 0 entonces el valor es -(1)^s (2^{e-1023})(0.m). Si e=0 y m=0 entonces el valor es 0. Todo lo demás, e = 2047, son casos que no representan números sino casos especiales (infinitos, NaN…). Y fuera de aquí tenemos overflows i underflows.

Y ya para acabar. En general , y a la luz del estándar, para que un sistema en coma flotante en un computador quede determinado, necesitamos especificar, como mínimo, el número total de bits que se utilizarán (uno será el signo), cuantos serán para la mantisa, cuantos se utilizarán para el exponente y cual será el exceso de su representación, como representaremos el cero y cual es el criterio de normalización.

 

En el artículo técnico What would a binary black hole merger look like?  se simula, mediante técnicas de ray tracing, como vería un observador externo el merge de dos agujeros negros (del artículo y del sitio web de los autores están sacadas prácticamente todas las imágenes).

Si estamos mirando hacia algún sitio, digamos:

ClockTower-400x300

y pasa un agujero negro frente a nosotros, lo primero que se nos viene a la cabeza es la siguiente imagen:

ClockTower-400x300b

ya que como de un agujero negro no puede escapar nada, ni la luz, pensamos que veríamos una simple esfera negra tapando un trozo de nuestra visión. Sin embargo, una imagen mas realista de lo que veríamos es:

ClockTower_BH-400x300debido a la curvatura que experimentan los rayos de luz por la curvatura del espacio-tiempo que genera el agujero negro: efecto de lente gravitacional.

Colocando una imagen de fondo más métrica, así es como se verían los espacios de Minkowski, Schwarzschild y Kerr.

analyticSpacetimesSi en lugar de un agujero negro tenemos un sistema binario de agujeros negros de igual masa, entonces tendríamos:

bbhSystem

Finalmente, una animación del merge:

Acabo de descubrir, en esta entrada de Tao, unos vídeos sobre los medallistas Fields del 2014.

Añado también una imagen curiosa sobre los medallistas de la IMO del 1995, donde lograron medalla de oro tanto Mirzakhani (7) como Avila (24). También aparece con medalla de plata Ben J. Green (46):

imo1995

Finalmente, un software creado por Hairer para Mac.

Dada una función f:[a,b] \rightarrow \mathbb{R} tal que f''(x)>0 entonces alcanza su valor máximo en uno de los dos extremos. La segunda derivada f''(x) no esta dando la pendiente de la función f'(x) en cada punto. Como ésta es positiva, nos indica que la derivada de f es estrictamente creciente, es decir, que la función no tiene ningún extremo relativo en el interior del intervalo.

De esta manera, diremos que una función que cumple f''>0 en un dominio A y que por este motivo alcanza su máximo en la frontera \partial A tiene un principio del máximo.

En una dimensión tenemos que si f:[a,b] \rightarrow \mathbb{R} continua y suficientemente derivable, entonces c \in (a,b) es un máximo relativo si u'(x) = 0 y u''(x) \leq 0. Si consideramos una función u tal que L:= u'' + g(x)u' >0, con g acotada, es imposible que cumpla las condiciones de máximo anteriores, por lo que se alcanzaran en los extremos.

Recordemos lo ya expuesto en este post: que en las coordenadas esferoidales prolatas (\mu, \nu, \varphi), las dos primeras (\mu, \nu) provienen de las coordenadas elípticas, donde \mu \in ]0,+\infty[ y \nu \in ]0,2\pi[, mientras que la última \varphi \in ]0,2\pi[ proviene de rotarlas alrededor del eje que une los focos.

Compactificamos la primera coordenada mediante \boxed{\mu = b \tan \frac{\pi \bar{\mu}}{2}}.

El Laplaciano y las fuentes, en estas coordenadas y con esta compactificación, utilizando una nueva función en Mathematica que nos lo calcula todo, quedan:

lap_ellComNor2

\boxed{\Delta \Theta_{X} = 6 \pi \mathcal{D}^j S^*_j}

s1_ellComNor2

\boxed{\Delta X^{i} = 8 \pi f^{ij} S^*_j - \frac{1}{3} \mathcal{D}^i \Theta_X}

s21_ellComNor2

\underline{\hat{A}^{ij} = \mathcal{D}^i X^j + \mathcal{D}^j X^i - \frac{2}{3} \mathcal{D}_k X^k f^{ij}}

A1x_ellComNor2

A2x_ellComNor2

A3x_ellComNor2

\boxed{\Delta \psi = -2 \pi E^* \psi^{-1} - \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-7} }

\boxed{\Delta (\alpha \psi) = [ 2 \pi (E^* + 2 S^*) \psi^{-7} + \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-8} ] (\alpha \psi) }

\boxed{\Delta \Theta_{\beta} = \frac{3}{4} \mathcal{D}_i \mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij} )}

\boxed{\Delta \beta^i = \mathcal{D}_j ( 2 \alpha \psi^{-6} \hat{A}^{ij} ) - \frac{1}{3} \mathcal{D}^i \Theta_{\beta} }

Recordemos lo ya expuesto en este post: que en las coordenadas biesféricas (\xi, \eta, \varphi), las dos primeras (\xi, \eta) provienen de las coordenadas bipolares, donde la primera indica el ángulo entre las dos rectas que unen nuestro punto con los dos focos que necesitamos para determinar las bipolares y la segundo es el logartimo del ratio entre la longitud de estas dos rectas, mientras que la última proviene de rotarlas alrededor del eje que une los focos.

Compactificamos la segunda coordenada mediante \boxed{\eta = \frac{b \bar{\eta}}{1 - \bar{\eta}}}.

El Laplaciano, en estas coordenadas y con esta compactificación, queda:

\Delta = \frac{(\cos \xi - \mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1-\bar{\eta}})^2}{a^2} \big [ \partial_{\xi \xi} + \csc \xi \frac{-1 + \cos \xi \, \mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1-\bar{\eta}}}{\mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1-\bar{\eta}} - \cos \xi} \partial_{\xi}

+\frac{(\bar{\eta} - 1)^4}{b^2} \partial_{\bar{\eta} \bar{\eta}} + \frac{(\bar{\eta} - 1)^2}{b^2} (2(\bar{\eta}-1) -\frac{b \, \mbox{\scriptsize sinh} \frac{b \bar{\eta}}{1 - \bar{\eta}}}{\mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1 - \bar{\eta}} - \cos \xi}) \partial_{\bar{\eta}} + \csc^2 \xi \partial_{\varphi} \big ],

las derivadas covariantes de covectores (1-formas):

CovDer_BiSphComNor1

y las fuentes:

\boxed{\Delta \Theta_{X} = 6 \pi \mathcal{D}^j S^*_j}

\Delta \Theta_X = 6 \pi f^{ji} \mathcal{D}_i S^*_j = 6 \pi ( \mathcal{D}_{\xi} S^*_{\xi} + \mathcal{D}_{\bar{\eta}} S^*_{\bar{\eta}} + \mathcal{D}_{\varphi} S^*_{\varphi} ) =

s1_biSphComNor1

\boxed{\Delta X^{i} = 8 \pi f^{ij} S^*_j - \frac{1}{3} \mathcal{D}^i \Theta_X}

Pasando la derivada contravariante a covariante mediante la métrica, queda:

\Delta X^{i} = 8 \pi f^{ij} S^*_j - \frac{1}{3} f^{ik} \mathcal{D}_k \Theta_X.

Definimos ahora

S_X^i := 8 \pi f^{ij} S^*_j - \frac{1}{3} f^{ik} \mathcal{D}_k \Theta_X,

de manera que:

S_X^{\xi} = 8 \pi f^{\xi j} S^*_j - \frac{1}{3} f^{\xi k}\mathcal{D}_{k} \Theta_X = 8 \pi S^*_{\xi} - \frac{1}{3} \mathcal{D}_{\xi} \Theta_X =

= 8 \pi S^*_{\xi} - \frac{\mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1 - \bar{\eta}} - \cos \xi}{a} \partial_{\xi} \Theta_X

S_X^{\bar{\eta}} = 8 \pi f^{\bar{\eta} j} S^*_j - \frac{1}{3} f^{\bar{\eta} k} \mathcal{D}_{k} \Theta_X = 8 \pi S^*_{\bar{\eta}} - \frac{1}{3} \mathcal{D}_{\bar{\eta}} \Theta_X =

= 8 \pi S^*_{\bar{\eta}} - \frac{\mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1 - \bar{\eta}} - \cos \xi}{a} \frac{(\bar{\eta} - 1)^2}{b} \partial_{\bar{\eta}} \Theta_X

S_X^{\varphi} = 8 \pi f^{\varphi j} S^*_j - \frac{1}{3} f^{\varphi k} \mathcal{D}^{k} \Theta_X = 8 \pi S^*_{\varphi} - \frac{1}{3} \mathcal{D}_{\varphi} \Theta_X =

= 8 \pi S^*_{\varphi} - \frac{\mbox{\scriptsize cosh} \frac{b \bar{\eta}}{1 - \bar{\eta}} - \cos \xi}{a} \csc \xi \partial_{\varphi} \Theta_X

En este punto tenemos que el vector

(S_X^{\xi}(\xi,\bar{\eta},\varphi),S_X^{\bar{\eta}}(\xi,\bar{\eta},\varphi),S_X^{\varphi}(\xi,\bar{\eta},\varphi))

expresado en la base que resulta de normalizar la base coordenada \{ \partial_{\xi}, \partial_{\bar{\eta}}, \partial_{\varphi} \}. Lo que hacemos ahora es expresar este vector en la nueva base \{ \partial_x, \partial_y, \partial_z \}, de manera que obtenemos

(S_X^{x}(\xi,\bar{\eta},\varphi),S_X^{y}(\xi,\bar{\eta},\varphi),S_X^{z}(\xi,\bar{\eta},\varphi)).

y como es esta base las ecuaciones están desacopladas y \Theta_X es un campo escalar, resolvemos independientemente:

\Delta X^{x} = S_X^{x},

\Delta X^{y} = S_X^{y},

\Delta X^{z} = S_X^{z}.

Finalmente, con el cambio de base inverso, calculamos a partir de (X^{x},X^{y},X^{z}) el vector (X^{\xi},X^{\bar{\eta}},X^\varphi) .

\underline{\hat{A}^{ij} = \mathcal{D}^i X^j + \mathcal{D}^j X^i - \frac{2}{3} \mathcal{D}_k X^k f^{ij}}

Necesitamos ahora la derivada covariante de un vector (hasta ahora habían coincidido las derivadas covariantes de vectores y covectores, pero en este caso no):

CovDer_BiSphComNor1_vec

volvemos a pasar las derivadas contravariantes a covariantes:

\hat{A}^{ij} = f^{im} \mathcal{D}_m X^j + f^{jn} \mathcal{D}_n X^i - \frac{2}{3} f^{ij} \mathcal{D}_k X^{k}

y obtenemos:

\hat{A}^{\xi \xi} = f^{\xi m} \mathcal{D}_m X^{\xi} + f^{\xi n} \mathcal{D}_n X^{\xi} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( 2 \mathcal{D}_{\xi} X^{\xi} - \mathcal{D}_{\bar{\eta}} X^{\bar{\eta}} - \mathcal{D}_{\varphi} X^{\varphi}) =

A11_biSphComNor1

\hat{A}^{\xi \bar{\eta}} = f^{\xi m} \mathcal{D}_m X^{\bar{\eta}} + f^{\bar{\eta} n} \mathcal{D}_n X^{\xi} = \mathcal{D}_{\xi} X^{\bar{\eta}} + \mathcal{D}_{\bar{\eta}} X^{\xi} =

A12_biSphComNor1

\hat{A}^{\xi \varphi} = f^{\xi m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\xi} = \mathcal{D}_{\xi} X^{\varphi} + \mathcal{D}_{\varphi} X^{\xi} =

A13_biSphComNor1

\hat{A}^{\bar{\eta} \bar{\eta}} = f^{\bar{\eta} m} \mathcal{D}_m X^{\bar{\eta}} + f^{\bar{\eta} n} \mathcal{D}_n X^{\bar{\eta}} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( - \mathcal{D}_{\xi} X^{\xi} + 2 \mathcal{D}_{\bar{\eta}} X^{\bar{\eta}} - \mathcal{D}_{\varphi} X^{\varphi}) =

A22_biSphComNor1

\hat{A}^{\bar{\eta} \varphi} = f^{\bar{\eta} m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\bar{\eta}} = \mathcal{D}_{\bar{\eta}} X^{\varphi} + \mathcal{D}_{\varphi} X^{\bar{\eta}} =

A23_biSphComNor1

\hat{A}^{\varphi \varphi} = f^{\varphi m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\varphi} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( - \mathcal{D}_{\bar{r}} X^{\bar{r}} - \mathcal{D}_{\theta} X^{\theta} +2 \mathcal{D}_{\varphi} X^{\varphi}) =

A33_biSphComNor1

Las dos ecuaciones no lineales correspondientes al factor conforme \psi y al lapse \alpha, como no contienen derivadas covariantes, quedan como las teniamos:

\boxed{\Delta \psi = -2 \pi E^* \psi^{-1} - \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-7} }

\boxed{\Delta (\alpha \psi) = [ 2 \pi (E^* + 2 S^*) \psi^{-7} + \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-8} ] (\alpha \psi) }

Finalmente, para el shift \beta y su ecuación auxiliar tenemos:

\boxed{\Delta \Theta_{\beta} = \frac{3}{4} \mathcal{D}_i \mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij} )}

\boxed{\Delta \beta^i = \mathcal{D}_j ( 2 \alpha \psi^{-6} \hat{A}^{ij} ) - \frac{1}{3} \mathcal{D}^i \Theta_{\beta} }

que trataremos en el siguiente post.

El círculo de Apolonio de focos F_1 y F_2 de razón r es el lugar geométrico de los puntos de plano P tales que:

\overline{PF_1}/\overline{PF_2} = r.

Definimos las coordenadas bipolares (\xi,\eta) de la siguiente manera:

x = a \frac{\mbox{\scriptsize sinh} \eta}{\mbox{\scriptsize cosh} \eta - cos \xi},

y = a \frac{\sin \eta}{\mbox{\scriptsize cosh} \eta - cos \xi},

donde en los puntos F_1=(-a,0) y F_2=(a,0) tenemos definidos los dos focos.

La interpretación geométrica es la siguiente: dado un punto P=(x,y) del plano, la coordenada \xi es el ángulo que forman los vectores u:=\overline{PF_1} y v:=\overline{PF_2}:

\xi = \arccos \frac{u \cdot v}{|u||v|}, de manera que \xi \in [-\pi,\pi],

mientras que la coordenada \eta es el rátio entre los módulos de éstos:

\eta = \log \frac{|u|}{|v|}, por lo que \eta \in ]0,\infty[,

que si las expresamos de manera que aparezcan explicitamente las coordenadas cartesianas, quedan:

\xi = \arccos \frac{x^2 + y^2 - a^2}{\sqrt{(a-x)^2 + y^2}\sqrt{(a+x)^2+y^2}},

\eta = \frac{\log[(a+x)^2+y^2] - \log[(a-x)^2+y^2]}{2}.

A continuación, una imagen con lineas coordenadas:

BiSphCoo

Para pasar a coordenadas en tres dimensiones, lo único que hacemos es añadir una nueva coordenada \varphi que nos indica un ángulo de rotación respecto de un eje. Si el eje es la recta que une los focos, obtenemos las coordenadas biesféricas, que son las que nos interesarán y utilizaremos en posteriores entradas. Si lo que hacemos es rotar respecto a la recta perpendicular a la anterior, la que separa los dos focos, obtenemos las coordenadas toroidales.

Compactificamos la primera coordenada mediante \boxed{r = \frac{a \bar{r}}{1 - \bar{r}}}.

El Laplaciano, con esta compactificación, queda:

\Delta = \frac{(a-\bar{r})^4}{a^4} \partial_{\bar{r} \bar{r}} + \frac{2}{\bar{r}} \frac{(a-\bar{r})^4}{a^4} \partial_{\bar{r}} + \frac{(a - \bar{r})^2}{a^2 \bar{r}^2} \partial_{\theta \theta} + \frac{(a-\bar{r})^2}{a^2 \bar{r}^2} \cot \theta \partial_{\theta} + \frac{(a - \bar{r})^2}{a^2 \bar{r}^2 } \csc \theta \partial_{\varphi \varphi}

y las fuentes:

\boxed{\Delta \Theta_{X} = 6 \pi \mathcal{D}^j S^*_j}

\Delta \Theta_X = 6 \pi f^{ji} \mathcal{D}_i S^*_j = 6 \pi ( \mathcal{D}_{\bar{r}} S^*_{\bar{r}} + \mathcal{D}_{\theta} S^*_{\theta} + \mathcal{D}_{\varphi} S^*_{\varphi} ) =

= 6 \pi \frac{1-\bar{r}}{a \bar{r}} ( (\bar{r}-\bar{r}^2) \partial_{\bar{r}} S^*_{\bar{r}} + 2 S^*_{\bar{r}} + \partial_{\theta} S^*_{\theta} + \cot \theta S^*_{\theta} + \csc \theta \partial_{\varphi} S^*_{\varphi} )

\boxed{\Delta X^{i} = 8 \pi f^{ij} S^*_j - \frac{1}{3} \mathcal{D}^i \Theta_X}

Pasando la derivada contravariante a covariante mediante la métrica, queda:

\Delta X^{i} = 8 \pi f^{ij} S^*_j - \frac{1}{3} f^{ik} \mathcal{D}_k \Theta_X.

Definimos ahora

S_X^i := 8 \pi f^{ij} S^*_j - \frac{1}{3} f^{ik} \mathcal{D}_k \Theta_X,

de manera que:

S_X^{\bar{r}} = 8 \pi f^{\bar{r} j} S^*_j - \frac{1}{3} f^{\bar{r} k}\mathcal{D}_{k} \Theta_X = 8 \pi S^*_{\bar{r}} - \frac{1}{3} \mathcal{D}_{\bar{r}} \Theta_X = 8 \pi S^*_{\bar{r}} - \frac{(1-\bar{r})^2}{a^2} \partial_{\bar{r}} \Theta_X,

S_X^{\theta} = 8 \pi f^{\theta j} S^*_j - \frac{1}{3} f^{\theta k} \mathcal{D}_{k} \Theta_X = 8 \pi S^*_{\theta} - \frac{1}{3} \mathcal{D}_{\theta} \Theta_X = 8 \pi S^*_{\theta} - \frac{1-\bar{r}}{a \bar{r}} \partial_{\theta} \Theta_X,

S_X^{\varphi} = 8 \pi f^{\varphi j} S^*_j - \frac{1}{3} f^{\varphi k} \mathcal{D}^{k} \Theta_X = 8 \pi S^*_{\varphi} - \frac{1}{3} \mathcal{D}_{\varphi} \Theta_X = 8 \pi S^*_{\varphi} - \frac{a-\bar{r}}{a \bar{r}} \csc \theta \partial_{\varphi} \Theta_X.

En este punto tenemos que el vector

(S_X^{\bar{r}}(\bar{r},\theta,\varphi),S_X^{\theta}(\bar{r},\theta,\varphi),S_X^{\varphi}(\bar{r},\partial_\theta,\partial_\varphi))

expresado en la base \{ \partial_{\bar{r}}, \theta, \varphi \}. Lo que hacemos ahora es expresar este vector en la nueva base \{ \partial_x, \partial_y, \partial_z \}, de manera que obtenemos

(S_X^{x}(\bar{r},\theta,\varphi),S_X^{y}(\bar{r},\theta,\varphi),S_X^{z}(\bar{r},\theta,\varphi)).

y como es esta base las ecuaciones están desacopladas y \Theta_X es un campo escalar, resolvemos independientemente:

\Delta X^{x} = S_X^{x},

\Delta X^{y} = S_X^{y},

\Delta X^{z} = S_X^{z}.

Finalmente, con el cambio de base inverso, calculamos a partir de (X^{x},X^{y},X^{z}) el vector (X^{\bar{r}},X^\theta,X^\varphi) .

\underline{\hat{A}^{ij} = \mathcal{D}^i X^j + \mathcal{D}^j X^i - \frac{2}{3} \mathcal{D}_k X^k f^{ij}}

volvemos a pasar las derivadas contravariantes a covariantes:

\hat{A}^{ij} = f^{im} \mathcal{D}_m X^j + f^{jn} \mathcal{D}_n X^i - \frac{2}{3} f^{ij} \mathcal{D}_k X^{k}

y obtenemos:

\hat{A}^{\bar{r} \bar{r}} = f^{\bar{r} m} \mathcal{D}_m X^{\bar{r}} + f^{\bar{r} n} \mathcal{D}_n X^{\bar{r}} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( 2 \mathcal{D}_{\bar{r}} X^{\bar{r}} - \mathcal{D}_{\theta} X^{\theta} - \mathcal{D}_{\varphi} X^{\varphi}) =

= \frac{2(1-\bar{r})}{3a\bar{r}} [ 2(\bar{r}-\bar{r}^2)\partial_{\bar{r}} X^{\bar{r}} - 2 X^{\bar{r}} - \partial_{\theta} X^{\theta} - \cot \theta X^{\theta} - \csc \theta \partial_{\varphi} X^{\varphi} ]

\hat{A}^{\bar{r} \theta} = f^{\bar{r} m} \mathcal{D}_m X^{\theta} + f^{\theta n} \mathcal{D}_n X^{\bar{r}} = \mathcal{D}_{\bar{r}} X^{\theta} + \mathcal{D}_{\theta} X^{\bar{r}} =

= \frac{1-\bar{r}}{a\bar{r}} [ (\bar{r}-\bar{r}^2) \partial_{\bar{r}} X^{\theta} - X^{\theta} + \partial_{\theta} X^{\bar{r}} ],

\hat{A}^{\bar{r} \varphi} = f^{\bar{r} m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\bar{r}} = \mathcal{D}_{\bar{r}} X^{\varphi} + \mathcal{D}_{\varphi} X^{\bar{r}} =

= \frac{1-\bar{r}}{a\bar{r}} [ (\bar{r}-\bar{r}^2) \partial_{\bar{r}} X^{\varphi} - X^{\varphi} + \csc \theta \partial_{\varphi} X^{\bar{r}} ],

\hat{A}^{\theta \theta} = f^{\theta m} \mathcal{D}_m X^{\theta} + f^{\theta n} \mathcal{D}_n X^{\theta} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( - \mathcal{D}_{\bar{r}} X^{\bar{r}} + 2 \mathcal{D}_{\theta} X^{\theta} - \mathcal{D}_{\varphi} X^{\varphi}) =

= \frac{2(1-\bar{r})}{3 a \bar{r}} [ -(\bar{r}-\bar{r}^2) \partial_{\bar{r}} X^{\bar{r}} + X^{\bar{r}} + 2 \partial_{\theta} X^{\theta} - \cot \theta X^{\theta} - \csc \partial_{\varphi} X^{\varphi} ]

\hat{A}^{\theta \varphi} = f^{\theta m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\theta} = \mathcal{D}_{\theta} X^{\varphi} + \mathcal{D}_{\varphi} X^{\theta} =

= \frac{1-\bar{r}}{a\bar{r}} [ \partial_{\theta} X^{\varphi} - \cot \theta X^{\varphi} + \csc \theta \partial_{\varphi} X^{\theta} ],

\hat{A}^{\varphi \varphi} = f^{\varphi m} \mathcal{D}_m X^{\varphi} + f^{\varphi n} \mathcal{D}_n X^{\varphi} - \frac{2}{3} \mathcal{D}_k X^{k} = \frac{2}{3}( - \mathcal{D}_{\bar{r}} X^{\bar{r}} - \mathcal{D}_{\theta} X^{\theta} +2 \mathcal{D}_{\varphi} X^{\varphi}) =

= \frac{2(1-\bar{r})}{3 a \bar{r}} [ -(\bar{r} - \bar{r}^2) \partial_{\bar{r}} X^{\bar{r}} + X^{\bar{r}} - \partial_{\theta} X^{\theta} + 2 \cot \theta X^{\theta} + 2 \csc \theta \partial_{\varphi} T^{\varphi} ]

Las dos ecuaciones no lineales correspondientes al factor conforme \psi y al lapse \alpha, como no contienen derivadas covariantes, quedan como las teniamos:

\boxed{\Delta \psi = -2 \pi E^* \psi^{-1} - \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-7} }

\boxed{\Delta (\alpha \psi) = [ 2 \pi (E^* + 2 S^*) \psi^{-7} + \frac{1}{8}(f_{il} f_{jm} \hat{A}^{lm} \hat{A}^{ij}) \psi^{-8} ] (\alpha \psi) }

Finalmente, para el shift \beta y su ecuación auxiliar tenemos:

\boxed{\Delta \Theta_{\beta} = \frac{3}{4} \mathcal{D}_i \mathcal{D}_j (2 \alpha \psi^{-6} \hat{A}^{ij} )}

\boxed{\Delta \beta^i = \mathcal{D}_j ( 2 \alpha \psi^{-6} \hat{A}^{ij} ) - \frac{1}{3} \mathcal{D}^i \Theta_{\beta} }

que lo tratamos en este post.

No hay que confundir los dos conceptos del título del post. Por un lado tenemos claramente definido el concepto de serie convergente: una serie es convergente si la sucesión de las sumas parciales tiene límite. Por otro, existen métodos de sumación, que pretenden extender, de manera consistente, la asignación de un valor de suma a una serie divergente en el sentido anterior.

Sin entrar en criterios de convergencia de series, existe una manera muy sencilla de saber si una serie diverge: si la sucesión de términos no tiende a cero, entonces la serie diverge:

\lim_{n \rightarrow \infty} a_n \neq 0 \Rightarrow \sum_{n=1}^{\infty} a_n = \infty.

Por ejemplo, las series:

\sum_{n=1}^{\infty} (-1)^{n-1}n = 1 - 2 + 3 - \ldots

\sum_{n=1}^{\infty} n = 1 + 2 + 3 + \ldots

son divergentes, ya que la sucesión de sus correspondientes sumas parciales (1,-1,2,-2\ldots; 1,3,6,10,\ldots ) no tiene límite o, por el criterio anterior, \lim_{n \rightarrow \infty} (-1)^{n-1} n \neq 0 y \lim_{n \rightarrow \infty} n \neq 0.

Sin embargo, podemos asignarle un valor a esta suma (jugando apropiadamente con los términos :-)). En el primer caso, ya Euler encontró que:

\sum_{n=1}^{\infty} n = 1 - 2 +3 -4 + \ldots = \frac{1}{4},

y posteriormente se establecieron métodos bien definidos para encontrar sumas generalizadas de series divergentes.

En el segundo caso, podemos utilizar la extensión analítica (de funciones analíticas y sus extensiones comenté algo en este post) de la función zeta de Riemann (una nota curiosa aquí):

\zeta(s):=\sum_{n=1}^{\infty} \frac{1}{n^s} = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \ldots,

ya que nuestra serie no es mas que \zeta(-1) y, como se puede demostrar que:

\zeta(-s) = \frac{B_{s+1}}{s+1},

donde B_n son los números de Bernoulli, entonces tenemos que:

\sum_{n=1}^{\infty} n = 1 + 2 + 3 + 4 + \ldots = -\frac{1}{12},

En este post, Tao da una interpretación matemáticamente consistente de lo que son estos valores.

Ya vimos aquí que a partir de la parametrización de la esfera S(r):

\Phi(\theta,\varphi) = r(\sin \theta \cos \varphi,\sin \theta \sin \varphi, cos \theta) con \theta \in [0,\pi] y \varphi \in [0,2\pi],

obtenemos la métrica:

g = r^2 d\theta \otimes d\theta + r^2 \sin^2 \theta d\varphi \otimes d\varphi,

donde la base para el espacio cotangente es \{ d\theta, d\varphi\} y para el tangente es \{ \partial_\theta, \partial_\varphi\}.

Supongamos que, como hicimos post, expresamos el cambio a esféricas mediante una carta, de manera que:

g_{ij} = dr \otimes dr + r^2 d\theta \otimes d\theta + r^2 \sin^2 \theta d\varphi \otimes d\varphi.

Ya vimos que el espacio tangente es un espacio vectorial, de manera que podemos definir diferentes bases para el mismo. Es por tanto normal preguntarse cómo varia la métrica al realizar estos cambios.

Supongamos que a la base coordenada \{ \partial_r, \partial_\theta, \partial_\varphi \} la llamamos e_i = \{e_1, e_2, e_3\}, y su base dual es \{ dr, d\theta, d\varphi\}. Supongamos dos nuevas bases:

\hat{e}_i := \{ \hat{e}_1, \hat{e}_2, \hat{e}_3\} = \{ \partial_r, \frac{1}{r} \partial_\theta, \frac{1}{r^2 \sin \theta} \partial_\varphi\},

\tilde{e}_i := \{ \tilde{e}_1, \tilde{e}_2, \tilde{e}_3\} = \{ \partial_r + \partial_\theta, \partial_r - \partial_\theta, -\partial_\varphi\}.

Como:

\hat{e}_1 = e_1, \hat{e}_2 = \frac{1}{r} e_2 y \hat{e}_3 = \frac{1}{r \sin \theta} e_3,

tenemos que:

A_{\hat{i}}^{i} = \left(  \begin{array}{ccc}  1 & 0 & 0 \\  0 & \frac{1}{r} & 0 \\  0 & 0 & \frac{\text{Csc}[\theta ]}{r}  \end{array}  \right) i A_{i}^{\hat{i}} = \left( \begin{array}{ccc}  1 & 0 & 0 \\  0 & r & 0 \\  0 & 0 & r \text{Sin}[\theta ]  \end{array}  \right),

de manera que, como la métrica es un tensor, cambia de base como estos:

g_{\hat{i} \hat{j}} = g_{ij} A_{\hat{i}}^{i} A_{\hat{i}}^{i} = d\hat{r} \otimes d\hat{r} + d\hat{\theta} \otimes d\hat{\theta} + d\hat{\varphi} \otimes d\hat{\varphi} = \eta_{ab}.

De la misma manera, tenemos:

\tilde{e}_1 = e_1 + e_2, \tilde{e}_2 = e_1-e_2 y \tilde{e}_3 = -e_3,

por lo que:

 A_{\tilde{i}}^{i} = \left(  \begin{array}{ccc}  1 & 1 & 0 \\  1 & -1 & 0 \\  0 & 0 & -1  \end{array}  \right) i A_{i}^{\tilde{i}} = \left(  \begin{array}{ccc}  \frac{1}{2} & \frac{1}{2} & 0 \\  \frac{1}{2} & -\frac{1}{2} & 0 \\  0 & 0 & -1  \end{array}  \right)

y ahora:

g_{\tilde{i} \tilde{j}} = g_{ij} A_{\tilde{i}}^{i} A_{\tilde{i}}^{i} = 2 d\tilde{r} \otimes d\tilde{r} + 2r^2 d\tilde{\theta} \otimes d\tilde{\theta} + r^2 \sin^2 \theta d\varphi \otimes d\varphi.

Seguimos con el post anterior. Vamos a utilizar ahora siete puntos separados de la siguiente manera:

x_0

x_1:=x_0+n^2 h

x_2:=x_0 + (n^2+n) h

x_3:=x_0 + (n^2+n+1) h

x_4:=x_0 + (n^2+n+2) h

x_5:=x_0 + (n^2+2n+2) h

x_6:=x_0 + (2n^2+2n+2) h

Dado cualquier número primo p, el conjunto \mathbb{Z}/p\mathbb{Z} de los enteros módulo p tiene estructura de cuerpo. Solemos escribir \mathbb{F}_p para referirnos al cuerpo finito de p elementos.

En particular, cuando p=2, tenemos el cuerpo \mathbb{F}_2 con la suma y el producto, que están en correspondencia, desde el punto de vista de las álgebras de Boole, con las operaciones XOR y AND.

Si \mathbb{K} es un cuerpo, se puede dotar fácilmente a \mathbb{K}^n de estructura de cuerpo.

Por lo tanto, (\mathbb{F}_2)^n es un cuerpo, también finito. Como todos los cuerpos finitos de q elementos son isomorfos entre ellos y el cuerpo anterior tiene 2^n, estamos tratando con \mathbb{F}_{2^n}. Podriamos haber construido el cuerpo de otra manera: buscar un polinomio irreducible f(x) de grado n con coeficientes en \mathbb{F}_2 y construir F_{2^n} como:

\mathbb{F}_{2^n} = \frac{\mathbb{F}_2[x]}{\langle f(x) \rangle}.

Se puede demostrar que cualquier cuerpo finito \mathbb{F}_p es un espacio vectorial sobre cierto \mathbb{Z}/q\mathbb{Z}. En particular, \mathbb{F}_{2^n} es un espacio vectorial sobre \mathbb{F}_2 de dimensión n.

Existe un teorema, que en nuestro nos lleva a algo que ya sabiamos, que nos dice que si V es un \mathbb{K}-ev de dimensión n entonces V es isomorfo a \mathbb{F}^n (en nuestro caso V es \mathbb{F}_{2^n} y \mathbb{K} es \mathbb{F}_2, por lo que \mathbb{F}_{2^n} \cong (\mathbb{F}_2)^n).

Dado un anillo R y sea A = R[x] el anilllo de polinomios sobre R. Entonces la derivada formal de f(x)=a_nx^n + \ldots + a_1 x + a_0 \in R[x] es f'(x) = Df(x) = n a_n x^{n-1} + \ldots + a_1.

 

 

Mi hermano, que es un manitas, se cambió la pantalla de su iPhone en un rato. Yo, que de manitas no tengo nada, me animé a hacer lo mismo. A continuación, una galería de fotos de mi experiencia:

No tiene nada que ver con los temas habituales del blog, pero bueno, la vida es algo mas que ciencia y tecnología ¿no?. Simplemente unas frases que me gustaron de “El lado bueno de las cosas”:

“Estoy convencido de ello. Tienes que hacer todo lo que puedes y esforzarte al máximo y, si mantienes el optimismo, siempre te quedará el lado bueno de las cosas.”

“Te lo advierto. Tienes que estar atento a las señales. Cuando la vida te brinda un momento como éste, es un pecado no aprovecharlo. Es un pecado.”

Interesantes 🙂

Estas Navidades me han regalado un muy interesante y entretenido libro: Física de lo imposible de Michio Kaku.

Al autor, el físico teórico Kaku, ya lo conocía con anterioridad de su libro Hiperespacio y de algunos documentales sobre ciencia. No se donde leí que sería el equivalente a nuestro Eduard Punset en cuanto a divulgación científica. En el  blog Francis (th)E mule Science’s News aparece en el post Los héroes de la ciencia del siglo XX como muñecos de Action Man dibujados. Además, en el libro cuenta algo de su vida:

“Para mi proyecto de ciencias en el instituto monté un colisionador de átomos en el garaje de mi madre.” (casi nada…) “Fui a la compañia Westing-house y reuní 200 kilos de chatarra procedente de un transformador. Durante las navidades bobiné 35 kilómetros” (vaya con el chaval…) “de cable de cobre en el campo de fútbol del instituto. Finalmente construí un betatrón de 2,5 millones de electones-voltio que consumía 6 kilovatios (toda la potencia eléctrica de mi casa) y generaba un campo magnético 20.000 veces mayor que el campo magnético de la Tierra. El objetivo era generar un haz de rayos gamma suficientemente potente para crear antimateria.”

“Mi primer contacto con Teller se remonta a la época en que yo estaba en el instituto Realicé una serie de experimentos sobre la naturaleza de la antimateria y gané el primer premio en la feria de la ciencia de San Francisco y un viaje a la Feria Nacional de la Ciencia en Albuquerque, Nuevo México. Aparecí en la televisión local con Teller, que estaba interesado en físicos jóvenes y brillantes. Con el tiempo, se me concedió la beca de ingeniería Hertz de Teller, que costeó mi educación universitaria en Harvard. Llegué a conocer bastante bien a su familia después de varias visitas a su casa en Berkeley.”

Con respecto al libro, aborda temas habituales en la ciencia ficción clasificándolos en

  • Imposibilidades de clase I, en palabras del propio Kaku: “Son tecnologías que hoy son imposibles pero que no violan las leyes de la física conocidas. Por ello, podrían ser posibles en este siglo, o en el próximo, de forma modificada”. Comenta: campos de fuerza, invisibilidad, fáseres y estrellas de la muerte, teletransporte, telepatía, psicoquinesia, robots, extraterrestres y ovnis, naves estelares, antimateria y antiuniversos.
  • Imposibilidades de clase II: “Son tecnologías situadas en el límite de nuestra comprensión del mundo físico”. Escribe sobre: más rápido que la luz, el viaje en el tiempo, universos paralelos.
  • Imposibilidades de clase III: “Son tecnologías que violan las leyes de la física conocidas”. Son: máquinas de movimiento perpetuo y precognición.

Elementos determinantes en multigrid son los esquemas de relajación y de transferencia entre mallas: restricciones y interpolaciones.

Restricciones:

Interpolaciones:

 

Queremos resolver el sistema Au =f. Multigrid utiliza métodos iterativos para su resolución, por lo que necesita de un aproximación inicial de la solución. Sin embargo, asegurada la convergencia del método, la aproximación inicial puede ser cualquier valor y solo determina la velocidad de convergencia (mayor cuanto mejor sea ésta).

La idea inicial es utilizar una estrategia que utiliza mallas mas gruesas para determinar un valor inicial para la siguiente malla mas fina, de manera que tendriamos algo asi como (denotamos con \Omega^{h} a una malla con un tamaño de intervalo h):

  • Relajación sobre A\boldsymbol{u} =\boldsymbol{ f} en una malla muy gruesa para determinar una primera aproximación para la siguiente malla mas fina.
  • Relajación sobre A\boldsymbol{u} =\boldsymbol{ f} en la malla \Omega^{4h} para obtener un valor inicial para la malla \Omega^{2h}.
  • Relajación sobre A\boldsymbol{u} =\boldsymbol{ f} en la malla \Omega^{2h} para obtener un valor inicial para la malla \Omega^{h}.
  • Relajación sobre A\boldsymbol{u} =\boldsymbol{ f} en la malla \Omega^{h} para obtener la solución final.

Esta es la aproximación por iteraciones anidadas, la mas sencilla. El principal problema es saber cuanto hemos de relajar en cada malla, pues hay que llegar a la última iteración de manera que, al terminar, el error sea del tipo apropiado. Puesto que sabemos que los métodos iterativos sufren en presencia de errores suaves, hay aprovechar esto para saltar de malla.

Una segunda idea incorpora la idea de utilizar la ecuación residual, A\boldsymbol{e} = \boldsymbol{r} = \boldsymbol{f} - A\boldsymbol{v}, para relajar sobre el error:

  • Relajación sobre A\boldsymbol{u} =\boldsymbol{ f} en la malla \Omega^{h} para obtener una aproximación de \boldsymbol{v}^h.
  • Calcular $\boldsymbol{r} = \boldsymbol{f} – A\boldsymbol{v}$.

La

Ya vimos que la discretización de las ecuaciones de conservación en forma Lagrangiana para SPH quedaba:

\frac{d}{dt}\boldsymbol{v}_a = - \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2}) \nabla_a W{ab} = \boldsymbol{F}_a

\frac{d}{dt} e_a = \frac{1}{2} \sum_b m_b (\frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2}) \boldsymbol{v}_{ab}\cdot \nabla_a W{ab} = E_a

con \rho_a = \sum_b m_b W_{ab} y P_a de las EoS.

Para pensar en terminos de ODE solvers, nos fijamos solo en:

\frac{d}{dt} \boldsymbol{v}_a = \boldsymbol{F}_a(t,\boldsymbol{v}_a)

\frac{d}{dt} e_a = E_a(t,e_a)

\frac{d}{dt} \boldsymbol{r}_a = \boldsymbol{v}_a(t,\boldsymbol{r}_a)

Para simplificar, suponiendo como metodo explicito un Euler para el predictor y el método del trapecio como implícito para el corrector, tendriamos:

\boldsymbol{\tilde{v}}_a^{n+1} = \boldsymbol{v}_a^{n} + h \boldsymbol{F}_a(t^n,\boldsymbol{v}_a^{n})

\tilde{e}_a^{n+1} = e_a^{n} + h E_a(t^n,e_a^{n})

\boldsymbol{\tilde{r}}_a^{n+1} = \boldsymbol{r}_a^{n} + h \boldsymbol{v}_a(t^n,\boldsymbol{r}_a^{n})

para la primera y:

\boldsymbol{v}_a^{n+1} = \boldsymbol{v}_a^n + \frac{1}{2}h(\boldsymbol{F}_a(t^n,\boldsymbol{v}_a^n)+\boldsymbol{F}_a(t^{n+1},\boldsymbol{\tilde{v}}_a^{n+1}))

e_a^{n+1} = e_a^n + \frac{1}{2}h(E_a(t^n,e_a^n)+E_a(t^{n+1},\tilde{e}_a^{n+1}))

\boldsymbol{r}_a^{n+1} = \boldsymbol{r}_a^n + \frac{1}{2}h(\boldsymbol{v}_a(t^n,\boldsymbol{r}_a^n)+\boldsymbol{v}_a(t^{n+1},\boldsymbol{\tilde{r}}_a^{n+1}))

para la segunda.

De esta manera, sustituyendo, tenemos:

\boldsymbol{\tilde{v}}_a^{n+1} = \boldsymbol{v}_a^n + h \boldsymbol{F}_a(t^n,\boldsymbol{v}_a^{n}) =

= \boldsymbol{v}_a^n - h \sum_b m_b (\frac{P_a^n}{(\rho_a^n)^2} + \frac{P_b^n}{(\rho_b^n)^2}) \nabla_a W(\boldsymbol{r}_a^n-\boldsymbol{r}_b^n,h)

junio 2017
L M X J V S D
« Feb    
 1234
567891011
12131415161718
19202122232425
2627282930