FireStats error : FireStats: Unknown commit strategy

Aunque siempre podemos hacer cambios de coordenadas, vamos a ver como quedan los esquemas de diferencias finitas en sistemas no rectangulares: coordenadas cilíndricas, $latex (rho,phi, z)$, y coordenadas esféricas, $latex (r,theta,phi)$. Nos centraremos en la ecuación de Poisson aunque la técnica se puede extender de manera inmediata a cualquier tipo de PDE.

En coordenadas cilíndricas podemos escribir:

$latex nabla cdot nabla u = frac{partial^2}{partial rho^2}u + frac{1}{rho}frac{partial}{partial rho}u + frac{1}{rho^2}frac{partial^2}{partial phi^2}u + frac{partial^2}{partial z^2} u = f$,

que podemos discretizar como:

$latex frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{(Delta rho)^2} + $

$latex + frac{1}{rho_{i,j,k}}frac{u_{i+1,j,k}-u_{i-1,j,k}}{2Delta rho} + $

$latex + frac{1}{rho_{i,j,k}^2} frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{(Delta phi)^2} + $

$latex + frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{(Delta z)^2} = f_{i,j,k}$

En coordenadas esféricas tenemos:

$latex nabla cdot nabla u = frac{partial^2}{partial r^2}u + frac{2}{r} frac{partial}{partial r}u + frac{1}{r^2}frac{partial^2}{partial theta^2}u + frac{1}{r^2sintheta} frac{partial}{partial theta}u + frac{1}{r^2 sin^2theta} frac{partial^2}{partial phi^2}u = f$

que podemos discretizar como:

$latex frac{u_{i-1,j,k-2u_{i,j,k}+u_{i+1,j,k}}}{(Delta r)^2} + $

$latex + frac{2}{r_{i,j,k}} frac{u_{i+1,j,k}+u_{i-1,j,k}}{2Delta r} + $

$latex + frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{(r_{i,j,k} Delta theta)^2} + $

$latex + frac{1}{r_{i,j,k}^2 sin phi_{i,j,k}} frac{u_{i,j+1,k}-u_{i,j-1,k}}{2 Delta phi} + $

$latex + frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{(r_{i,j,k} sin phi_{i,j,k} Delta phi)^2} = f_{i,j,k}$

Tags: , , , , ,

En $latex n$ dimensiones, el operador Laplaciano queda como:

$latex Delta u= sum_{i=1}^n frac{partial^2}{partial x_i^2}u$

en coordenadas cartesianas, y como:

$latex Delta u = frac{partial}{partial r^2}u + frac{n-1}{r}frac{partial}{partial r}u + frac{1}{r^2}Delta_{S^{n-1}}u$

en esféricas, donde $latex Delta_{S^{n-1}}$ es el operador de Laplace-Beltrami, una generalización del Laplaciano para funciones definidas sobre variedades,  en la $latex (n-1)$-esfera ($latex S^{n-1}$), el operador Laplaciano esférico.

Un punto es un tensor sin índices, un vector es un tensor con $latex 1$ índice, una matriz es un tensor con $latex 2$ índices, etc. Cuando discreticemos una PDE en $latex n$ dimensiones, llegaremos a un tensor con $latex n$ índices y $latex 2n$ tensores con $latex n-1$ índices para las condiciones en las fronteras.

Tags: , , , , ,

Vamos a suponer $latex n=3$ para reducir el tamaño de las matrices.

Empezamos suponiendo que conocemos:

$latex frac{partial}{partial x}|_{0,0,}u, frac{partial}{partial x}|_{0,1}u, frac{partial}{partial x}|_{0,2}u$

$latex frac{partial}{partial y}|_{0,0}u, frac{partial}{partial y}|_{1,0}u$

$latex frac{partial}{partial y}|_{0,2}u, frac{partial}{partial y}|_{1,2}u$

$latex u|_{2,0}, u|_{2,1}, u|_{2,2}$

Discretizamos:

$latex frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{h^2} + frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{h^2} = f_{0,0}$

$latex frac{u_{-1,1}-2u_{0,1}+u_{1,1}}{h^2} + frac{u_{0,0}-2u_{0,1}+u_{0,2}}{h^2} = f_{0,1}$

$latex frac{u_{-1,2}-2u_{0,2}+u_{1,2}}{h^2} + frac{u_{0,1}-2u_{0,2}+u_{0,3}}{h^2} = f_{0,2}$

$latex frac{u_{0,0}-2u_{1,0}+u_{2,0}}{h^2} + frac{u_{1,-1}-2u_{1,0}+u_{1,1}}{h^2} = f_{1,0}$

$latex frac{u_{0,1}-2u_{1,1}+u_{2,1}}{h^2} + frac{u_{1,0}-2u_{1,1}+u_{1,2}}{h^2} = f_{1,1}$

$latex frac{u_{0,2}-2u_{1,2}+u_{2,2}}{h^2} + frac{u_{1,1}-2u_{1,2}+u_{1,3}}{h^2} = f_{1,2}$

En las fronteras, sabemos que:

$latex frac{u_{1,0}-u_{-1,0}}{2h} = frac{partial}{partial x}|_{0,0}u Leftrightarrow u_{-1,0}=u_{1,0}-2h frac{partial}{partial x}|_{0,0}u$

$latex frac{u_{1,1}-u_{-1,1}}{2h} = frac{partial}{partial x}|_{0,1}u Leftrightarrow u_{-1,1}=u_{1,1}-2h frac{partial}{partial x}|_{0,1}u$

$latex frac{u_{1,2}-u_{-1,2}}{2h} = frac{partial}{partial x}|_{0,2}u Leftrightarrow u_{-1,2}=u_{1,2}-2h frac{partial}{partial x}|_{0,2}u$

$latex frac{u_{0,1}-u_{0,-1}}{2h} = frac{partial}{partial y}|_{0,0}u Leftrightarrow u_{0,-1}=u_{0,1}-2h frac{partial}{partial y}|_{0,0}u$

$latex frac{u_{1,1}-u_{1,-1}}{2h} = frac{partial}{partial y}|_{1,0}u Leftrightarrow u_{1,-1}=u_{1,1}-2h frac{partial}{partial y}|_{1,0}u$

$latex frac{u_{0,3}-u_{0,1}}{2h} = frac{partial}{partial y}|_{0,2}u Leftrightarrow u_{0,3}=u_{0,1}+2h frac{partial}{partial y}|_{0,2}u$

$latex frac{u_{1,3}-u_{1,1}}{2h} = frac{partial}{partial y}|_{1,2}u Leftrightarrow u_{1,3}=u_{1,1}+2h frac{partial}{partial y}|_{1,2}u$

La matriz queda:

$latex left(
begin{array}{ccc|ccc}
-4 & 2 & 0 & 2 & 0 & 0 \
1 & -4 & 1 & 0 & 2 & 0 \
0 & 2 & -4 & 0 & 0 & 2 \ hline
1 & 0 & 0 & -4 & 2 & 0 \
0 & 1 & 0 & 1 & -4 & 1 \
0 & 0 & 1 & 0 & 2 & -4
end{array}
right)$

Simetrizable como:

$latex left(
begin{array}{ccc|ccc}
-1 & frac{1}{2} & 0 & frac{1}{2} & 0 & 0 \
frac{1}{2} & -2 & frac{1}{2} & 0 & 1 & 0 \
0 & frac{1}{2} & -1 & 0 & 0 & frac{1}{2} \ hline
frac{1}{2} & 0 & 0 & -2 & 1 & 0 \
0 & 1 & 0 & 1 & -4 & 1 \
0 & 0 & frac{1}{2} & 0 & 1 & -2
end{array}
right)$

Tenemos $latex 6$ ecuaciones con $latex 6$ incognitas y la matriz tiene rango $latex 6$, por lo que la solución es única.

En el segundo caso, suponemos que todas las fronteras son Neumann:

$latex frac{partial}{partial x}|_{0,0}u, frac{partial}{partial x}|_{0,1}u, frac{partial}{partial x}|_{0,2}u$

$latex frac{partial}{partial y}|_{0,0}u, frac{partial}{partial y}|_{1,0}u, frac{partial}{partial y}|_{2,0}u$

$latex frac{partial}{partial y}|_{0,2}u, frac{partial}{partial y}|_{1,2}u, frac{partial}{partial y}|_{2,2}u$

$latex frac{partial}{partial x}|_{2,0}u, frac{partial}{partial x}|_{2,1}u, frac{partial}{partial x}|_{2,2}u$

Si discretizamos:

$latex frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{h^2} + frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{h^2} = f_{0,0}$

$latex frac{u_{-1,1}-2u_{0,1}+u_{1,1}}{h^2} + frac{u_{0,0}-2u_{0,1}+u_{0,2}}{h^2} = f_{0,1}$

$latex frac{u_{-1,2}-2u_{0,2}+u_{1,2}}{h^2} + frac{u_{0,1}-2u_{0,2}+u_{0,3}}{h^2} = f_{0,2}$

$latex frac{u_{0,0}-2u_{1,0}+u_{2,0}}{h^2} + frac{u_{1,-1}-2u_{1,0}+u_{1,1}}{h^2} = f_{1,0}$

$latex frac{u_{0,1}-2u_{1,1}+u_{2,1}}{h^2} + frac{u_{1,0}-2u_{1,1}+u_{1,2}}{h^2} = f_{1,1}$

$latex frac{u_{0,2}-2u_{1,2}+u_{2,2}}{h^2} + frac{u_{1,1}-2u_{1,2}+u_{1,3}}{h^2} = f_{1,2}$

$latex frac{u_{1,0}-2u_{2,0}+u_{3,0}}{h^2} + frac{u_{2,-1}-2u_{2,0}+u_{2,1}}{h^2} = f_{2,0}$

$latex frac{u_{1,1}-2u_{2,1}+u_{3,1}}{h^2} + frac{u_{2,0}-2u_{2,1}+u_{2,2}}{h^2} = f_{2,1}$

$latex frac{u_{1,2}-2u_{2,2}+u_{3,2}}{h^2} + frac{u_{2,1}-2u_{2,2}+u_{2,3}}{h^2} = f_{2,2}$

En las fronteras, sabemos que:

$latex frac{u_{1,0}-u_{-1,0}}{2h} = frac{partial}{partial x}|_{0,0}u Leftrightarrow u_{-1,0}=u_{1,0}-2h frac{partial}{partial x}|_{0,0}u$

$latex frac{u_{1,1}-u_{-1,1}}{2h} = frac{partial}{partial x}|_{0,1}u Leftrightarrow u_{-1,1}=u_{1,1}-2h frac{partial}{partial x}|_{0,1}u$

$latex frac{u_{1,2}-u_{-1,2}}{2h} = frac{partial}{partial x}|_{0,2}u Leftrightarrow u_{-1,2}=u_{1,2}-2h frac{partial}{partial x}|_{0,2}u$

$latex frac{u_{0,1}-u_{0,-1}}{2h} = frac{partial}{partial y}|_{0,0}u Leftrightarrow u_{0,-1}=u_{0,1}-2h frac{partial}{partial y}|_{0,0}u$

$latex frac{u_{1,1}-u_{1,-1}}{2h} = frac{partial}{partial y}|_{1,0}u Leftrightarrow u_{1,-1}=u_{1,1}-2h frac{partial}{partial y}|_{1,0}u$

$latex frac{u_{2,1}-u_{2,-1}}{2h} = frac{partial}{partial y}|_{2,0}u Leftrightarrow u_{2,-1}=u_{2,1}-2h frac{partial}{partial y}|_{2,0}u$

$latex frac{u_{0,3}-u_{0,1}}{2h} = frac{partial}{partial y}|_{0,2}u Leftrightarrow u_{0,3}=u_{0,1}+2h frac{partial}{partial y}|_{0,2}u$

$latex frac{u_{1,3}-u_{1,1}}{2h} = frac{partial}{partial y}|_{1,2}u Leftrightarrow u_{1,3}=u_{1,1}+2h frac{partial}{partial y}|_{1,2}u$

$latex frac{u_{2,3}-u_{2,1}}{2h} = frac{partial}{partial y}|_{2,2}u Leftrightarrow u_{2,3}=u_{2,1}+2h frac{partial}{partial y}|_{2,2}u$

$latex frac{u_{3,0}-u_{1,0}}{2h} = frac{partial}{partial x}|_{2,0}u Leftrightarrow u_{3,0}=u_{1,0}+2h frac{partial}{partial x}|_{2,0}u$

$latex frac{u_{3,1}-u_{1,1}}{2h} = frac{partial}{partial x}|_{2,1}u Leftrightarrow u_{3,1}=u_{1,1}+2h frac{partial}{partial x}|_{2,1}u$

$latex frac{u_{3,2}-u_{1,2}}{2h} = frac{partial}{partial x}|_{2,2}u Leftrightarrow u_{3,2}=u_{1,2}+2h frac{partial}{partial x}|_{2,2}u$

La matriz, por tanto, queda:

$latex text{A6}=left(
begin{array}{ccc|ccc|ccc}
-4 & 2 & 0 & 2 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 \
0 & 2 & -4 & 0 & 0 & 2 & 0 & 0 & 0 \ hline
1 & 0 & 0 & -4 & 2 & 0 & 1 & 0 & 0 \
0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \
0 & 0 & 1 & 0 & 2 & -4 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 2 & 0 & 0 & -4 & 2 & 0 \
0 & 0 & 0 & 0 & 2 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 2 & 0 & 2 & -4
end{array}
right)$

Simetrizable como:

$latex text{A6s}=left(
begin{array}{ccc|ccc|ccc}
-1 & 1/2 & 0 & 1/2 & 0 & 0 & 0 & 0 & 0 \
1/2 & -2 & 1/2 & 0 & 1 & 0 & 0 & 0 & 0 \
0 & 1/2 & -1 & 0 & 0 & 1/2 & 0 & 0 & 0 \ hline
1/2 & 0 & 0 & -2 & 1 & 0 & 1/2 & 0 & 0 \
0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \
0 & 0 & 1/2 & 0 & 1 & -2 & 0 & 0 & 1/2 \ hline
0 & 0 & 0 & 1/2 & 0 & 0 & -1 & 1/2 & 0 \
0 & 0 & 0 & 0 & 1 & 0 & 1/2 & -2 & 1/2 \
0 & 0 & 0 & 0 & 0 & 1/2 & 0 & 1/2 & -1
end{array}
right)$

En este caso, tenemos $latex 9$ ecuaciones con $latex 9$ incognitas pero la matriz tiene rango $latex 8$, por lo que tenemos infinitas soluciones. Hay que conservar.

Tags: , , , , ,

Suponemos $latex Delta u = f$ en $latex 2D$, es decir,

$latex frac{partial^2}{partial x^2}u(x,y) + frac{partial^2}{partial y^2}u(x,y) = f(x,y)$.

Miraremos como queda la matriz del sistema al discretizar, como simetrizarla y su rango en tres casos: condición Neuman respecto $latex x$ en una frontera, condición Neumann respecto $latex y$ en una frontera y condición Neumann respecto $latex x$ e $latex y$ en dos fronteras.

Discretizamos con $latex n=5$. Si todas las condiciones fueran Dirichlet, la matriz quedaría:

$latex A_1 = left(
begin{array}{ccc|ccc|ccc}
-4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \
0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \ hline
1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \
0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \
0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4
end{array}
right) $.

En este caso, $latex A_1 in mathcal{M}(9 times 9)$ y simétrica, lo que permite tratar de manera conjunta los problemas de existencia y unicidad de solución. Si calculamos su rango obtenemos $latex 9$ por lo que existe solución y es única. Desde el punto de vista algebraico, es el punto $latex (u_{1,1},u_{1,2},u_{1,3},u_{2,1},u_{2,2},u_{2,3},u_{3,1},u_{3,2},u_{3,3})$ intersección de $latex 9$ hiperplanos

$latex -4x_{1,1} + x_{1,2} + x_{2,1} = f_{1,1}$,

$latex x_{1,1}-4x_{1,2}+x_{1,3} + x_{2,2} = f_{1,2}$,

$latex ldots$

en el espacio $latex mathbb{R}^9$.

Si condiremos conocidos $latex frac{partial}{partial x}|_{0,1}u, frac{partial}{partial x}|_{0,2}u, frac{partial}{partial x}|_{0,3}u$ en lugar de $latex u_{0,1}, u_{0,2}, u_{0,3}$ ($latex u_{0,0}$ y $latex u_{0,4}$ son conocidos por las otras fronteras que son Dirichelt), tenemos:

$latex A_2 = left(
begin{array}{ccc|ccc|ccc|ccc}
-4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & -4 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \ hline
1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \ hline
0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4
end{array}
right)$

de manera que $latex A_2 in mathcal{M}(12 times 12)$ y no es simétrica. Sin embargo es facilmente simetrizable dividiendo las tres primera filas (hacemos lo mismo en el termino independiente) por $latex 2$:

$latex A_2 = left(
begin{array}{ccc|ccc|ccc|ccc}
-2 & frac{1}{2} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
frac{1}{2} & -2 & frac{1}{2} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & frac{1}{2} & -2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \ hline
1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \ hline
0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4
end{array}
right)$

Tenemos $latex 12$ incognitas ($latex u_{i,j}$ con $latex i=0..3$ y $latex j=1..3$) y el rango de $latex A_2$ es $latex 12$, por lo que la solución, nuevamente, es única.

Para el caso en el que conocemos $latex frac{partial}{partial y}|_{1,0}u, frac{partial}{partial y}|_{2,0}u, frac{partial}{partial y}|_{3,0}u$ en lugar de $latex u_{1,0}, u_{2,0}, u_{3,0}$, si el orden que tomamos es el contrario al tomado anteriormente llegaremos a la misma estructura de antes. Sin embargo, como en el siguiente caso nos veremos obligados a seleccionar uno de los dos, vamos a ver como queda este caso utilizando el mismo orden que antes:

$latex A_3 = left(
begin{array}{cccc|cccc|cccc}
-4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ hline
1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 \
0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \
0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \
0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4
end{array}
right)$

que podemos simetrizar facilmente y queda:

$latex A_3 = left(
begin{array}{cccc|cccc|cccc}
-2 & 1 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ hline
frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 \
0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \
0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \
0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4
end{array}
right)$

Tenemos $latex 12$ ecuaciones con $latex 12$ incognitas ($latex u_{i,j}$ con $latex i=1..3$ y $latex j=0..3$) y el rango de $latex A_3$ es $latex 12$, por lo que la solución es única.

Finalmente, suponemos conocidos $latex frac{partial}{partial x}|_{0,0}u, frac{partial}{partial x}|_{0,1}u, frac{partial}{partial x}|_{0,2}u, frac{partial}{partial x}|_{0,3}u, frac{partial}{partial y}|_{0,0}u, frac{partial}{partial y}|_{1,0}u, frac{partial}{partial y}|_{2,0}u, frac{partial}{partial y}|_{3,0}u$ que incorpora $latex 7$ ecuaciones mas a las $latex 9$ que ya teniamos por lo que nos queda una matrix $latex A_4 in mathcal{M}(16 times 16)$:

$latex left(
begin{array}{cccc|cccc|cccc|cccc}
-4 & 2 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
1 & -4 & 1 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & -4 & 1 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & -4 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ hline
1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ hline
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4
end{array}
right)$,

simetrizable dividiendo la fila correspondiente a $latex u_{0,0}$ por $latex 4$, y las correspondientes a $latex u_{0,1}, u_{0,2}, u_{0,3}, u_{1,0},u_{2,0}, u_{3,0}$  por $latex 2$, quedando:

$latex left(
begin{array}{cccc|cccc|cccc|cccc}
-1 & frac{1}{2} & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
frac{1}{2} & -2 & frac{1}{2} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & frac{1}{2} & -2 & frac{1}{2} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & frac{1}{2} & -2 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ hline
frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \ hline
0 & 0 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \ hline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4
end{array}
right)$,

con lo que el sistema vuelve a ser compatible y determinado.

Tags: , , , , , ,

Suponemos $latex n=5$. En el caso de tener todas las fronteras con condiciones Dirichlet:

$latex frac{u_{0,1} -2u_{1,1} + u_{2,1}}{h^2} + frac{u_{1,0} -2u_{1,1} + u_{1,2}}{h^2} = f_{1,1}$ para $latex i,j=1,1$,

$latex frac{u_{0,2} -2u_{1,2} + u_{2,2}}{h^2} + frac{u_{1,1} -2u_{1,2} + u_{1,3}}{h^2} = f_{1,2}$ para $latex i,j=1,2$,

$latex frac{u_{0,3} -2u_{1,3} + u_{2,3}}{h^2} + frac{u_{1,2} -2u_{1,3} + u_{1,4}}{h^2} = f_{1,3}$ para $latex i,j=1,3$,

$latex frac{u_{1,1} -2u_{2,1} + u_{3,1}}{h^2} + frac{u_{2,0} -2u_{2,1} + u_{2,2}}{h^2} = f_{2,1}$ para $latex i,j=2,1$,

$latex frac{u_{1,2} -2u_{2,2} + u_{3,2}}{h^2} + frac{u_{2,1} -2u_{2,2} + u_{2,3}}{h^2} = f_{2,2}$ para $latex i,j=2,2$,

$latex frac{u_{1,3} -2u_{2,3} + u_{3,3}}{h^2} + frac{u_{2,2} -2u_{2,3} + u_{2,4}}{h^2} = f_{2,3}$ para $latex i,j=2,3$,

$latex frac{u_{2,1} -2u_{3,1} + u_{4,1}}{h^2} + frac{u_{3,0} -2u_{3,1} + u_{3,2}}{h^2} = f_{3,1}$ para $latex i,j=3,1$,

$latex frac{u_{2,2} -2u_{3,2} + u_{4,2}}{h^2} + frac{u_{3,1} -2u_{3,2} + u_{3,3}}{h^2} = f_{3,2}$ para $latex i,j=3,2$,

$latex frac{u_{2,3} -2u_{3,3} + u_{4,3}}{h^2} + frac{u_{3,2} -2u_{3,3} + u_{3,4}}{h^2} = f_{3,3}$ para $latex i,j=3,3$,

de donde:

$latex begin{bmatrix} f_{1,1} -frac{u_{1,0} + u_{0,1}}{h^2} & f_{1,2} – frac{u_{0,2}}{h^2} & f_{1,3} – frac{u_{0,3}+u_{1,4}}{h^2} \ f_{2,1} -frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} – frac{u_{2,4}}{h^2} \ f_{3,1} – frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} – frac{u_{4,2}}{h^2} & f_{3,3} – frac{u_{4,3}+u_{3,4}}{h^2} end{bmatrix}$

En forma de matriz por bloques (para pensar en la simetrización):

$latex frac{1}{h^2} begin{bmatrix} -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \ 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \ 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \ 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \ 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \ 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \ 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \ 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 end{bmatrix} u_{i,j} = begin{bmatrix} f_{1,1} -frac{u_{1,0} + u_{0,1}}{h^2} \ f_{1,2} – frac{u_{0,2}}{h^2} \ f_{1,3} – frac{u_{0,3}+u_{1,4}}{h^2} \ f_{2,1} -frac{u_{2,0}}{h^2} \ f_{2,2} \ f_{2,3} – frac{u_{2,4}}{h^2} \ f_{3,1} – frac{u_{3,0}+u_{4,1}}{h^2} \ f_{3,2} – frac{u_{4,2}}{h^2} \ f_{3,3} – frac{u_{4,3}+u_{3,4}}{h^2} end{bmatrix}$

¿Qué pasa ahora si en lugar de conocer $latex u_{0,1}, u_{0,2}, u_{0,3}$ conocemos $latex frac{partial}{partial x}|_{0,1}u, frac{partial}{partial x}|_{0,2}u, frac{partial}{partial x}u|_{0,3}$? Necesitamos tres ecuaciones mas:

$latex frac{u_{-1,1} -2u_{0,1} + u_{1,1}}{h^2} + frac{u_{0,0} -2u_{0,1} + u_{0,2}}{h^2} = f_{0,1}$ para $latex i,j=0,1$

$latex frac{u_{-1,2} -2u_{0,2} + u_{1,2}}{h^2} + frac{u_{0,1} -2u_{0,2} + u_{0,3}}{h^2} = f_{0,2}$ para $latex i,j=0,2$

$latex frac{u_{-1,3} -2u_{0,3} + u_{1,3}}{h^2} + frac{u_{0,2} -2u_{0,3} + u_{0,4}}{h^2} = f_{0,3}$ para $latex i,j=0,3$

y

$latex frac{u_{1,1}-u_{-1,1}}{2h} = frac{partial}{partial x}|_{0,1}u Leftrightarrow u_{-1,1} = u_{1,1} – 2h , frac{partial}{partial x}|_{0,1}u$

$latex frac{u_{1,2}-u_{-1,2}}{2h} = frac{partial}{partial x}|_{0,2}u Leftrightarrow u_{-1,2} = u_{1,2} – 2h , frac{partial}{partial x}|_{0,2}u$

$latex frac{u_{1,3}-u_{-1,3}}{2h} = frac{partial}{partial x}|_{0,3}u Leftrightarrow u_{-1,3} = u_{1,3} – 2h , frac{partial}{partial x}|_{0,3}u$

por lo que:

$latex begin{bmatrix} f_{0,1} +frac{2h , frac{partial}{partial x}|_{0,1}u – u_{0,0}}{h^2} & f_{0,2} + frac{2h , frac{partial}{partial x}|_{0,2}u}{h^2} & f_{0,3} + frac{ 2h , frac{partial}{partial x}|_{0,3}u – u_{0,4}}{h^2} \ f_{1,1} -frac{u_{1,0}}{h^2} & f_{1,2} & f_{1,3} – frac{u_{1,4}}{h^2} \ f_{2,1} -frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} – frac{u_{2,4}}{h^2} \ f_{3,1} – frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} – frac{u_{4,2}}{h^2} & f_{3,3} – frac{u_{4,3}+u_{3,4}}{h^2} end{bmatrix}$

La matriz queda:

$latex frac{1}{h^2} begin{bmatrix} -4 & 1 & 0 & 2 & 0 & 0 & ldots \ 1 & -4 & 1 & 0 & 2 & 0 & ldots \ 0 & 1 & -4 & 0 & 0 & 2 & ldots \ 1 & 0 & 0 & -4 & 1 & 0 & ldots \ 0 & 1 & 0 & 1 & -4 & 1 & ldots \ 0 & 0 & 1 & 0 & 1 & -4 & ldots \ vdots & vdots & vdots & vdots & vdots & vdots & ddots end{bmatrix} $

Que podemos simetrizar:

$latex frac{1}{h^2} begin{bmatrix} -2 & frac{1}{2} & 0 & 1 & 0 & 0 & ldots \ frac{1}{2} & -2 & frac{1}{2} & 0 & 1 & 0 & ldots \ 0 & frac{1}{2} & -2 & 0 & 0 & 1 & ldots \ 1 & 0 & 0 & -4 & 1 & 0 & ldots \ 0 & 1 & 0 & 1 & -4 & 1 & ldots \ 0 & 0 & 1 & 0 & 1 & -4 & ldots \ vdots & vdots & vdots & vdots & vdots & vdots & ddots end{bmatrix} $

con:

$latex begin{bmatrix} frac{1}{2}(f_{0,1} +frac{2h , frac{partial}{partial x}|_{0,1}u – u_{0,0}}{h^2}) & frac{1}{2}(f_{0,2} + frac{2h , frac{partial}{partial x}|_{0,2}u}{h^2}) & frac{1}{2}(f_{0,3} + frac{ 2h , frac{partial}{partial x}|_{0,3}u – u_{0,4}}{h^2}) \ f_{1,1} -frac{u_{1,0}}{h^2} & f_{1,2} & f_{1,3} – frac{u_{1,4}}{h^2} \ f_{2,1} -frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} – frac{u_{2,4}}{h^2} \ f_{3,1} – frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} – frac{u_{4,2}}{h^2} & f_{3,3} – frac{u_{4,3}+u_{3,4}}{h^2} end{bmatrix}$

Si las condiciones las tenemos sobre la derivada en el extremo opuesto llegaremos a la misma estructura pero en la parte inferior de la frontera y de la matriz.

Si las condiciones las tenemos sobre derivadas en la otra dirección, podemos llegar también a estas estructuras tomando el orden de variables donde tiene prioridad la variable contraria a la tomada en los casos anteriores.

Tags: , , , , ,

En el post anterior hablamos sobre condiciones de frontera y su transferencia entre mallas pero no comentamos en el caso de que las condición haga referencia al valor de la derivada y no al de la función: condición de Neumann.

En $latex 1D$ supongamos que ahora tenemos $latex frac{partial^2}{partial x^2}u = f$ en $latex [a,b]$ con $latex u(a)=u_a$ pero $latex frac{partial}{partial x} = du_b$. Suponiendo de nuevo $latex n=8$, las ecuaciones nos quedan:

$latex frac{u_0 -2u_1 + u_2}{h^2} = f_1$ para $latex i=1$,

$latex frac{u_1 -2u_2 + u_3}{h^2} = f_2$ para $latex i=2$,

$latex frac{u_2 -2u_3 + u_4}{h^2} = f_3$ para $latex i=3$,

$latex frac{u_3 -2u_4 + u_5}{h^2} = f_4$ para $latex i=4$,

$latex frac{u_4 -2u_5 + u_6}{h^2} = f_5$ para $latex i=5$,

$latex frac{u_5 -2u_6 + u_7}{h^2} = f_6$ para $latex i=6$,

$latex frac{u_6 -2u_7 + u_8}{h^2} = f_7$ para $latex i=7$,

La única diferencia con respecto al caso anterior es que, en la primera ecuación, desconocemos el valor de $latex u_0$ pero  conocemos el de su primera derivada. Sabemos que:

$latex frac{u_1 – u_{-1}}{2h} = frac{d}{dx}u_{0} = du_0$,

que, despejando, nos da:

$latex u_1 – u_{-1} = 2h , du_0 Leftrightarrow u_{-1} = u_1 – 2h , du_0$,

Como tenemos una incognita mas por determinar, añadimos una nueva ecuación:

$latex frac{u_{-1} -2u_0 + u_1}{h^2} = f_0$ para $latex i=0$,

donde reescribimos el valor de $latex u_{-1}$ según acabamos de determinar:

$latex frac{ u_1 – 2h , du_0 -2u_0 + u_1}{h^2} = frac{ -2u_0 + 2u_1 – 2h , du_0}{h^2} = f_0$.

Por lo tanto,  en forma matricial tenemos:

$latex frac{1}{h^2} begin{bmatrix} -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \ 1 & -2 & 1 & 0 & 0 & 0 & 0 &0 \ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\ 0 & 0 & 1 & -2 & 1 & 0 & 0 &0 \ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 end{bmatrix} begin{bmatrix} u_0 \ u_1 \ u_2 \ u_3 \ u_4 \ u_5 \ u_6 \u_7 end{bmatrix} = begin{bmatrix} frac{1}{2} (f_0 + frac{2h , du_0}{h^2}) \ f_1 \ f_2 \ f_3 \f_4 \ f_5 \ f_6 \ f_7 – frac{u_8}{h^2}end{bmatrix}$.

De la misma manera, en el caso por el otro extremo, llegariamos a:

$latex frac{1}{h^2} begin{bmatrix} -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \ 1 & -2 & 1 & 0 & 0 & 0 & 0 &0 \ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\ 0 & 0 & 1 & -2 & 1 & 0 & 0 &0 \ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 end{bmatrix} begin{bmatrix} u_1 \ u_2 \ u_3 \ u_4 \ u_5 \ u_6 \ u_7 \u_8 end{bmatrix} = begin{bmatrix} f_1 – frac{u_0}{h^2} \ f_2 \ f_3 \ f_4 \f_5 \ f_6 \ f_7 \ frac{1}{2}(f_8 + frac{2h , du_8}{h^2})end{bmatrix}$.

En resumen, básicamente hay que hacer dos trabajos: en primer lugar, construir el termino independiente de manera apropiada para incorporar la información de las fronteras; en segundo, llegados a los extremos, escoger entre $latex -2$ y $latex -1$ en la diagonal en función de si es Dirichlet o Neumann.

Tags: , , , , ,

A ver si nos aclaramos sobre como van las condiciones frontera en las transiciones entre mallas…

Empezamos en $latex 1D$. Tenemos $latex Delta u = f$, o lo que es lo mismo, $latex u_{xx} = f$ (en este caso es una ODE pero bueno…) definida en $latex [a,b]$ con $latex u(a) = u_a$ y $latex u(b) = u_b$. Vamos a suponer una discretización en una malla con $latex n+1$ nodos con $latex u_i$ con $latex i=1..(n-1)$ puntos interiores y $latex u_0 = u_a$ y $latex u_n = u_b$. En la discretización tenemos:

$latex frac{u_{i-1} -2u_i + u_{i+1}}{h^2} = f_i$ donde $latex h = frac{1}{n}$.

Si escribimos todas las ecuaciones para todos los puntos interiores (si $latex n=8$ entonces $latex i=1..7$) tenemos :

$latex frac{u_0 -2u_1 + u_2}{h^2} = f_1$ para $latex i=1$,

$latex frac{u_1 -2u_2 + u_3}{h^2} = f_2$ para $latex i=2$,

$latex frac{u_2 -2u_3 + u_4}{h^2} = f_3$ para $latex i=3$,

$latex frac{u_3 -2u_4 + u_5}{h^2} = f_4$ para $latex i=4$,

$latex frac{u_4 -2u_5 + u_6}{h^2} = f_5$ para $latex i=5$,

$latex frac{u_5 -2u_6 + u_7}{h^2} = f_6$ para $latex i=6$,

$latex frac{u_6 -2u_7 + u_8}{h^2} = f_7$ para $latex i=7$,

que en forma matricial y despejando $latex u_0$ y $latex u_8$ que son conocidos queda:

$latex frac{1}{h^2} begin{bmatrix} -2 & 1 & 0 & 0 & 0 & 0 & 0 \ 1 & -2 & 1 & 0 & 0 & 0 & 0 \ 0 & 1 & -2 & 1 & 0 & 0 & 0 \ 0 & 0 & 1 & -2 & 1 & 0 & 0 \ 0 & 0 & 0 & 1 & -2 & 1 & 0 \ 0 & 0 & 0 & 0 & 1 & -2 & 1 \ 0 & 0 & 0 & 0 & 0 & 1 & -2 end{bmatrix} begin{bmatrix} u_1 \ u_2 \ u_3 \ u_4 \ u_5 \ u_6 \u_7 end{bmatrix} = begin{bmatrix} f_1 – frac{u_0}{h^2} \ f_2 \ f_3 \f_4 \ f_5 \ f_6 \ f_7 – frac{u_8}{h^2}end{bmatrix}$.

Cuando pasamos a una malla de $latex n=4$ tenemos que la matriz es $latex 3 times 3$ y tiene la misma estructura pero con $latex frac{1}{(2h)^2}$. En este caso, si restringimos la $latex vec{f}$  nos queda:

$latex frac{1}{4} begin{bmatrix} 1 & 2 & 1 & 0 & 0 & 0 & 0 \ 0 & 0 & 1 & 2 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & 2 & 1 end{bmatrix} begin{bmatrix} f_1 – frac{u_0}{h^2} \ f_2 \ f_3 \ f_4 \ f_5 \ f_6 \ f_7 – frac{u_8}{h^2} end{bmatrix} = begin{bmatrix} frac{f_1 – frac{u_0}{h^2} +2 f_2 + f_3}{4} \ frac{f_3+2 f_4 + f_5}{4} \ frac{f_5 + 2 f_6 + f7 – frac{u_8}{h^2}}{4}end{bmatrix}$.

¿Qué pasa en este caso si restringimos por separado las fuentes y los valores en la frontera? Por un lado tenemos:

$latex frac{1}{4} begin{bmatrix} 1 & 2 & 1 & 0 & 0 & 0 & 0 \ 0 & 0 & 1 & 2 & 1 & 0 & 0 \ 0 & 0 & 0 & 0 & 1 & 2 & 1 end{bmatrix} begin{bmatrix} f_1 \ f_2 \ f_3 \ f_4 \ f_5 \ f_6 \ f_7 end{bmatrix} = begin{bmatrix} frac{f_1 +2 f_2 + f_3}{4} \ frac{f_3+2 f_4 + f_5}{4} \ frac{f_5 + 2 f_6 + f7 }{4}end{bmatrix}$,

y por otro, como $latex u_0$ y $latex u_8$ no cambian, al despejar nos quedan $latex -frac{u_0}{(2h)^2}$ y $latex -frac{u_8}{(2h)^2}$, por lo que tenemos:

$latex begin{bmatrix} frac{f_1 +2 f_2 + f_3}{4} – frac{u_0}{(2h)^2} \ frac{f_3+2 f_4 + f_5}{4} \ frac{f_5 + 2 f_6 + f7 }{4} – frac{u_8}{(2h)^2} end{bmatrix}$,

que es equivalente a lo encontrado anteriormente.

¿Que pasará en $latex 2D$? Vamos a verlo. Tenemos ahora:

$latex Delta u = f$

como:

$latex frac{partial^2}{partial x^2} u(x,y) + frac{partial^2}{partial y^2} u(x,y) = f(x.y)$.

Suponemos $latex n=8$. Por un lado tenemos que las fuentes menos las fronteras nos da:

$latex begin{bmatrix} f_{1,1}-frac{u_{1,0}+u_{0,1}}{h^2} & f_{1,2}-frac{u_{0,2}}{h^2} & f_{1,3}-frac{u_{0,3}}{h^2} & f_{1,4}-frac{u_{0,4}}{h^2} & f_{1,5}-frac{u_{0,5}}{h^2} & f_{1,6}-frac{u_{0,6}}{h^2} & f_{1,7}-frac{u_{1,8}+u_{0,7}}{h^2} \ f_{2,1}-frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} & f_{2,4} & f_{2,5} & f_{2,6} & f_{2,7}-frac{u_{2,8}}{h^2} \ f_{3,1}-frac{u_{3,0}}{h^2} & f_{3,2} & f_{3,3} & f_{3,4} & f_{3,5} & f_{3,6} & f_{3,7}-frac{u_{3,8}}{h^2} \ f_{4,1}-frac{u_{4,0}}{h^2} & f_{4,2} & f_{4,3} & f_{4,4} & f_{4,5} & f_{4,6} & f_{4,7}-frac{u_{4,8}}{h^2} \ f_{5,1}-frac{u_{5,0}}{h^2} & f_{5,2} & f_{5,3} & f_{5,4} & f_{5,5} & f_{5,6} & f_{5,7}-frac{u_{5,8}}{h^2} \ f_{6,1}-frac{u_{6,0}}{h^2} & f_{6,2} & f_{6,3} & f_{6,4} & f_{6,5} & f_{6,6} & f_{6,7}-frac{u_{6,8}}{h^2} \ f_{7,1}-frac{u_{7,0}+u_{8,1}}{h^2} & f_{7,2}-frac{u_{8,2}}{h^2} & f_{7,3}-frac{u_{8,3}}{h^2} & f_{7,4}-frac{u_{8,4}}{h^2} & f_{7,5}-frac{u_{8,5}}{h^2} & f_{7,6}-frac{u_{8,6}}{h^2} & f_{7,7}-frac{u_{8,7}+u_{7,8}}{h^2} end{bmatrix}$.

Vamos a calcular uno a uno los elementos de la nueva malla mediante los dos métodos (restricción directa sobre la matriz anterior o restricción sobre las fronteras)  ya que con que salga alguno distinto ya podremos concluir su no equivalencia. Empezamos:

$latex frac{u_{1,2}^{2h}+u_{2,1}^{2h}-4u_{1,1}^{2h}}{(2h)^2} = frac{1}{16} [ f_{1,1}-frac{u_{1,0}+u_{0,1}}{h^2}+f_{1,3}-frac{u_{0,3}}{h^2}+f_{3,1}-frac{u_{3,0}}{h^2}+f_{3,3} + $

$latex + 2(f_{1,2}-frac{u_{0,2}}{h^2} + f_{2,1}-frac{u_{2,0}}{h^2} +f_{2,3} + f_{4,2} ) + 4 f_{2,2} ]$

Si primero aplicamos la restricción a las fronteras nos quedan:

$latex begin{bmatrix} frac{u_{0,1} + 2u_{0,2} + u_{0,3}}{4} & frac{u_{0,3} + 2u_{0,4} + u_{0,5}}{4} & frac{u_{0,5} + 2u_{0,6} + u_{0,7}}{4} end{bmatrix}$,

$latex begin{bmatrix} frac{u_{1,0} + 2u_{2,0} + u_{3,0}}{4} \ frac{u_{3,0} + 2u_{4,0} + u_{5,0}}{4} \ frac{u_{5,0} + 2u_{6,0} + u_{7,0}}{4} end{bmatrix}$ $latex ,,,,,,,$ $latex begin{bmatrix} frac{u_{1,8} + 2u_{2,8} + u_{3,8}}{4} \ frac{u_{3,8} + 2u_{4,8} + u_{5,8}}{4} \ frac{u_{5,8} + 2u_{6,8} + u_{7,8}}{4} end{bmatrix}$

$latex frac{1}{4}begin{bmatrix} u_{8,1} + 2u_{8,2} + u_{8,3} & u_{8,3} + 2u_{8,4} + u_{8,5} & u_{8,5} + 2u_{8,6} + u_{8,7} end{bmatrix}$,

y al calcular el primer término:

$latex frac{1}{16} [ f_{1,1} + f_{1,3} + f_{3,1} + f_{3,3} + 2(f_{1,2} + f_{2,1} + f_{2,3} + f_{4,2}) + 4 f_{2,2} ]$

que combinado con las fronteras, tenemos:

$latex frac{u_{1,2}^{2h}+u_{2,1}^{2h}-4u_{1,1}^{2h}}{(2h)^2} = frac{1}{16} [ f_{1,1} + f_{1,3} + f_{3,1} + f_{3,3} +$

$latex + 2(f_{1,2} + f_{2,1} + f_{2,3} + f_{4,2}) + 4 f_{2,2} ] – $

$latex – frac{1}{4}frac{u_{0,1} + 2u_{0,2} + u_{0,3} + u_{1,0} + 2u_{2,0} + u_{3,0}}{(2h)^2}$

y que es lo mismo que habíamos obtenido…

Llego ya un tiempo implementando los diferentes esquemas multigrid en C++ y, tras un periodo de larga oscuridad, por fin aparece un poco de luz: ayer , por primera vez, pareció funcionar todo (en una dimensión y con fronteras Dirichlet, claro… Aunque todo lo demás ya está también implementado).

En el primer nivel necesitamos un esquema de relajación. Esta implementado un weighted Jacobi, por su, a mi humilde entender, natural paralelización en CUDA (y que cuando $latex w=1$ tenemos un Jacobi), un SOR (sobrerelajación, que con $latex w=1$ deriva en un Gauss-Seidel), que suele ser la opción de siempre y un weighted red-black Gauss-Seidel, de cara a paralelizaciones tradicionales (OpenMP, MPI).

En el segundo nivel tenemos los algoritmos de corrección con esquemas V-Cycle, en iterativo, y $latex mu-$Cycle, en recursivo.

Finalmente, tenemos el FMG en el último nivel, tambien en iterativo y en recursivo. En total, 12 combinaciones posibles (en realidad, doce para una dimensión, doce para dos dimensiones y doce para tres dimensiones)

A ver que tal se dan hoy las ejecuciones en dos dimensiones.

Tags: , , , , , , ,

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.

Tags: ,

Por unos motivos o por otros, aunque esencialmente por el tiempo que me robaba escribir  un post a diario, mi querido blog ha quedado inactivo durante una temporada.

Dado que es de uso casi personal, no creo que mis lectores esten muy molestos. Sin embargo, vamos a retomar esta buena costumbre de poner por escrito aquellas cosas que he ido trabajando, aunque no se pueda diariamente.

Para empezar, un reblog de la última entrada de Terry Tao, el genio matemático actual por excelencia, el Mozart de las matemáticas. En éste, formaliza matemáticamente, de dos maneras diferentes, el análisis dimensional de la física. Últimamente a dedicado una pequeña parcela de su tiempo a la física, lo cual es de agradecer cuando se comenta que una de las ocupaciones de los principales matemáticos de cada área es tratar de interesar a Tao en su área.

Y para terminar Jason Padgett, un savant matemático accidental (después de sufrir golpes en la cabeza durante un atraco) que ahora ve la realidad a lo fracta, de esta manera. Y, si las fuentes son ciertas, no es cuento, pues se le han realizado escaneres cerebrales que muestran actividad cognitiva en areas que normalmente no utilizamos y que en su caso compensan la actividad de las areas dañadas. Curiosos todos su diagramas: $latex pi$, dualidad onda-partícula, relatividad, etc.

Tags: , , , , ,

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

Restricciones:

Interpolaciones:

 

Queremos resolver el sistema $latex 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 $latex Omega^{h}$ a una malla con un tamaño de intervalo $latex h$):

  • Relajación sobre $latex Aboldsymbol{u} =boldsymbol{ f}$ en una malla muy gruesa para determinar una primera aproximación para la siguiente malla mas fina.
  • Relajación sobre $latex Aboldsymbol{u} =boldsymbol{ f}$ en la malla $latex Omega^{4h}$ para obtener un valor inicial para la malla $latex Omega^{2h}$.
  • Relajación sobre $latex Aboldsymbol{u} =boldsymbol{ f}$ en la malla $latex Omega^{2h}$ para obtener un valor inicial para la malla $latex Omega^{h}$.
  • Relajación sobre $latex Aboldsymbol{u} =boldsymbol{ f}$ en la malla $latex 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, $latex Aboldsymbol{e} = boldsymbol{r} = boldsymbol{f} – Aboldsymbol{v}$, para relajar sobre el error:

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

La

Ya comentamos en el post anterior que formato tienen los sistemas $latex Au=f$ que surgen a partir de la discretización de un problema de Poisson y que tendremos que resolver. Esta resolución la podriamos hacer de manera directa pero son mucho mas eficaces los métodos iterativos.

Se puede demostrar que la manera mas sencilla para obtener los algoritmos de los diferentes métodos iterativos es pensar sencillamente en la aproximación por diferencias finitas, despejando la variable central i calcularla en función de los vecinos del paso de tiempo anterior para Jacobi, introduciendo un peso para ponderar el paso anterior de la misma variable con el calculado por Jacobi para relajación y, finalmente, utilizar no las varialbles del paso anterior sino tambien las ya calculadas del paso actual para Gauss-Seidel.

Además, el hecho de trabajar con $latex n$ dimensiones solo significa colocar nuestra expresión en $latex n$ bucles anidados. Así de sencillo…

En el caso de $latex 1D$, tendriamos:

[sourcecode languaje=”cpp”]
for (int i=1; i<nx-1; i++) {

//Jacobi
vj[t+1][i] = 0.5 * ( vj[t][i-1] + vj[t][i+1] + h*h*f[i] );

//weighted Jacobi
vwj[t+1][i] = vwj[t][i] + w * ( vwj[t][i] – vj[t+1][i] );

//Gauss-Seidel
vgs[i] = 0.5 * ( vgs[i-1] + vgs[i+1] + h*h*f[i] );

//weighted Gauss-Seidel
vwgs[i] = vwgs[i] + w * ( vwgs[i] – vgs[i] );

}
[/sourcecode]

Tags: , , ,

En $latex 1D$ tenemos:

$latex u_{xx} = f(x)$,

con $latex 0 < x < 1$ y $latex u(0)= u(1)=0$. Particionamos el dominio en $latex n$ subintervalos introduciendo los puntos $latex x_i = ih$ con $latex h=frac{1}{n}$. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas (a partir de Taylor alrededor de $latex x_i$):

$latex frac{partial^2}{partial x^2} u(x_i) = frac{u_{i-1} -2 u_i + u_{i+1}}{h^2} – frac{h^2}{12} frac{partial^4}{partial x^4}u(xi_i)$

con $latex xi_i in (x_{i-1}, x_{i+1})$ por lo que:

$latex frac{v_{i-1}-2v_i+v_{i+1}}{h^2} = f_i$.

con un error de truncamiento local del orden de $latex O(h^2)$ y donde utilizamos la $latex v$ para hablar de la aproximación a la $latex u$. Si lo escribimos de manera matricial, tenemos:

$latex frac{1}{h^2} begin{bmatrix} -2 & 1 & & & & & \ 1 & -2 & 1 & & & & \ & 1 & -2 & 1 & & & \ & & cdot & cdot & cdot & & \ & & & cdot & cdot & cdot & \ & & & & 1 & -2 & 1 \ & & & & & 1 & -2end{bmatrix} begin{bmatrix} v_1 \ v_2 \ v_3 \ cdot \ cdot \ v_{n-2} \ v_{n-1} end{bmatrix} = begin{bmatrix} f_1 \ f_2 \ f_3 \ cdot \ cdot \ f_{n-2} \ f_{n-1} end{bmatrix}$

Existen otros esquemas de mayor orden para la aproximación, por ejemplo el esquema central de cuarto orden:

$latex frac{partial^2}{partial x^2} u(x_i) = frac{-u_{i-2}+16u_{i-1}-30u_{i}+16u_{i+1}-u_{i+2}}{12 h^2} + frac{h^4}{90} frac{partial^6}{partial x^6}u(xi_i)$

con $latex xi_i in (x_{i-2},x_{i+2})$.

En $latex 2D$ tenemos:

$latex u_{xx}+u_{yy} = f(x,y)$,

con $latex 0 < x < 1$, $latex 0 < y < 1$ y $latex u=0$ en el cuadrado unidad. Particionamos el dominio en $latex n_x times n_y$ subintervalos introduciendo los puntos $latex (x_i,y_j) = (ih_x, jh_y)$ con $latex h_x=frac{1}{n_x}$ y $latex h_y = frac{1}{n_y}$. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas:

$latex frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{h_x^2}+frac{v_{i,j-1}-2v_{i,j}+v_{i,j+1}}{h_y^2} = f_{i,j}$

con $latex v_{0,j} = v_{n_x,j} = v_{i,0} = u_{i,n_y} = 0$. Tenemos un error local de truncamiento del orden de $latex O(h_x^2+h_y^2)$.

Si hacemos $latex h_x = h_y = h$ nos queda una plantilla, que representa un operador que interactua localmente en la malla, como esta:

$latex A = frac{1}{h^2} begin{bmatrix} & 1 & \ 1 & -4 & 1 \ & 1 & end{bmatrix}$

En $latex 3D$ tenemos:

$latex u_{xx} + u_{yy} + u_{zz} = f(x,y,z)$,

con $latex 0 < x, y, z < 1$ y $latex u=0$ en la frontera del cubo unitario. Procediento de la misma manera que antes tenemos:

$latex frac{v_{i-1,j,k}-2v_{i,j,k}+v_{i+1,j,k}}{h_x^2}+frac{v_{i,j-1,k}-2v_{i,j,k}+v_{i,j+1,k}}{h_y^2} + frac{v_{i,j,k-1}-2v_{i,j,k}+v_{i,j,k+1}}{h_z^2}= f_{i,j,k}$

En este caso, la manera mas sencilla de representar el operador local es mediante el tensor:

$latex A = frac{1}{h^2} a_{ijk}$ con $latex a_{0jk}=begin{bmatrix} & & \ & 1 & \ & & end{bmatrix}, a_{1jk}=begin{bmatrix} & 1 & \ 1 & -6 & 1 \ & 1 & end{bmatrix}$ y $latex a_{2jk} = begin{bmatrix} & & \ & 1 & \ & & end{bmatrix}$

si $latex h_x = h_y = h_z = h$.

Tags: , ,

Aquí está el artículo donde aparece el nuevo esquema en el que el sistema se desacopla de manera jerárquica:

(1) Conocidas las cantidades hidrodinámicas conservadas, resolver:

$latex Delta X^i + frac{1}{3} mathcal{D}^i mathcal{D}_j X^j = 8 pi f^{ij} S_j^*$

para encontrar

$latex hat{A}^{ij} approx (LX)^{ij} = mathcal{D}^i X^j + mathcal{D}^j X^i – frac{2}{3} mathcal{D}_k X^k f^{ij}$.

(2) Resolver la ecuación:

$latex Delta psi = -2 pi psi^{-1} E^{*} – psi^{-7} frac{ f_{il} f_{jm} hat{A}^{lm} hat{A}^{ij} }{8}$

para encontrar $latex psi$, donde la unicidad local ahora esta garantizada. Podemos calcular $latex S^*$ de manear consistente.

(3) Resolver la ecuación:

$latex Delta(psi N) = 2 pi N psi^{-1} (E^* + 2 S^*) + N psi^{-7} frac{7 f_{il} f_{jm} hat{A}^{lm} hat{A}^{ij} }{8}$

para $latex N psi$, una ecuación lineal donde podemos aplicar el principio del máximo con lo que, con las codiciones de contorno apropiadas, se sigue la unicidad y existencia.

(4) Finalmente, resolver:

$latex Delta beta^i + frac{1}{3} mathcal{D}^i ( mathcal{D}_j beta^j ) = D_j( 2 N psi^{-6} hat{A}^{ij} )$

para encontrar $latex beta^i$.

Además, en este otro artículo, presentan una manera de reducir una ecuación elíptica vectorial, un complicado sistema acoplado de PDEs, a un conjunto de ecuaciones Poisson escalares desacopladas. Para el caso del shift, la $latex beta$ anterior, por ejemplo, en coordenadas esféricas, tendríamos:

(1) Resolver ecuación:

$latex Delta mu = mu_S$

que corresponde a la parte toroidal, para la resolución de la parte angular se introducen un potencial poloidal $latex eta$ y  un potencial toroidal $latex mu$ de manera que $latex boldsymbol{beta} = $, y está desacoplada del resto para obtener $latex mu$.

(2) Resolver la también desacoplada ecuación para la divergencia (de $latex boldsymbol{beta}$ respecto de la conexión plana $latex mathcal{D}$):

$latex Delta Theta = frac{3}{4} mathcal{D}_{hat{k}} S(boldsymbol{beta}^{hat{k}})$.

(3) Obtener $latex beta^r$ a partir de una de las siguiente ecuaciones:

(i) $latex frac{partial^2 beta^r}{partial r^2} + frac{4}{r} frac{partial beta^r}{partial r} + frac{2 beta^r}{r^2} + frac{1}{r^2}Delta_{theta varphi} beta^r = S(boldsymbol{beta})^r – frac{1}{3} frac{partial Theta}{partial r} + frac{2}{r} Theta$

(ii) $latex Delta chi = r S(boldsymbol{beta})^r – frac{r}{3} frac{partial Theta}{partial r} + 2 Theta $, donde $latex chi = r beta^r$

(4) Deducir $latex eta$ de una de las siguientes ecuaciones:

(i) $latex Delta_{theta varphi} eta = r Theta – r frac{partial beta^r}{delta r} – 2 beta^r$, que tiene la ventaja de que solo requiere una división por $latex -l (l+1)$ de los coeficientes de la expansión por armónicos esféricos pero la desventaja de que utiliza la derivada radial de $latex beta^r$ que puede tener problemas con el orden.

(ii) $latex Delta eta = eta_S – frac{2 beta^r}{r^2} – frac{1}{3} frac{Theta}{r}$, que requiere la resolución de otra ecuación de Poisson adicional.

Tags: , , , , ,

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

$latex 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$

$latex 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 $latex rho_a = sum_b m_b W_{ab}$ y $latex P_a$ de las EoS.

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

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

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

$latex 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:

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

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

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

para la primera y:

$latex 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}))$

$latex 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}))$

$latex 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:

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

$latex = 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)$

Obviamente, en los kernels cúbico, cuártico y quíntico, al llegar a derivadas de orden superior, tenemos expresiones del tipo $latex alpha(boldsymbol{q}) . beta(boldsymbol{q})$ y lo que hemos hecho, de manera incorrecta, es:

$latex frac{partial}{partial q_j} big ( alpha(boldsymbol{q}) . beta(boldsymbol{q}) big ) = alpha(boldsymbol{q}) . frac{partial}{partial q_j} beta(boldsymbol{q})$,

por lo que nos falta sumarle:

$latex frac{partial}{partial q_j} alpha(boldsymbol{q}) . beta(boldsymbol{q})$.

Así pues, para el kernel Gaussiano si es correcta la formula genérica, pero no para el resto. De todas maneras, si es válida para las primeras derivadas.

A continuación una animación con:

$latex W_G(boldsymbol{q}), frac{partial}{partial q_1}W_G(boldsymbol{q}), frac{partial}{partial q_2}W_G(boldsymbol{q}), frac{partial^2}{partial q_1^2}W_G(boldsymbol{q}), frac{partial^2}{partial q_1 partial q_2}W_G(boldsymbol{q}), frac{partial^2}{partial q_2^2}W_G(boldsymbol{q}),$

$latex W_3(boldsymbol{q}), frac{partial}{partial q_1}W_3(boldsymbol{q}), frac{partial}{partial q_2}W_3(boldsymbol{q}),$

$latex W_4(boldsymbol{q}), frac{partial}{partial q_1}W_4(boldsymbol{q}), frac{partial}{partial q_2}W_4(boldsymbol{q}),$

$latex W_5(boldsymbol{q}), frac{partial}{partial q_1}W_5(boldsymbol{q}), frac{partial}{partial q_2}W_5(boldsymbol{q})$

[youtube http://www.youtube.com/watch?v=LEKKtUjuk7o?rel=0&w=420&h=315]

Tags: , , ,

Vamos a suponer que estamos en dimensión $latex n$ y miraremos como quedan los diferentes kernels y sus derivadas. Denotaremos con $latex boldsymbol{q}$ a $latex (q_1,ldots,q_n)$ y $latex q = sqrt{sum_{i=1}^n q_i^2}$ será su módulo.

Para poder hablar de las derivadas parciales en espacios de esta dimensión vamos a introducir la notación multi-índice de Schwartz. Un multi-índice de dimensión $latex n$ es una $latex n$-tupla de enteros no negativos $latex alpha = (alpha_1,ldots,alpha_n)$ con $latex alpha_j in mathbb{Z}_+$ para $latex j= 1ldots n$, de manera que $latex alpha in mathbb{Z}_+^n$. Llamamos  módulo de $latex alpha$ a:

$latex |alpha| = sum_{i=1}^n alpha_i$

Además, podemos definir

$latex alpha + beta := (alpha_1 + beta_1, ldots ,alpha_n + beta_n) $

y

$latex q^alpha:=Pi_{i=1}^n q_i^{alpha_i}$.

Finalmente, definimos el operador diferencial para $latex alpha in mathbb{Z}_+^n$ como:

$latex D^alpha:=frac{partial^{|alpha|}}{partial q_1^{alpha_1} partial q_2^{alpha_2} ldots partial q_n^{alpha_n}}$

kernel Gaussiano:

$latex W(boldsymbol{q},h) = frac{sigma}{h^n} e^{-(frac{q}{h})^2}$

$latex frac{partial}{partial q_j}W(boldsymbol{q},h) = D^{e_j}W(boldsymbol{q},h) = frac{sigma}{h^n}frac{-2}{h^2} e^{-(frac{q}{h})^2}q_j$

$latex frac{partial^2}{partial q_j partial q_k}W(boldsymbol{q},h)=D^{e_j+e_k}W= frac{sigma}{h^n}(frac{-2}{h^2})^2 e^{-(frac{q}{h})^2} q^{e_j+e_k}=frac{sigma}{h^n}(frac{-2}{h^2})^2 e^{-(frac{q}{h})^2}q_jq_k$

$latex D^alpha W(boldsymbol{q},h) = frac{sigma}{h^n}(frac{-2}{h^2})^{|alpha|} e^{-(frac{q}{h})^2}q^alpha$

kernel cúbico:

$latex W(boldsymbol{q},h) = frac{sigma}{h^n} begin{cases} frac{1}{4}(2-q)^3-(1-q)^3 mbox{ si } 0 leq q < 1 \ frac{1}{4}(2-q)^3 mbox{ si } 1 leq q < 2 \ 0 mbox{ si } q geq 2 end{cases}$

Recordemos que $latex q = sqrt{ sum_{i=1}^n q_i^2 }$, por lo que $latex frac{partial}{partial q_j}q = frac{2 q_j}{2 sqrt{sum_{i=1}^{n} q_i^2}} = frac{q_j}{q}$. Así pues:

$latex frac{partial}{partial q_j}W(boldsymbol{q},h) = frac{sigma}{h^n}(-1)frac{3}{q} q_j begin{cases} frac{1}{4}(2-q)^2-(1-q)^2 mbox{ si } 0 leq q < 1 \ frac{1}{4}(2-q)^2 mbox{ si } 1 leq q < 2 \ 0 mbox{ si } q geq 2 end{cases}$

$latex frac{partial^2}{partial q_j partial q_k}W(boldsymbol{q},h) = frac{sigma}{h^n}frac{3.2}{q^2} q_j q_k begin{cases} frac{1}{4}(2-q)-(1-q) mbox{ si } 0 leq q < 1 \ frac{1}{4}(2-q) mbox{ si } 1 leq q < 2 \ 0 mbox{ si } q geq 2 end{cases}$

$latex frac{partial^3}{partial q_j partial q_k partial q_l}W(boldsymbol{q},h) = frac{sigma}{h^n}(-1)frac{3.2.1}{q^3} q_j q_k q_l begin{cases} frac{1}{4}-1 mbox{ si } 0 leq q < 1 \ frac{1}{4} mbox{ si } 1 leq q < 2 \ 0 mbox{ si } q geq 2 end{cases}$

Por lo tanto, para $latex 0 leq |alpha| leq 3$ (si $latex |alpha|>3$ entonces $latex W equiv 0$) podemos escribir:

$latex D^alpha W(boldsymbol{q},h) = $

$latex =frac{sigma}{h^n}(-1)^{|alpha|}frac{frac{3!}{(3-|alpha|)!}}{q^{|alpha|}} q^alpha begin{cases} frac{1}{4}(2-q)^{3-|alpha|}-(1-q)^{3-|alpha|} mbox{ si } 0 leq q < 1 \ frac{1}{4}(2-q)^{3-|alpha|} mbox{ si } 1 leq q < 2 \ 0 mbox{ si } q geq 2 end{cases}$

kernel cuártico:

 $latex W(boldsymbol{q},h) = frac{sigma}{h^n} begin{cases} (frac{5}{2}-q)^4-5(frac{3}{2}-q)^4+10(frac{1}{2}-q)^4 mbox{ si } 0 leq q < frac{1}{2} \ (frac{5}{2}-q)^4-5(frac{3}{2}-q)^4 mbox{ si } frac{1}{2} leq q < frac{3}{2} \ (frac{5}{2}-q)^4 mbox{ si } frac{3}{2} leq q < frac{5}{2} \ 0 mbox{ si } q geq frac{5}{2} end{cases}$

Por tanto, para $latex |alpha| = 0 ldots 4$:

$latex D^alpha W(boldsymbol{q},h) =$

$latex =frac{sigma}{h^n}(-1)^{|alpha|}frac{frac{4!}{(4-|alpha|)!}}{q^{|alpha|}} q^alpha begin{cases} (frac{5}{2}-q)^{4-|alpha|}-5(frac{3}{2}-q)^{4-|alpha|}+10(frac{1}{2}-q)^{4-|alpha|} mbox{ si } 0 leq q < frac{1}{2} \ (frac{5}{2}-q)^{4-|alpha|}-5(frac{3}{2}-q)^{4-|alpha|} mbox{ si } frac{1}{2} leq q < frac{3}{2} \ (frac{5}{2}-q)^{4-|alpha|} mbox{ si } frac{3}{2} leq q < frac{5}{2} \ 0 mbox{ si } q geq frac{5}{2} end{cases}$

kernel quíntico:

$latex W(boldsymbol{q},h) = frac{sigma}{h^n} begin{cases} (3-q)^5-6(2-q)^5+15(1-q)^5 mbox{ si } 0 leq q < 1 \ (3-q)^5-6(2-q)^5 mbox{ si } 1 leq q < 2 \ (3-q)^5 mbox{ si } 2 leq q < 3 \ 0 mbox{ si } q geq 3 end{cases}$

Por tanto, para $latex |alpha| = 0 ldots 5$:

$latex D^alpha W(boldsymbol{q},h) =$

$latex =frac{sigma}{h^n}(-1)^{|alpha|}frac{frac{5!}{(5-|alpha|)!}}{q^{|alpha|}} q^alpha begin{cases} (3-q)^{5-|alpha|}-6(2-q)^{5-|alpha|}+15(1-q)^{5-|alpha|} mbox{ si } 0 leq q < 1 \ (3-q)^{5-|alpha|}-6(2-q)^{5-|alpha|} mbox{ si } 1 leq q < 2 \ (3-q)^{5-|alpha|} mbox{ si } 2 leq q < 3 \ 0 mbox{ si } q geq 3 end{cases}$

Tags: , , , , ,

Volvemos a tratar los kernels pero teniendo en cuenta que nuesta anterior $latex q$ de hecho es $latex frac{boldsymbol{||r||}}{h}$, por lo que en $latex 1D$, con $latex h=1$, tenemos $latex sqrt{x^2} = |x|$, en $latex 2D$ tenemos $latex sqrt{x^2+y^2}$ y en $latex 3D$ tenemos $latex sqrt{x^2+y^2+z^2}$.  En realidad, como lo que tenemos es $latex ||r-r’||$, habrá cosas del estilo de:

$latex sqrt{(x-x’)^2+(y-y’)^2+(z-z’)^2}$.

En los siguientes enlaces se encuentran información sobre cada uno de los kernels especificados:

  1. funciones definidas a trozos y sus gráficas
  2. kernel Gaussiano
  3. kernel cúbico
  4. kernel cuártico
  5. kernel quíntico

Dada una partícula $latex a$ situada en $latex r_a = (x_a,y_a)$ tenemos allí el kernel $latex W_a(r-r_a,h)$. Dada otra partícula vecina $latex b$ de posición $latex r_b = (x_b, y_b)$, el punto

$latex (x_b, y_b,W_a(r_b-r_a,h) = W_{ab} = W_{ba})$

queda sobre el kernel:

Tags: , , , ,

En palabra de Terry Tao:


En la física, el espacio de fases es un concepto que unifica la mecánica clásica (Hamiltoniana) con la mecánica cuántica; en matemáticas, el espacio de fases es un concepto que unifica la geometría simpléctica con el análisis armónico y las PDE.

En mecánica clásica, el espacio de fases es el espacio de todas las posibles configuraciones de un sistema: no solo las posiciones $latex q$ de todos los objetos del sistema, sino también sus momentos $latex p$. Matemáticamente, el espacio de configuraciones puede definirse como una variedad $latex M$ de manera que, para cada posicion $latex q in M$, los momentos $latex p$ toman valores en el espacio  cotangente $latex T_q^*M$. De esta manera, el espacio de fases puede verse de manera natural como el fibrado cotangente:

$latex T^*M:=bigsqcup_{q in M} T_q^*M$,

y, como ya vimos en ese mismo post, si $latex dim M = n$ entonces $latex dim T^*M = dim TM = 2n$, es decir, que el espacio de fases siempre va a tener dimensión par.

Tags: , , ,

We try to write the formulas in its most general form and fix

First of all, in the simplest form of the HD equations we have

$latex frac{dboldsymbol{v}_a}{dt} = -sum_b m_b (frac{P_a}{rho_a^2} + frac{P_b}{rho_b^2}) nabla_a W_{ab}$ (momentum)

where:

$latex rho_a = sum_b m_b W_{ab}$ (density),

$latex frac{de_a}{dt} = 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}$ (density energy),

$latex P_a$ and $latex P_b$

are obtained from the EoS and $latex frac{dboldsymbol{r}_a}{dt} = boldsymbol{v}_a$.

In SHD we have:

$latex d$

Finally, the GHD equations are:

$latex d$

Tags: , , ,

A continuación se muestra el fragmento de código de la app RadSol que resuelve por radicales la ecuación cuártica (pulsar para ampliar):

La idea, como se puede apreciar en el código, es que resuelve una ecuación cuadrática en la que aparece la solución real de la cúbica resolvente.

Para la resolvente, dada la ecuación de tercer grado $latex x^3 + bx^2 + cx+d = 0$, la fórmula de las soluciones es:

$latex -frac{b}{3}+frac{t}{3}+frac{b^2-3c}{3t}$ con $latex t = E^{frac{1}{3}}$

donde:

$latex E=frac{9bc – 2b^3 -27d + sqrt{D}}{2}$ y $latex D = (9bc -2b^3-27d)^2 + 4(3c-b^2)^3$

Tags: , , , ,

Construimos una malla de $latex 100times 100$ puntos de nuestros kernels en $latex 2D$. La primera animación contiene el kernel Gaussiano en el primer frame, el cúbico en el segundo, el cuártico en el tercero y el quíntico en el último:

[youtube http://www.youtube.com/watch?v=toa1Asoo600?rel=0&w=420&h=315]

En la siguiente animación mostramos los puntos de la primera derivada de los kernels:

[youtube http://www.youtube.com/watch?v=VMnyetAVyN0?rel=0&w=420&h=315]

Finalmente, la animación con puntos de las segundas derivadas de los kernels de nuestra aplicación:

[youtube http://www.youtube.com/watch?v=YZiuoCBrfOQ?rel=0&w=420&h=315]

Tags: , , , , ,

En $latex 2D$ nos encontramos las siguientes gráficas.

En primer lugar, dos imagenes, cualitativas y parciales, del Kernel Gaussiano y su derivada (en blanco opaco y hueco respectivamente) junto al Kernel cúbico y su derivada (en azul), y a continuación dos del Kernel Gaussiano y su derivada junto al Kernel quíntico y su derivada (en rojo):

Tags: ,

A continuación tenemos diferentes funciones kernel, el linea contínua, junto a su primera y segunda derivada, en lineas discontinuas de mayor y menor longitud respectivamente.

En color negro el kernel gaussiano, su primera derivada y su segunda derivada:

en azul el cúbico:

en verde el cuártico:

y en rojo el quíntico:

$latex y$ representa $latex f(q)$, $latex x$ representa $latex q$ con $latex r/h$ y estamos en $latex 1D$.

Tags: ,

Desde un punto de vista muy computacional, y de las ideas, donde las partículas en los métodos sin malla no son mas que el equivalente a los puntos de interpolación de los métodos mallados, unos libres y otros fijos, y dado un volumen en el que pretendemos modelar diferentes entidades como fluidos y campos, parece lógico utilizar puntos de interpolación libres para los primeros pero puntos de interpolación fijos para los segundos, pues mientras que habitualmente éstos no ocuparan todo el volumen modelado, aquellos si lo haran…

Parece lógico, o igual no tanto…  Lo ideal sería poder combinar ambos, ¿no?

Tags: , , ,

Como ya comentamos, de la tesis de Bauswein, adoptando la foliación $latex 3+1$ del espacio-tiempo la métrica queda:

$latex ds^2 = (- alpha^2 + beta_i beta^i) dt^2 + 2 beta_i dx^i dt + gamma_{ij} dx^i dx^j$

En la aproximación CFC resolvemos repetidamente el problema de valor inicial. De acuerdo con esta aproximación, la parte espacial de la métrica se puede escribir como:

$latex gamma_{ij} = psi^4 delta_{ij}$

donde $latex psi$ es el factor conforme (una transformación conforme preserva los ángulos. En geometría Riemanniana, dos métricas de Riemann $latex g$ y $latex h$ sobre una variedad $latex M$ son conformemente equivalentes si $latex g=uh$ para alguna función positiva $latex u$ sobre $latex M$. La función $latex u$ es el factor conforme).

De esta manera, las ecuaciones de Einstein, asumiendo $latex K := tr(K_{ij}) = K_i^i =0$, se reducen al sistema de cinco PDE elipticas no lineales acopladas:

$latex Delta psi = -2 pi psi^5 E – frac{1}{8} psi^5 K_{ij}K^{ij}$

$latex Delta(alpha psi) = 2 pi alpha psi^5 (E + 2S) + frac{7}{8} alpha psi^5 K_{ij}K^{ij}$

$latex Delta beta^i + frac{1}{3}partial^i partial_j beta^j = 16 pi alpha rho W + 2 psi^{10} K^{ij} partial_j (frac{alpha}{psi^6}) =: S_beta$

donde $latex E = rho h W^2 – P$, $latex S = rho h (W^2 -1) + 3P$ y

$latex K_{ij} = frac{psi^4}{2 alpha} (delta_{il} partial_j beta_l + delta_{jl} partial_i beta^l – frac{2}{3} delta_{ij} partial_k beta^k )$

que podemos escribir de manera mas compacta como:

$latex Delta B^i = S_beta$

$latex Delta chi = partial_i B^i$

si definimos $latex beta^i = B^i – frac{1}{4} partial_i chi$ y que es un sistema tipo Poisson que puede ser resuelto iterativamente hasta la convergencia con un método multigrid.

Las condiciones en la frontera se dan mediante desarrollo multipolar () de los terminos fuente, que son no compactas, hasta el armónico quadrupolar.

Tags: , , , ,

Finalmente, tenemos ya varias clases abstractas que instanciamos mediante la creación de algunas de sus clases concretas hijas en función de un fichero de configuración.

Algunas de estas clases abstractas son: OdeSolver, SpatialDecomposition, PdeSolver, Kernel, EoS, Equation, Display que pueden ser concretadas en alguna de sus clases herederas (por ejemplo SpatialDecomposition podria ser instanciada como AllPair, como LinkedList o como cualquiera que implementemos posteriormente).

A nivel de código, lo que nos encontramos es lo siguiente:

  1. En el fichero de configuración aparece algo como SpatialDecomposition=LinkedList o <SpatialDecomposition>LinkedList</SpatialDecomposition>, en función del formato que tenga nuestro fichero de configuración
  2. Este fichero lo lee una clase Configuracion que posteriormente es capaz de suministrar esta información.
  3. En algun lugar de nuestro código tenemos SpatialDecomposition *spatialDecomposition y una composición alternativa de acciones donde en una de las ramas, la que se corresponde a la información que nos suministra Configuración, con spatialDecomposition = new LinkedList(...)
  4. En la clase abstracta SpatialDecomposition tenemos un método virtual getNNP que implementan todos sus descendientes, de manera que en LinkedList tenemos su correspondiente implementación para esta clases concreta.
  5. Cuando realizamos la invocación del método spatialDecomposition->getNNP(...) (-> en lugar de . porque tenemos un puntero a la clase y no la clase) este siempre tendrá este formato independientemente de la clase concreta que tengamos en un momento determinado por debajo, por lo que si posteriormente decidimos utilizar AllPair, la llamada anterior queda de la misma manera.

Esto permite que exista una capa de abstracción de manera que el conjunto principal de clases trabaja invocando a métodos virtuales de clases abstractas que se transforman en llamadas a métodos concretos de la clases descendiente instanciada de manera trasparente.

Tags: , , , ,

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
...

Tags: , , , , , , ,

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 $latex h$.

En la siguiente animación se utiliza una LinkedList:

Tags: , , ,

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

Supongamos que tenemos $latex N$ partículas. Si $latex 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 $latex h$ (all-pair search). La complejidad del algorítmo es $latex 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 $latex 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 $latex h$ es constante para todas las partículas. La complejidad, con un número suficientemente pequeño de partículas por celda, es $latex O(N)$.

Finalmente, podemos trabajar con árboles. Los árboles permiten complejidades genéricas en búsquedas de $latex O(N log N)$ y permiten trabajar con $latex 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).

Tags: , , , , , , ,

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

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

La ley de Newton $latex 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:

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

que es:

$latex m_ax_a” = -k_1 (1-frac{l}{sqrt{x_a^2+y_a^2}})x_a – k_2frac{x’_a}{sqrt{x_a’^2+y_a’^2}}$

$latex m_ay_a” = -k_1 (1-frac{l}{sqrt{x_a^2+y_a^2}})y_a – k_2frac{y’_a}{sqrt{x_a’^2+y_a’^2}}$

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

$latex 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_2frac{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_2frac{v_a}{sqrt{u_a^2+v_a^2}} – m_ag \
end{array}
right .
$

En la siguiente simulación tomamos $latex m_a$ aleatoria para cada partícula $latex a$, por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos $latex l=1$, $latex k_1 = 10000$ y $latex 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, $latex u_{a_0}=0$, $latex v_{a_0} = -0.7$ y $latex v_{z_0} = 0$. El intervalo de tiempo es desde $latex t_i = 0$ hasta $latex t_f=10$ que se divide en $latex 400$ subintervalos.

El resultado es el siguiente:

[youtube http://www.youtube.com/watch?v=emirYBIds6o?rel=0&w=420&h=315]

Tags: , , ,

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.

Tags: , , , ,

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.

Tags: ,

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:

[youtube http://www.youtube.com/watch?v=-bgdirYegrc?rel=0&w=420&h=315]

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.

Tags: , ,

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

Tags: ,

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!

Tags: ,

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:

$latex rho_i = sum_j m_j W_{ij}$

con $latex W_{ij}:=W(boldsymbol{x}_i – boldsymbol{x}_j)$, o también:

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

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

Ellos proponen utilizar,

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

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

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

Tags: , ,

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:

$latex langle f rangle(boldsymbol{x}) := int f(boldsymbol{x’})W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

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

$latex langle frangle(boldsymbol{x}_i) – f(boldsymbol{x}_i) = frac{h^2_{ef}}{4} nabla^2 f(boldsymbol{x}) + O(h^4_{ef})$.

donde

$latex h^2_{ef}:=2 int x^2 W(boldsymbol{x},h)dboldsymbol{x}$,

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

También se tiene:

$latex langle frac{partial f}{partial x} rangle (boldsymbol{x}) = int frac{partial f(boldsymbol{x’})}{partial x’} W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

$latex = – int f(boldsymbol{x’}) frac{partial}{partial x’} W(boldsymbol{x’}-boldsymbol{x},h)dboldsymbol{x’}$

$latex = frac{partial}{partial x} int f(x’)W(x’-x,h)dx’$

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

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

entonces:

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

$latex 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:

$latex f_i:= langle f rangle (boldsymbol{x}_i) = int f(boldsymbol{x’})W(boldsymbol{x’}-boldsymbol{x}_i,h)dboldsymbol{x’}$

$latex = 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)dboldsymbol{x}’$

$latex sum_j Delta f_{i,j}$

donde:

$latex 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)dboldsymbol{x}’$.

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

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

y por otra, con:

$latex 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:

$latex ddot{boldsymbol{x}}:=int frac{dboldsymbol{v}(boldsymbol{x})}{dt}W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

$latex = – int frac{1}{rho(boldsymbol{x})} nabla P(x) W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

y para la energía

$latex dot{u}_i := int frac{du(boldsymbol{x})}{dt} W(boldsymbol{x}-boldsymbol{x}_i,h)dboldsymbol{x}$

$latex = 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) ] dboldsymbol{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.

Tags: , , , ,

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 $latex mathbb{R}^n$ o $latex 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.

Tags: , ,

Teorema de Laurent.

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

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

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

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

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

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

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

$latex 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 $latex f:U-A longrightarrow mathbb{C}$ una función analítica con $latex U subset mathbb{C}$ un abierto y $latex A$ el conjunto de singularidades aisladas de $latex f$. Diremos que $latex f$ es una función meromorfa en $latex U$.

Teorema (de los residuos)

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

Tags: , , , , , , ,

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 $latex t$ es su estado cuántico, que se representa por un vector (ket) de norma $latex 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 $latex A$ sobre un estado $latex |psi rangle$ solo puede ser uno de sus autovalores. La probabilidad de que al medir $latex A$, de descomposicón espectral $latex A=sum_i lambda_i Pi_i,, (lambda_i neq lambda_j)$, sobre un estado $latex |psi rangle$ obtengamos el resultado $latex lambda_i$ es:

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

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

Dualidad onda-particula

Incertidumbre de Heisenberg

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

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

$latex frac{Pi_i | phi rangle}{||Pi_i | phi rangle||}$

con probabilidad $latex ||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

$latex ihbar frac{partial}{partial t} | varphi(t) rangle = H | varphi (t) rangle$

donde $latex 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:

$latex [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 $latex X_i$ y $latex P_i$ evolucionan según las leyes clásicas pero con fuerzas medias.

Tags: , , , , , , , , , , ,

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.

Tags: , , , , , , , , , , , , , , , , ,

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 $latex Omega subset mathbb{C}$ un abierto simplemente conexo ($latex I(gamma,z) = 0$ si $latex z notin Omega$), $latex f in mathcal{H}(Omega)$ y $latex gamma$ un ciclo tal que $latex gamma^* subset Omega$. Entonces, para $latex z in Omega-gamma^*$:

$latex 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 $latex Omega subset mathbb{C}$ un abierto simplemente conexo ($latex I(gamma,z) = 0$ si $latex z notin Omega$), $latex f in mathcal{H}(Omega)$ y $latex gamma$ un ciclo tal que $latex gamma^* subset Omega$. Entonces, para $latex z in Omega-gamma^*$:

$latex int_gamma f(u)du = 0$

demostración:

Fórmula integral de Cauchy para las derivadas.

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

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

demostración:

Tags: , , ,

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 $latex H$ sobre un cuerpo $latex mathbb{K} = mathbb{R}, mathbb{C}$ en el que tenemos definido un producto escalar.

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

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

tal que es:

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

Ejemplos

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

Vamos a demostrar este último ejemplo. Probaremos, en primer lugar, que $latex ell^2$ es un espacio vectorial, y posteriormente, que $latex 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:

$latex langle cdot | cdot rangle : $

Espacios de Banach

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

$latex ||cdot||_E : E longrightarrow [0,+infty[$

que satisface las propiedades de:

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

Tags: , ,

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 $latex f$ una función analítica, es decir, $latex f = sum_{n=0}^infty a_n(z-z_0)^n$ con $latex z in D(z_0,R)$. Sea $latex 0 < r < R$ y

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

Entonces

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

demostración:

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

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

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

$latex 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:

$latex 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 $latex |int_gamma f(z) dz| leq L(gamma) , max { |f(z)| : z in gamma([a,b]) }$, nos queda:

$latex |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 $latex |frac{1}{2 pi i}| = frac{1}{2 pi}$ y

$latex 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:

$latex |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 $latex f in mathcal{H}(mathbb{C})$, $latex f$ no constante. Entonces $latex f$ no está acotada.

demostración:

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

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

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

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

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

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

$latex |a_n| leq frac{C}{r^n}$.

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

Lo que nos está diciendo el teorema de Liouville, por ejemplo, es que la función $latex 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 $latex p(z) = a_0 + a_1 z + a_2 z^2 + ldots + a_n z^n$ un polinomio no constante. Entonces la ecuación $latex p(z)=0$ tiene solución.

demostración:

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

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

$latex 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 $latex |a+b| geq |a|-|b|$, nos queda:

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

$latex ,Box$

Proposición:

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

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

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

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

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

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

Tags: , ,

Fórmula integral de Cauchy para las derivadas.

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

$latex 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 $latex f(xi) = frac{1}{2 pi i} int_{C(z_0,R)} frac{f(u)}{u-xi}du$ siempre que $latex |xi – z_0|<R$.

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

$latex 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:

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

Aplicaciones:

1. Calcular

$latex 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 $latex 0$. Tenemos que aislarla, por lo que escribimos la integral de la siguiente forma:

$latex 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 $latex f(z):=frac{e^z}{(1-z)^2}$. Ahora tenemos que $latex f(z)$ es holomorfa en todos los puntos excepto en el $latex 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

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

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

2. Calcular

$latex 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 $latex 1$ y tenemos que aislarlo para poder aplicar el teorema integral para derivadas. Escribimos

$latex 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 $latex f'(z) = frac{ze^z – e^z}{z^2}$ y $latex f'(1)=0$, donde $latex f(z):=frac{e^z}{z}$, nos queda que la integral vale $latex 0$.

3. Calcular

$latex 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 $latex frac{1}{z(z-1)^2}$ en fracciones simples:

$latex 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 $latex A=1=C, B=-1$, por lo que podemos escribir:

$latex 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 $latex f(z):=e^z$ que es holomorfa en todo el plano, a cada integral:

$latex I = 2 pi i (f(0) – f(1) + f'(1)) = 2 pi i (1 – e + e) = 2 pi i$.

Tags: ,

En el cáculo de este post, como $latex mathcal{L}_X dalpha = d mathcal{L}_X alpha$, y en este caso conocemos $latex alpha$, también podemos hacer:

$latex mathcal{L}_{[X,Y]}d alpha = dmathcal{L}_Z alpha = d big [ mathcal{L}_{(-4 frac{partial}{partial x} – 12y frac{partial}{partial z})} (2xy,dz – z , dx big ])$

Calculamos, en primer lugar, $latex mathcal{L}_Z alpha$:

$latex mathcal{L}_Z (2xy , dz) – mathcal{L}_Z (z , dx)$

Por una parte, tenemos:

$latex mathcal{L}_Z (2xy , dz) = (mathcal{L}_Z 2xy) , dz + 2xy , mathcal{L}_Z dz =$

$latex = (-4 frac{partial}{partial x} 2xy – 12y frac{partial}{partial z} 2xy ) , dz + 2xy , d(-4 frac{partial}{partial x} z – 12y frac{partial}{partial z} z) =$

$latex = -24xy , dy – 8y , dz$

Y por otra:

$latex mathcal{L}_Z (z , dx) = mathcal{L}_Z z , dx + z , mathcal{L}_Z dx =$

$latex = (-4 frac{partial}{partial x} z – 12y frac{partial}{partial z} z) , dx + z , d(-4 frac{partial}{partial x} x – 12y frac{partial}{partial z} x) = $

$latex = -12y , dx$

Por lo tanto, tenemos:

$latex 12y , dx -24xy ,dy – 8y , dz$

Finalmente, solo queda calcular:

$latex d big [ 12y , dx -24xy ,dy – 8y , dz big ] = $

$latex = 12 , dy wedge dx – 24 , d(xy) wedge dy – 8 , dy wedge dz = $

$latex 12 , dy wedge dx – 24 , (y , dx + x , dy) wedge dy – 8 , dy wedge dz = $

$latex -12(2y + 1) , dx wedge dy – 8 , dy wedge dz$

ya que $latex -24x , dy wedge dy = 0$ por la propiedad de que $latex d^2 = 0$, obteniendo así el mismo resultado que anteriormente.

Tags:

Ya nos apareció la derivada de Lie. Con los datos de este post, ¿Cómo calcularíamos $latex mathcal{L}_{[X,Y]} beta$ con $latex beta := d alpha$?

Primero necesitamos calcular el corchete de Lie de los campos dados:

$latex Z:=[X,Y] = -4 frac{partial}{partial x} -12y frac{partial}{partial z}$

y que es otro campo vectorial, a continuación necesitamos  la $latex 2$-forma resultante de calcular la diferencial exterior de la $latex 1$-forma::

$latex beta = 2y+1,dx wedge dz + 2x,dywedge dz$

y finalmente calcular la derivada de Lie de la forma respecto del campo.

Todos los cálculos se reduciran a saber aplicar la derivada de Lie a funciones, campos vectoriales y a la diferencial exterior de $latex 1-$formas sabiendo que es una derivación:

  1. $latex mathcal{L}_Xh = X(h)$
  2. $latex mathcal{L}_XY = [X,Y]$
  3. $latex mathcal{L}_X dalpha = d mathcal{L}_X (alpha)$

de manera que:

$latex mathcal{L}_Z beta = mathcal{L}_Z (2y+1,dx wedge dz + 2x,dywedge dz) = $

$latex mathcal{L}_Z (2y+1 , dx wedge dz) + mathcal{L}_Z (2x , dy wedge dz) = $

$latex mathcal{L}_Z(2y+1) , dx wedge dz + (2y +1) mathcal{L}_Z(dx) wedge dz + (2y+1) , dx wedge mathcal{L}_Z(dz) + $

$latex mathcal{L}_Z(2x) , dy wedge dz + 2x mathcal{L}_Z(dy) wedge dz + 2x , dy wedge mathcal{L}_Z(dz))$.

Para evitar errores, calculamos separadamente cada derivada de Lie:

$latex mathcal{L}_Z (2y+1) = Z(2y+1) = [X,Y](2y+1) = (-4 frac{partial}{partial x} -12y frac{partial}{partial z})(2y+1) = $

$latex = -4 frac{partial}{partial x}(2y+1) -12y frac{partial}{partial z}(2y+1) = 0$

$latex mathcal{L}_Z (dx) = d mathcal{L}_Z x = d([X,Y](x)) = d(-4 frac{partial}{partial x}x -12y frac{partial}{partial z}x) = d(-4)=0$

$latex mathcal{L}_Z(dz) = d mathcal{L}_Z z = d([X,Y](z) = d(-4 frac{partial}{partial x}z -12y frac{partial}{partial z}z)) = -12 dy$

$latex mathcal{L}_Z (2x) = [X,Y](2x) = -4 frac{partial}{partial x}2x -12y frac{partial}{partial z}2x = -8$

$latex mathcal{L}_Z (dy) = d mathcal{L}_Z y = d([X,Y](y)) = d(-4 frac{partial}{partial x}y -12y frac{partial}{partial z}y) = 0$

$latex mathcal{L}_Z(dz) = -12 dy$

Por lo que, finalmente, tenemos:

$latex -12(2y+1) , dx wedge dy – 24 , dy wedge dy – 8 , dy wedge dz$

y como $latex d^2 = 0$, nos queda:

$latex -12(2y+1) , dx wedge dy -8 , dy wedge dz$

Tags:

Realizaremos algunos cálculos típicos en variedades diferenciables. Para ello, vamos a suponer que tenemos el campo vectorial

$latex X = -4y frac{partial}{partial x} + 9x frac{partial}{partial y} + frac{partial}{partial z} in mathfrak{X}(mathbb{R}^3)$

y el campo vectorial

$latex Y=-frac{partial}{partial y} + 3x frac{partial}{partial z} in mathfrak{X}(mathbb{R}^3)$.

Tenemos tambien la $latex 1$-forma $latex alpha = 2xy dz – z dx$. En este caso trabajamos con la variedad $latex M = mathbb{R}^3$ (en física los escribiríamos $latex X(x,y,z) = (-4y,9x,1)$ o $latex vec{Y}=(0,-1,3x)$).

Empezaremos calculando un corchete de Lie. ¿Cómo queda el cálculo de $latex [X,Y]$? Si escribimos el corchere con los campos concretos nos queda:

$latex [X,Y] = [-4y frac{partial}{partial x} + 9x frac{partial}{partial y} + frac{partial}{partial z},-frac{partial}{partial y} + 3x frac{partial}{partial z}]$.

Lo primero que hacemos es aplicar que el corchete de Lie es bilineal, por lo que la expresión anterior queda:

$latex [-4y frac{partial}{partial x},-1 frac{partial}{partial y} ] + [-4y frac{partial}{partial x}, 3x frac{partial}{partial z}] + $

$latex + [9x frac{partial}{partial y},-1 frac{partial}{partial y}] + [9x frac{partial}{partial y}, 3x frac{partial}{partial z}] + $

$latex + [1 frac{partial}{partial z},-1 frac{partial}{partial y}] + [1 frac{partial}{partial z}, 3x frac{partial}{partial z}]$

Se puede demostrar que si $latex g,h in C^infty(M)$ y $latex X,Y in mathfrak{X}(M)$ entonces:

$latex [gX,hY] = gX(h)Y – hY(g)X + gh[X,Y]$,

por lo que, si la aplicamos a nuestro caso concreto, nos queda:

$latex -4y frac{partial(-1)}{partial x} frac{partial}{partial y} + 1 frac{partial (-4y)}{partial y} frac{partial}{partial x} + (-4y)(-1)[frac{partial}{partial x},frac{partial }{partial y}] + $

$latex -4y frac{partial(3x)}{partial x} frac{partial}{partial z} – 3x frac{partial (-4y)}{partial z} frac{partial}{partial x} + (-4y)(3x)[frac{partial}{partial x},frac{partial }{partial z}] + $

$latex 9x frac{partial(-1)}{partial y} frac{partial}{partial y} + 1 frac{partial (9x)}{partial y} frac{partial}{partial y} + (9x)(-1)[frac{partial}{partial y},frac{partial }{partial y}] + $

$latex 9x frac{partial(3x)}{partial y} frac{partial}{partial z} – 3x frac{partial (9x)}{partial z} frac{partial}{partial y} + (9x)(3x)[frac{partial}{partial y},frac{partial }{partial z}] + $

$latex 1 frac{partial(-1)}{partial z} frac{partial}{partial y} + 1 frac{partial (1)}{partial y} frac{partial}{partial z} + (1)(-1)[frac{partial}{partial z},frac{partial }{partial y}] + $

$latex 1 frac{partial(3x)}{partial z} frac{partial}{partial z} – 3x frac{partial (1)}{partial z} frac{partial}{partial z} + (1)(3x)[frac{partial}{partial z},frac{partial }{partial z}]$.

Realizando las derivadas parciales indicadas ($latex frac{partial (-1)}{partial x} =0, frac{partial (-4y)}{partial y} = 4$, …) y teniendo en cuenta que el corchete de Lie de campos coordenados es nula, nos queda:

$latex [X,Y] = -4 frac{partial}{partial x} -12y frac{partial}{partial z} in mathfrak{X}(mathbb{R}^3)$.

Vamos a calculara ahora $latex dalpha$. Para empezar, aprovechamos que es lineal, por lo que:

$latex dalpha = d(2xydz – zdx) = d(2xydz) – d(zdx)$

y ahora aplicamos la definición del operador diferencial exterior a cada sumando:

$latex d(2xy)wedge dz – d(z) wedge dx$

finalmente, como la diferencial exterior sobre funciones es la diferencial ordinaria, aplicando las propiedades distributiva, antisimétrica y $latex d^2=0$, nos queda:

$latex (2ydx + 2xdy) wedge dz – dz wedge dx = $

$latex 2y , dx wedge dz + 2x , dy wedge dz + dx wedge dz$

Por lo, finalmente, nos queda:

$latex dalpha = 2y+1,dx wedge dz + 2x,dywedge dz$

Tags: ,

En el artículo “MAGMA: a three-dimensional, Lagrangian magnetohydrodynamic code for merger applications” de S. Rosswog y D. Price comentan como introducir campos magnéticos en SPH. Las ecuaciones ya discretizadas quedan:

Ecuación de densidad:

$latex rho = sum_b m_b W(r-r_b,h)$

Ecuación del momento (conh: “grad-h” term, mag: magnetic force term, g: self-gravity and gravitational softening term)

$latex frac{d}{dt}v_{a,MHD} = frac{d}{dt} (v_{a,h}+v_{a,h,dis}+v_{a,g}+v_{a,mag}+v_{a,mag,dis})$

donde

$latex frac{d}{dt}v_{a,h} = -sum_b m_b big ( frac{P_a}{Omega_a rho_a^2} nabla_a W_{ab}(h_a)+ frac{P_b}{Omega_b rho_b^2} nabla_a W_{ab}(h_b) big )$

$latex frac{d}{dt}v_{a,h,dis} = $

$latex frac{d}{dt}v_{a,g} = -G sum_b m_b big [ frac{phi’_{ab}(h_a) + frac{phi’_{ab}(h_b)}{}}{2} big ]hat{e}_{ab}$

$latex -frac{G}{2} sum_b m_b big [ frac{zeta_a}{Omega_a} nabla_a W_{ab}(h_a) + frac{zeta_b}{Omega_b} nabla_a W_{ab}(h_b) big ]$

con

$latex phi’_{ab} = frac{partial phi}{partial|r_a-r_b|}$, $latex zeta_k := frac{partial h_k}{partial rho_k} sum_b m_b frac{partial phi_{kb}(h_k)}{partial h_k}$

y

$latex Omega_a := big ( 1- frac{partial h_a}{partial rho_a} cdot sum_b m_b frac{partial}{partial h_a} W_{ab}(h_a) big )$

$latex frac{d}{dt}v_{a,mag} = – sum_b frac{m_b}{mu_0} big { frac{B_a^2 / 2}{Omega_a rho_a^2} nabla_a W_{ab}(h_a) + frac{B_b^2 / 2}{Omega_b rho_b^2} nabla_a W_{ab}(h_b) }$

$latex + sum_b frac{m_b}{mu_0} big { frac{B_a(B_a cdot overline{nabla_a W_{ab}}) – B_b(B_b cdot overline{nabla_a W_{ab}})}{rho_a rho_b} big }$

con

$latex overline{nabla_a W_{ab}} = frac{1}{2} big [ frac{1}{Omega_a} nabla_a W_{ab}(h_a) + frac{1}{Omega_b} nabla_a W_{ab}(hb) big ]$

$latex frac{d}{dt}v_{a,mag,dis} = $

Ecuación de la energía (con h: “grad-h” term, AV: Artificial Viscosity term y C: Condutivity):

$latex frac{d}{dt}u_{a,MDH} = frac{d}{dt} (du_{a,h} + du_{a,AV} + du_{a,C})$

donde

$latex frac{d}{dt}u_{a,h} = frac{1}{Omega_a} frac{P_a}{rho_a^2} sum_b m_b v_{ab} cdot nabla_a W_{ab}(h_a) $

$latex frac{d}{dt}u_{a,AV} = $

$latex frac{d}{dt}u_{a,C} = $

Tags: ,

Integración en el plano complejo.

Definición (Camino, camino opuesto y caminos equivalentes): Un camino es una aplicación

$latex gamma:[a,b] longrightarrow mathbb{C}$

de clase $latex C^1$ a trozos. Si además $latex gamma(a) = gamma(b)$ entonces tenemos un camino cerrado.

Se llama camino opuesto de $latex gamma$ a la aplicación

$latex (-gamma):[-b,-a] longrightarrow mathbb{C} ,/, t mapsto (-gamma)(t):=gamma(-t)$

o, de manera equivalente, para mantener el mismo dominio

$latex (-gamma):[a,b] longrightarrow mathbb{C} ,/, t mapsto (-gamma)(t):=gamma(a+b-t)$.

Dados dos caminos $latex alpha:[a,b] longrightarrow mathbb{C}$ y $latex beta:[c,d] longrightarrow mathbb{C}$ son equivalentes si existe $latex varphi:[c,d] longrightarrow [a,b]$ biyectiva de clase $latex C^1$a trozos creciente de manera que $beta = alpha circ varphi$.

Definición (Integración a lo largo de un camino): Sea $latex f:U subset mathbb{C} longrightarrow mathbb{C}$ una función continua y sea $latex gamma:[a,b] longrightarrow mathbb{C}$ un camino $latex C^1$ tal que $latex gamma^* := gamma([a,b]) subset U$. Entonces:

$latex int_gamma f(z) dz = int_a^b f(gamma(t)) gamma'(t) dt$

Propiedades: Dada $latex f:A subset mathbb{C} longrightarrow mathbb{C}$ es sencillo demostrar

  1. $latex int_gamma a f(z)+b g(z) dz = a int_gamma f(z)dz + b int_gamma f(z)dz$.
  2. Si $latex alpha$ y $latex beta$ son equivalentes en $latex A$ entonces $latex int_alpha f(z)dz = int_beta f(z)dz$.
  3. $latex int_gamma f(z)dz = -int_{-gamma} f(z)dz$.
  4. Si $latex gamma = alpha cup beta$ entonces $latex int_gamma f(z)dz = int_alpha f(z)dz + int_beta f(z)dz$.
  5. $latex |int_gamma f(z)dz| leq M L(gamma)$ donde $latex M:= max { |f(z)|: z in gamma^* }$ y $latex L(gamma)$ es la longitud de la curva $latex gamma$.
  6. Teorema fundamental de cálculo: $latex int_gamma f(z) dz = F(gamma(b)) – F(gamma(a))$, siendo $latex F$ una primitiva de $latex f$. En particular, si $latex gamma$ es cerrado y $latex gamma(a) = gamma(b)$ entonces $latex int_gamma f(z) dz = 0$.

Teorema de Cauchy.

Teorema (de Cauchy)

Corolario (Existencia de primitiva de una función analítica)

Teorema (de Morera)

Fórmula de Cauchy.

Fórmula integral de Cauchy: Sea $latex A$ un abierto, sea $latex overline{D(z_0,R)} subset A$ y $latex f in mathcal{H}(A)$. Entonces

$latex f(z) = frac{1}{2 pi i} int_{C(z_o,R)} frac{f(u)}{u-z} du$ siempre que $latex |z-z_0|<R$

Fórmula integral de Cauchy para derivadas: Sea $latex f in mathcal{H}(A)$, $latex overline{D(z_0,R)} subset A$. Entonces

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

Teorema (de Liouville):

Teorema (fundamental del álgebra):

Teorema (de Weierstrass):

Definición (Indice de un punto respecto de un camino): Sea $latex gamma$ un camino y sea $latex A = mathbb{C} – gamma^*$. Entonces se define el índice de $latex z in A$ respecto de $latex gamma$ al número:

$latex ind_gamma(z) = frac{1}{2pi i} int_gamma frac{du}{u-z}$

Tags: , , , , , , ,

El Dr. Treb Lae y la Dra. Aini Vadh simulan colisiones de objetos compactos. Para poder hacerlas es para lo que me utilizan. Mi nombre es Tnala y soy una supercomputadora. Pero no una cualquiera, no, soy una de las que aparece en el TOP500.

—Vamos a hacer otra simulación —dice Treb—. ¿Qué colocamos esta vez?

—Probemos con una estrella de neutrones y un agujero negro. Pero vamos a cambiar algunas cosillas —contesta Aini.

Coge el ratón, se inclina sobre la pantalla y empieza a trabajar. Para empezar, distribuye los objetos sobre el escenario de la simulación, a continuación, define los parámetros que caracterizan a cada uno de los objetos compactos y, finalmente, modifica algunas cuestiones referentes a la propia simulación. Pero cuando ya parece que lo tiene, se reclina sobre la silla y se queda pensando. No lo tiene claro. Casi siempre, antes de empezar cualquier simulación, cambian una y otra vez las cosas. Es por eso que, hasta que no empiezo a simular, no pierdo el tiempo creado una y otra vez los diferentes objetos. Simplemente me guardo la información para luego.

—Creo que será mejor trabajar solo con estrellas —afirma finalmente Aini con determinación—. Probemos con dos estrellas con estas ecuaciones de estado, ¿te parece?

Treb se lo piensa. Tras el silencio esboza una leve sonrisa. Acaba de entender lo que pretende Aini.

—Perfecto. Creo que esta vez si que lo conseguimos —contesta feliz—. ¡Eres realmente brillante!

Aini sonríe contenta e inmediatamente lanza la simulación.

—¡Ya está! —exclama—. Vamonos. A ver si mañana tenemos buenas noticias.

Se marchan a casa a dencansar y me dejan. Yo no necesito descansar. Ahora es cuando me pongo manos a la obra.

Lo primero que hago es crear el escenario y colocar los objetos tal y como me lo han definido, para lo que utilizo la información que me he guardado previamente. A partir de ahora si que habrá cuerpos compactos creados de verdad. Inicializo el campo gravitatorio (una malla con nodos para la parte elíptica), añadiendo los posibles BH, e inicializó los fluidos de las NS, SS o DW (las partículas para la parte hiperbólica) junto con la descomposición espacial para las relaciones de partículas vecinas. Me guardo el estado inicial.

Ahora si que empiezo con la simulación, que es un vaivén entre la parte elíptica y la hiperbólica, es decir, miraré como el fluido me deforma el espacio-tiempo, desplazaré el fluido según esa deformación y volveré buscar como esa nueva distribución del fluido me modifica nuevamente el espacio-tiempo y vuelta a empezar.

Técnicamente, primero resuelvo el sistema elíptico para el paso de tiempo actual con lo que obtengo una métrica y una curvatura que necesitaré conocer para desplazar las partículas del fluido, que seguirán geodésicas. A continuación tengo que desplazar el fluido. Necesito conocer la densidad, la presión y la velocidad que recuperaré de la densidad, el momento y la energía. Para cada partícula calculo el valor de sus variables conservadas, para lo que necesito sus vecinas, el kernel, las eos y el solve rk4. Lo ideal sería tener un sistema que me resuelva esto para todas las partículas de forma simultánea (¿se puede?) y que en un futuro paralelizaré. Una vez resuelto, desplazo las partículas, actualizo sus vecinos y muevo los agujeros. Almaceno el estado. Finalmente, y antes de volver a empezar, revisó si el paso de tiempo es el apropiado (la simulación de la perdida de momento y de energía fruto de la radiación gravitacional en forma de ondas gravitatorias, así como la viscosidad y la conductividad térmica ya están incluidas en las ecuaciones discretizadas que usamos?).

Hago esto tantas veces como me han configurado, de manera que, en cuanto acabo, tengo almacenados el estado de la simulación en cada instante de tiempo, por lo que Treb i Aini podran visualizarlo y estudiarlo. Creo que el resultado si que es el que se esperaban. Me alegro por ellos (si es que puede una máquina alegrarse :-)).

La siguiente figura muestra un diagrama en el que se ve como se han comunicado los diferentes objetos:

A proposito, ¿sabeis a quién hacen honor los nombres de Treb Lae, Aini Vadh y Tnala? Envia un comentario con la respuesta correcta y participarás en el sorteo… 🙂

Tags: , ,

Definición (Función exponencial): Sea $latex z in mathbb{C}$, definimos:

$latex exp(z) = e^z := sum_{n geq 0} frac{z^n}{n!} in C^infty(mathbb{C})$.

Propiedades: Es sencillo verificar que $latex (e^z)’ = e^z$, $latex e^{z+w} = e^z e^w$ y $latex e^z neq 0, e^{-z} = frac{1}{e^z}$. Además:

$latex e^z = 1 Leftrightarrow z = 2pi n i$ con $latex n in mathbb{N}$.

Por lo que la función exponencial compleja es periódica de periodo imaginario $latex 2pi$.

Definición (Funciones trigonométricas e hiperbólicas): Sea $latex z in mathbb{C}$, definimos:

$latex sin(z):=frac{e^{iz}-e^{-iz}}{2i}$, por lo que $latex sin(z) = sum_{n geq 0} (-1)^n frac{z^{2n+1}}{(2n+1)!}$

$latex cos(z):=frac{e^{iz}+e^{-iz}}{2}$, por lo que $latex cos(z) = sum_{n geq 0} (-1)^n frac{z^{2n}}{(2n)!}$

$latex sinh(z):=frac{e^{z}-e^{-z}}{2}$ y $latex cosh(z):=frac{e^{z}+e^{-z}}{2}$

Propiedades: Las funciones trigonométricas y las hiperbólicas son holomorfas y de clase $latex C^infty(mathbb{C})$. Además, se cumple:

$latex sin(z pm w) = sin(z) cos(w) pm cos(z) sin(w)$,

$latex cos(z pm w) = cos(z) cos(w) mp sin(z) sin(w)$,

$latex sinh(z pm w) = sinh(z) cosh(w) pm cosh(z) sinh(w)$,

$latex cosh(z pm w) = cosh(z) cosh(w) pm sinh(z) sinh(w)$.

Además, $latex sin^2 z + cos^2 z = 1$, son $latex 2pi$-periódicas, $latex sin(z) = 0 Leftrightarrow z = n pi$ y $latex cos(z)=0 Leftrightarrow z = frac{pi}{2} + n pi$ con $latex n in mathbb{N}$.

Definición (Argumento): Sea $latex z in mathbb{C} – { 0}$, diremos que $latex alpha$ es un argumento de $latex z$ si $latex z = |z|(cos alpha + i sin alpha)$. Definimos:

$latex arg z := { alpha in mathbb{R}: z=|z|(cos alpha + i sin alpha)}$,

y si $latex alpha_0 in arg z$ entonces:

$latex arg z := { alpha_0 + 2 pi n : n in mathbb{Z} }$.

Definición (Logaritmo y ramas del logaritmo): Sea $latex z in mathbb{C}-{0}$. Diremos que $latex w in mathbb{C}$ es un logaritmo de $latex z$ si $latex e^w = z$. Así pues, definiremos:

$latex log z := { w in mathbb{C}: e^w = z}$.

Además, si $latex w=x+iy$ entonces $latex w = log z = ln |z| + i arg z$

Definición: Sea $latex alpha in mathbb{R}$ y $latex H_alpha={ -r(cos alpha + i sin alpha): r geq 0 }$. Definimos $latex arg_alpha: mathbb{C} – H_alpha longrightarrow ]alpha – pi, alpha + pi[$ como el único argumento de $latex z$ en el intervalo $latex ]alpha – pi, alpha + pi[$.

Definición: Sea $latex alpha in mathbb{R}$ y $latex H_alpha = { -r (cos alpha + i sin alpha): r geq 0 }$. Entonces:

$latex log_alpha: mathbb{C} – H_alpha longrightarrow mathbb{C} / z mapsto log_alpha(z) = ln |z| + i arg_alpha z$

Teorema: La función $latex log_alpha$ es derivable en $latex mathbb{C} – H_alpha$ y la derivada es $latex (log_alpha z)’ = frac {1}{z}$. Además, $latex log_0 (1+z) = sum_{n geq o} (-1)^n frac{z^{n+1}}{n+1}$ donde $latex log_0$ es el logaritmo principal ($latex alpha = 0$ y $latex H_0 = { -r: r geq 0 }$ ). El logaritmo no se puede extender continuamente.

Tags: , ,

El software DualSPHysics es la versión paralela con CPUs y GPUs de SPHysics. El código se ha implementado esta vez en C++ y CUDA.

A continuación mostramos dos diagramas de colaboración de DualSPHysics generados automáticamente por doxygen,

Tags: , , ,

Definición (Función holomorfa): Sean $latex U subset mathbb{C}$ abierto, $latex f: U longrightarrow mathbb {C}$ y $latex z_0 in U$. Decimos que $latex f$ es holomorfa (derivable o entera) en $latex z_0$ si existe:

$latex lim_{z rightarrow z_0} frac{f(z)-f(z_0)}{z-z_0} = f'(z_0)$.

Diremos que es derivable en $latex U$ si lo es todo sus puntos. Representaremos por $latex mathcal{H}(U)$ al conjunto de todas las funciones holomorfas en $latex U$. Si $latex f in mathcal{H}(U)$ entonces también es continua en $latex U$.

Teorema (Reglas de derivación): es sencillo deducir las reglas de derivación compleja para la suma, producto, cociente y la regla de la cadena.

Teorema: Sea $latex f(z) = sum_{n geq 0} c_n (z-z_0)^n$ una serie de potencias con radio de convergencia $latex R>0$. Entonces $latex f$ es holomorfa en $latex D(z_0,R)$ y $latex f'(z) = sum_{n geq 1} n c_n z^{n-1}$ con radio de convergencia $latex R$.

Corolario: Si $latex f$ es analítica en $latex U$ entonces $latex f$ es holomorfa en $latex U$ y $latex f in C^{infty}(U)$ (el recíproco veremos que también es cierto).

Definición (Ecuaciones de Cauchy-Riemann): Si $latex f(z) = u + iv$ es derivable en $latex z_0 = x_0 + i y_0$, entonces existen las derivadas parciales $latex u_x, u_y, v_x, v_y$ y verifican las ecuaciones de Cauchy-Riemann:

$latex u_x = v_y$ y $latex u_y = -v_x$.

De esta manera, las funciones derivables complejas no son mas que funciones diferenciables de $latex mathbb{R}^2$ cumpliendo unas ecuaciones adicionales.

Teorema (Condición suficiente de derivabilidad): Sea $latex f: U longrightarrow mathbb{C}$ continua con $latex U subset mathbb{C}$ abierto y $latex f(z)=u+iv$. Si las cuatro derivadas parciales $latex u_x, u_y, v_x, v_y$ existen, son continuas en $latex U$ y verifican las ecuaciones de Cauchy-Riemann, entonces $latex f$ es derivable en $latex U$.

Tags: ,

Mostramos como, comentando correctamente nuestro diseño, Visual Paradigm for UML nos genera reports en PDF o HTML con documentación sobre nuestra aplicación.

A continuación mostramos un pequeño fragmento de como queda el HTML:

Tags: , ,

Consideraremos el plano complejo $latex mathbb{C}$ dotado de la topología inducida a partir del módulo con la distancia $latex d(z,w) = |z-w|$. Llamaremos disco abierto de centro $latex z_0$ y radio $latex r$ al conjunto:

$latex D(z_0,r):= { z in mathbb{C}: |z-z_0|<r }$.

Dada una sucesión de números complejos $latex { z_n }_{n geq 0}$, diremos que es convergente a un número complejo $latex z$ cuando $latex |z-z_0| rightarrow 0$ si $latex n rightarrow infty$.

Una serie $latex sum_{n geq 0} z_n$ es convergente a $latex z in mathbb{C}$ cuando la sucesión de sumas parciales $latex { sum_{n=0}^m z_n}_{m geq 0}$ es convergente a $latex z$.

La serie $latex sum_{n geq 0} z_n$ es absolutamente convergente cuando $latex sum_{n geq 0} |z_n|$ converge.
Teorema (Criterio de la raíz o Cauchy, de Dirichlet, de sumación parcial de Abel, de Dirichlet y de Abel): son una serie de criterios que permiten determinar cuando una serie es convergente, absolutamente convergente o divergente.

Una serie de potencias centrada en $latex z_0 in mathbb{C}$, es una serie de la forma:

$latex sum_{n geq 0} c_n (z-z_0)^n$

con $latex c_n, z in mathbb{C}$. Al punto $latex z_0$ se le llama centro de la serie y a la sucesión $latex { c_n }$ coeficientes de la serie de potencias.

Definición (Radio de convergencia): Definimos el radio de convergencia como:

$latex R = (limsup_{n rightarrow infty} sqrt[n]{|c_n|})^{-1} in [0,+infty[$,

donde asumiremos que $latex frac{1}{+infty}=0$ y $latex frac{1}{0} = +infty$ y que $latex R = +infty Rightarrow D(z_0,R) = mathbb{C}$.

Teorema (Cauchy-Hadamard): Sea $latex R$ el radio de convergencia de $latex sum_{n geq 0} c_n (z-z_0)^n$, entonces la serie converge absolutamente en $latex |z-z_0|<R$ y diverge en $latex |z-z_0|>R$. Además, la serie converge uniformemente en $latex |z-z_0| leq r$ con $latex 0 leq r < R$ (subconjuntos compactos de $latex D(z_0,R)$).

Corolario: la serie de potencias $latex sum_{n geq 0} c_n (z-z_0)^n$ define una función continua en el interior del disco de convergencia $latex |z-z_0| < R$.

Teorema (Principio de los ceros aislados): Sea

$latex f(z)= sum_{n geq 0} c_n (z-z_0)^n$

con radio de convergencia $latex R > 0$. Entonces si $latex f(z_0) = 0$ y $latex f$ no es idénticamente nula, existe $latex delta > 0$ tal que $latex f(z) neq 0$ si $latex 0< |z-z_0|< delta$.

Lo que nos dice este principio es, sencillamente, que los ceros de las funciones analíticas no son puntos de acumulación, es decir, que son puntos aislados.

Definición (Función analítica): Sea $latex U subset mathbb{C}$ un abierto, $latex f:U longrightarrow mathbb{C}$ y $latex z_0 in U$. Diremos que $latex f$ es analítica en $latex z_0$ si existe $latex D(z_0,R) subset U$ tal que $latex f$ podemos escribirla mediante una serie de potencias centrada en $latex z_0$ en $latex D(z_0,R)$. Diremos que $latex f$ es analítica en $latex U$ si lo es en cada punto de $latex U$.

Ejemplo: los polinomios y las series de potencias centradas en $latex z_0$ son funciones analíticas en $latex z_0$.

Definición (Prolongación analítica): Si $latex h: V longrightarrow mathbb{C}$ y $latex f: U longrightarrow mathbb{C}$ son dos funciones analíticas con $latex V subset U subset mathbb{C}$ abiertos de manera que $latex f(z) = h(z)$ si $latex z in V$, o sea, $latex f|_V = h$, diremos que $latex f$ es una prolongación analítica de $latex h$.

Teorema (Principio de prolongación analítica): Sea $latex U$ un abierto conexo del plano. Sea $latex V subset U$ un conjunto de puntos que tiene puntos de acumulación en $latex U$, o sea, $latex ac_U(V) neq emptyset$. Entonces, si $latex f$ es analítica en $latex U$ y la restricción de $latex f$ sobre $latex V$ es $latex f|_V$ = 0, entonces $latex f=0$ en $latex U$.

Corolario: Sean $latex f, g$ analíticas en un abierto conexo $latex U$ y sea $latex V subset U$ tal que $latex ac_U(V) neq emptyset$ y $latex f|_V = g|_V$. Entonces $latex f = g$ en $latex U$.

Con todo esto, resulta que las extensiones analíticas son únicas si el dominio de éstas es conexo, es decir, si $latex f: U longrightarrow mathbb{C}$ y $latex g: U longrightarrow mathbb{C}$ con $latex U subset mathbb{C}$ conexo son dos prolongaciones analíticas de $latex h: V longrightarrow mathbb{C}$, entonces $latex f = g$ en todos los puntos. Esto se debe a que $latex f-g: U longrightarrow mathbb{C}$ se anula en $latex V$ que es un conjunto no vacío de $latex U$, y una función analítica que se anula en un conjunto no vacío debe anularse en todo su dominio si este es conexo, por lo que debe ser $latex f-g equiv 0$.

Tags: , , ,

Cuando exigimos que las funciones de transición del atlas que recubre una variedad diferenciable sean holomorfas podemos decir que tenemos una variedad compleja. En particular, toda variedad compleja de dimensión $latex n$ será una variedad diferenciable de dimensión $latex 2n$ dotada de una orientación natural. Las superfícies de Riemann, o los grupos de Lie con operaciones de grupo holomorfas son ejemplos de variedades complejas.

Las variedades lorentzianas son importantes en relatividad general. La mecánica cuántica tiene su base matemática en los espacios de Hilbert complejos separables ($latex L^2(mathbb{R})$ o $latex mathbb{C}^{2n+1}$).

Existen diferentes maneras de definir el cuerpo de los números complejos. La mas intuitiva es hacer tomar el conjunto $latex mathbb{R} times mathbb{R}$ y definir las operaciones internas:

  • $latex (x_1,y_1) + (x_2,y_2) = (x_1+x_2,y_1+y_2)$ .
  • $latex (x_1,y_1) cdot (x_2,y_2) = (x_1 x_2 – y_1 y_2, x_1 y_2 + x_2 y_1)$.

con $latex (x_1,y_1), (x_2, y_2) in mathbb{R} times mathbb{R}$. Es sencillo comprobar que $latex mathbb{C} := (mathbb{R} times mathbb{R}, +, cdot)$ tiene estructura de cuerpo.

Otra manera mas elegante es ver que $latex mathbb{C} cong frac{mathbb{R}[x]}{(x^2+1)}$ donde $latex (x^2+1)$ es el ideal generado por el polinomio irreducible $latex x^2+1 in mathbb{R}[x]$, pues existe una propiedad que nos dice que un anillo conciente de este tipo tiene estructura de cuerpo.

El análisis complejo se encarga del estudio de las funciones de variable compleja. En próximos posts exponemos de manera breve algunos de los contenidos básicos del área: funciones analíticas, funciones holomorfas, funciones elementales, integración en el plano complejo con formula y teorema de Cauchy y teoremas de Laurent y de los residuos con algunas aplicaciones.

Tags: , , , , , , ,

Nada que se mueva en el tejido del espacio-tiempo puede superar la velocidad de la luz. Sin embargo, el propio tejido no esta sometido a esta restricción (o al menos eso sugiere la teoría de la inflación cósmica para resolver el problema del horizonte).

La materia provoca deformación en el espacio-tiempo. Lo que propone M. Alcubierre en su artículo “The Warp Drive: Hyper-fast Transluminic within General Relativity”, es crear una distorsión local del espacio-tiempo de manera que produzca una expansión detrás de la nave espacial, y una contracción opuesta por delante de ella. De esta manera, es el propio espacio-tiempo el que empuja lejos de la Tierra a la nave y la atrae hacia la estrella distante a la que se pretende viajar.

En el formalismo $latex 3+1$, que permite una clara interpretación de los resultados, la métrica se escribe.

$latex ds^2 = -dtau^2 = g_{alphabeta}dx^alpha dy^beta = $

$latex = -(alpha^2 – beta_i beta^i) dt^2 + 2 beta_i dx^i dt + gamma_{ij} dx^i dx^j$

Una manera de conseguir lo pretendido es hacer:

$latex alpha = 1, beta^x = -v_s(t)f(r_s(t)), beta^y = beta^z = 0, gamma_{ij}=delta_{ij}$

con

$latex v_s(t)=frac{dx_s(t)}{dt}$ y $latex r_s(t)=[(x-x_s(t))^2+y^2+z^2]^{1/2}$

y

$latex f(r_s) = frac{tanh(sigma(r_s+R))-tanh(sigma(r_s-R))}{2 tanh (sigma R)}$

con $latex R>0$ y $latex sigma >0$ arbitrarios.

La siguiente imagen del artículo original de Alcubierre muestra como queda deformado el espacio-tiempo con esta métrica (corresponde a los valores $latex sigma=8$ y $latex R = v_s =1$ en la métrica) consiguiendo la burbuja warp donde se sitúa la nave:

El Prof. Dr. Daniel Weiskopf que trabaja, entre otras cosas, en Special and general relativistic visualization, ya editó en el año 2000 un pequeño video en el que simula, mediante técnicas de Ray tracing, como veriamos diferentes cuerpos del Sistema Solar si por el espacio-tiempo que hay entre ellos y nosotros pasara una nave viajando dentro de una burbuja warp a diferentes velocidades, sub y superlumínicas, cercanas a la de la luz:

[youtube http://www.youtube.com/watch?v=aFyHOUCfplc?rel=0&w=420&h=315]

Dentro de la página que hemos enlazado anteriormente hay otros vídeos mas actuales.

Tags: , , , , , , , ,

Otro diagrama disponible en UML 2 es el diagrama de comunicación.

El diagrama de comunicación permite modelar la interacción entre los diferentes objetos que se produce mediante mensajes en secuencia, es decir, muestra que mensajes se pasan los objetos entre si y en que orden. Es un diagrama muy útil, pues muestra tanto información estática, tomada del diagrama de clases, como información dinámica, tomada del diagrama de casos de uso y del  diagrama de secuencia.

A continuación muestro el diagrama de comunicación que he realizado mediante Visual Paradigm for UML de SPHysics:

Técnicamente, tendríamos que hablar de un pseudodiagrama de comunicación, pués el lenguaje de programción es fortran que no trabaja con orientación a objetos, por lo que en lugar de clases  tenemos módulos. Sin embargo, y salvando las distancias, cumple sobradamente su propósito semántico.

Una vez un profesor nos dijo algo así como: “No programeis lo que ya esta programado”, una especie de versión light de la reinvención de la rueda. Lo que quería decir era, sencillamente, que antes de empezar a diseñar algo desde cero para resolver un problema teniamos que investigar si ya existia un software capaz de hacerlo. Si la respuesta era negativa, adelante, pero si la respuesta era positiva, teniamos que plantearnos una serie de preguntas: ¿Es gratuito o hay que pagar por él? ¿Se adapta perfectamente a mis necesidades o necesita cierto grado de adaptación? ¿Tengo acceso al código fuente y tengo permiso para modificarlo a mi medida? ¿Que piensa el cliente al respecto? etc.

En el caso que nos ocupa, un sofware que promete es SPHysics, un esfuerzo conjunto de investigadores en la Johns Hopkins University (U.S.A.), la University of Vigo (Spain), la University of Manchester (U.K.) y la University of Rome La Sapienza (Italy). Evidentemente, se trata de hidrodinámica newtoniana y no magnetohidrodinámica relativista, por lo que hay un salto cualitativo en la física considerada, además de que se trata de programación modular y no orientada a objetos. Sin embargo, al ser de código abierto y disponer de versión paralela sobre CUDA, nos puede servir de mucha ayuda a la hora de implementar muchas cosas.

A continuación se enlazan dos simulaciones con 100M y 226M partículas en las que se puede apreciar la potencia del software:

[youtube http://www.youtube.com/watch?v=o6_iPdIQmyc?rel=0&w=420&h=315]

[youtube http://www.youtube.com/watch?v=rHq6rcdqHzM?rel=0&w=420&h=315]

Tags: , , , , , , ,

En el artículo “Introducción a la relatividad numérica” de M. Alcubierre, también explica el formalismo $latex 3+1$.

Resolver las ecuaciones de campo de Einstein en la práctica es complicado, ya que son un sistema de $latex 10$ EDPs en $latex 4D$ acopladas y no lineales con muchísimos términos. Se conocen soluciones exactas en situaciones muy concretas con un alto grado de simetría espacial o temporal (simetría esférica, simetría axial, soluciones estáticas, homogéneas, isotrópicas, etc.), pero en la mayoría de situaciones interesantes en astrofísica no se dan estas condiciones y tenemos que resolverlas utilizando aproximaciones numéricas mediante complicados programas.

En relatividad numérica, la idea es separar las ecuaciones de campo de Einstein de forma que podamos dar ciertas condiciones iniciales y, a partir de ellas, obtener la evolución del campo gravitacional. Existen diferentes maneras de hacerlo y el formalismo $latex 3+1$, que es el mas ampliamente utilizado, lo hace separando las tres componentes espaciales de la temporal.

Para estudiar la evolución en el tiempo, lo primero que se hace es formular un problema de Cauchy. En las ecuaciones de Einstein, el espacio y el tiempo son simétricos, por lo que primero debemos separar uno de otro. Lo que queremos es un espacio-tiempo orientable temporalmente de manera que podamos elegir de manera contínua a través del espacio-tiempo la mitad del cono de luz que contituye la dirección futura y de la mitad que corresponde a la dirección pasada. A la formulación de la relatividad general (GR) resultante de esta separación es el formalismo $latex 3+1$.

Un conjunto abierto $latex U$ de un espacio-tiempo es globalmente hiperbólico sii:

  1. Para cualquier par de puntos $latex p$ y $latex q$, el conjunto $latex gamma^+(p) cap gamma^-(q) subset U$ y es compacto, donde $latex gamma^pm(S)$ son el futuro y pasado causal de una region $latex S$.
  2. No existen curvas espacio-temporales cerradas que pasen por $latex U$, prohibiendo los viajes al pasado en el tiempo y cumpliendose de esta manera el pricipio de causalidad en el abierto.

Una $latex p$-foliación de una variedad $latex M$ de dimension $latex n$ consiste en una partición de ésta en subvariedades diferenciables $latex {N_i}_{i in I}$ de dimensión $latex dim N_i = p<;m,,forall iin I$, por lo que localmente $latex M$ tiene estrutura topologica de variedad producto. Por ejemplo, una $latex 2$-foliación del espacio $latex mathbb{R}^3$ es $latex { mathbb{R}^2_z}_{z in mathbb{R}}$. La foliación de una variedad no es única (por ejemplo, $latex { mathbb{R}^2_y}_{y in mathbb{R}}$ es otra foliación posible de $latex mathbb{R}^3$, asumiendo que sabemos de que hablamos cuando nos referimos a $latex z$ o $latex y$: planos con vector normal $latex (0,0,1)$ en el primer caso y $latex (0,1,0)$ en el segundo).

Todo espacio-tiempo globalmente hiperbólico puede ser foliado, es decir, separado en cortes tridimensionales apilados, de tal forma que cada hoja es una hipersuperficie de Cauchy espacial. Sea $latex (t,x^i)$ un sistema de coordenadas tal que la función $latex t$ es de gradiente temporal. Entonces las superficies $latex t=cte$ (un tiempo “universal” que no tiene porque ser el tiempo propio de nadie) definen una foliación del espacio-tiempo. Denotaremos por $latex Sigma_t$ a la hipersuperfície de Cauchy de tiempo $latex t$.

Dada una foliación $latex Sigma_t$ definida para la función $latex t$, podemos encontrar campos vectoriales $latex xi$ de manera que $latex mathcal{L}_xi t = 1$. Llamamos base de evolución a la pareja $latex (xi,t)$. Podemos descomponer $latex xi$ de manera relativa a un observador euleriano:

$latex xi = alpha n + beta$

donde $latex n = frac{dt}{|dt|}$ con $latex g(n,n)=-1$ es la normal a la foliación, $latex alpha$ es la función de paso y $latex beta$ el vector desplazamiento de la base de evolución. A cada base de evolución podemos asociarle unas coordenadas adaptadas de manera que $latex xi = partial_t$ y $latex (x^i)$ son coordenadas en $latex Sigma_t$, de manera que $latex beta = beta^i partial_i$.

Considerar diferentes $latex alpha$ es considerar diferentes foliaciones, por lo que el paso entre $latex Sigma_0$ y $latex Sigma_t$ depende del punto de $latex Sigma_0$ considerado: $latex alpha = alpha(t,x^i)$.

El vector desplazamiento $latex beta$ determina $latex xi$, define el difeomorfismo entre $latex Sigma_0$ y $latex Sigma_t$: si $latex beta=0$ entonces $latex varphi_t(P_0) = bar{P_t}$ con $latex P_0 in Sigma_0$ y $latex bar{P_t}in Sigma_t$ ambos sobre la misma curva integral de $latex n$, por lo que tenemos evolución sin desplazamiento. Si $latex beta neq 0$ entonces $latex varphi_t(P_0) = P_t$ en $latex Sigma_t$ desplazado respecto $latex bar{P_t}$.

Por lo tanto, tenemos:

$latex g_{munu} = begin{pmatrix} -alpha^2 + beta_k beta^k & beta_i \ beta_j & gamma_{ij} end{pmatrix}$

y

$latex n^mu = frac{1}{alpha}(1,-beta^i)$, $latex n_mu = (-alpha,0)$

Tags: , , , , ,

Ya hemos comentado que podemos ver un campo tensorial diferenciable como una generalización de funciones, campos vectoriales y $latex 1$-formas. Estudiaremos ahora el caso de las métricas.

Como comentamos en un post anterior, una métrica de Riemann:

$latex g_m: T_mM times T_mM longrightarrow mathbb{R}$

podemos verla como un campo tensorial dos veces covariante, de tipo $latex (0,2)$. Efectivamente, ya que en cada $latex m in M$ tenemos definido:

$latex g_m in mathcal{L}(T_mM times T_mM, mathbb{R}) cong otimes^2 T_m^*M = T_m^{(0,2)}M$.

Por lo tanto, tenemos que $latex g: M longrightarrow T^{(0,2)}M$ define una métrica sobre la variedad $latex M$.

Por ejemplo, la métrica de Schwarzschild en coordenadas de Schwarzschild $latex (r,theta,varphi,tau)$ es:

$latex ds^2 = frac{1}{1-frac{2M}{r}}dr^2+ r^2 (dtheta^2 + sin^2 theta dvarphi^2)-(1-frac{2M}{r})dtau^2$

que en notación de productos tensoriales queda:

$latex g = frac{1}{1-frac{2M}{r}}dr otimes dr + r^2 (dtheta otimes dtheta + sin^2 theta dvarphi otimes dvarphi)-(1-frac{2M}{r})dtau otimes dtau$

Acabamos de ver que $latex g_m in mathcal{L}(T_mM times T_mM, mathbb{R})$. En nuestro caso, una base de $latex T_mM$ es:

$latex { frac{partial}{partial r}|_m, frac{partial}{partial theta}|_m, frac{partial}{partial varphi}|_m, frac{partial}{partial tau}|_m}$

por lo que $latex dim T_mM = 4$ y su base dual:

$latex { dr_m, dtheta_m, dvarphi_m, dtau_m }$

es una base de $latex T_m^*M$ con $latex dim T_m^*M = 4$. Como:

$latex mathcal{L}(T_mM times T_mM, mathbb{R}) cong otimes^2 T_m^*M = T_m^{(0,2)}M$

tenemos que:

$latex dim T_m^*M otimes T_m^*M = dim T_m^*M cdot dim T_m^*M = 4 cdot 4 = 16$

y una base de $latex T_m^*M otimes T_m^*M$ es:

$latex { dr_m otimes dr_m, dr_m otimes dtheta_m, dr_m otimes dvarphi_m, dr_m otimes dtau_m,$

$latex dtheta_m otimes dr_m, dtheta_m otimes dtheta_m, dtheta_m otimes dvarphi_m, dtheta_m otimes dtau_m,$

$latex dvarphi_m otimes dr_m, dvarphi_m otimes dtheta_m, dvarphi_m otimes dvarphi_m, dvarphi_m otimes dtau_m,$

$latex dtau_m otimes dr_m, dtau_m otimes dtheta_m, dtau_m otimes dvarphi_m, dtau_m otimes dtau_m }$

Las componentes de nuestra métrica en esta base son:

$latex g_{11} = frac{1}{1-frac{2M}{r}}$, $latex g_{22} = r^2$, $latex g_{33} = r^2 sin^2 theta$, $latex g_{44} = -(1-frac{2M}{r})$ y $latex g_{ij} = 0$ si $latex i neq j$.

En general, dada una variedad $latex M$ de dimensión $latex dim M = n$ y un punto $latex m in M$, entonces $latex dim T_mM = dim T_m^*M = n$ y $latex dim T_m^*M otimes T_m^*M = n^2$, por lo que $latex g_m in mathcal{M}_n(mathbb{R})$.

Si llamamos $latex x^1$ a la coordenada $latex r$, $latex x^2$ a la coordenada $latex theta$, $latex x^3$ a la coordenada $latex varphi$ y $latex x^4$ a la coordenada $latex tau$, entonces podemos referirnos a la métrica $latex g$ de una forma mas compacta:

$latex g = sum_{alpha,beta=1}^4 g_{alphabeta}dx^alpha otimes dx^beta$

que en física y siguiendo el convenio de suma de Einstein, con índices griegos variando de $latex 1$ a $latex 4$ y indices latinos haciendolo entre $latex 1$ y $latex 3$, queda:

$latex g_{alpha beta}dx^alpha dx^beta$

Si calculamos la inversa de la matriz $latex g_{alpha gamma}g^{gamma beta} = delta_alpha^beta$ obtenemos las componentes contravariantes de la métrica:

$latex g = sum_{alpha,beta=1}^4 g^{alpha beta} frac{partial}{partial x^alpha} otimes frac{partial}{partial x^beta} equiv g^{alpha beta} frac{partial}{partial x^alpha} frac{partial}{partial beta}$

que en el caso que nos ocupa son:

$latex g^{11} = 1-frac{2M}{r} $, $latex g^{22}= frac{1}{r^2}$, $latex g^{33}=frac{csc theta}{r^2}$ y $latex g^{44}=-frac{1}{1-frac{2M}{r}}$

por lo que nos queda:

$latex g= 1-frac{2M}{r} frac{partial}{partial r} otimes frac{partial}{partial r} + frac{1}{r^2}frac{partial}{partial theta} otimes frac{partial}{partial theta} + frac{csc theta}{r^2}frac{partial}{partial varphi} otimes frac{partial}{partial varphi} -frac{1}{1-frac{2M}{r}}frac{partial}{partial tau} otimes frac{partial}{partial tau}$

que también puede escribirse:

$latex g= 1-frac{2M}{r} partial_r^2 + frac{1}{r^2} partial_theta^2 + frac{csc theta}{r^2} partial_varphi^2 -frac{1}{1-frac{2M}{r}} partial_tau^2$

Tags: , , , , , , ,

Hemos hablado mucho de las ecuaciones de campo de Einstein pero aún no han aparecido de manera explícita. En el artículo “Introducción a la relatividad numérica” de M. Alcubierre, éste habla sobre ellas.

Las ecuaciones de campo de Einstein, derivadas buscando una generalización relativista y consistente de la ley de gravitación de Newton, como lo hizo Einstein, o de manera formal a partir de un principio variacional partiendo de un Lagrangiano adecuado, como lo hizo Hilbert, se escriben en su forma mas compacta como (signatura $latex (-,+,+,+)$ y $latex G=c=1$):

$latex G_{munu} = 8 pi T_{munu}$

donde $latex G_{munu}$ es el tensor de curvatura de Einstein que representa la geometría del espacio-tiempo, $latex 8 pi$ es un factor de normalización para obtener el límite Newtoniano correcto y $latex T_{munu}$ es el tensor de energia-momento que representa la distribución de materia y energía. Como $latex G_{mu nu}, T_{mu nu} in mathcal{M}_{16}(mathbb{R})$, tenemos $latex 16$ ecuaciones que se reducen a $latex 10$ por ser simétricos los dos tensores en sus dos índices. Son $latex 10$ PDEs acopladas en $latex 4D$.

El tensor de Einstein se define como:

$latex G_{mu nu} := R_{mu nu} – frac{1}{2} g_{mu nu}R$

donde $latex R_{mu nu}:=R^lambda_{mu lambda nu}$ es el tensor de Ricci ($latex R_{mu nu} in mathcal{M}_{16}(mathbb{R})$) que se obtiene contrayendo dos índices libres del tensor de curvatura de Riemann y $latex R:=g^{mu nu}R_{mu nu}$ es la traza del tensor de Ricci o la curvatura escalar.

El tensor curvatura está definido para toda variedad dotada de una conexión $latex nabla$:

$latex R(u,v)w = nabla_u nabla_v w – nabla_v nabla_u w – nabla_{[u,v]}w$

y nos permite hablar de transporte paralelo, nos dice el cambio que sufre un vector al transportalo paralelamente. En una variedad de Riemann siempre podemos definir una conexión libre de torsión, la conexión de Levi-Civita, que expresada en componentes queda:

$latex R^{rho}_{sigma mu nu} = partial_mu Gamma^rho_{sigma nu} – partial_nu Gamma^rho_{sigma mu} + Gamma^alpha_{sigma nu} Gamma^rho_{alpha mu} – Gamma^alpha_{sigma mu} Gamma^rho_{alpha nu}$

y que con $latex 4$ índices en $latex n$ dimensiones tiene $latex n^4$ componentes, de las que solo $latex 20$ (si $latex n=4$ y $latex 4^4=256$), al tener en cuenta simetrías, son independientes. Se puede demostrar que $latex R=0 Leftrightarrow$ variedad plana.

Recordar que se pueden subir y bajar índices contrayendo con el tensor métrico o su inverso:

$latex v_alpha = g_{alpha beta} v^beta$

$latex v^alpha = g^{alpha beta} v_beta$

$latex R_{rho sigma mu nu} = g_{rho alpha} R^{alpha}_{sigma mu nu}$

En el último ejemplo obtenemos la versión de la curvatura de Riemann totalmente covariante, un tensor de tipo $latex (0,4)$ (los elementos de la base pasan de ser de la forma $latex frac{partial}{partial_{x^alpha}} otimes dx^beta otimes dx^gamma otimes dx^delta$ de un tensor de tipo $latex (1,3)$ a ser de la forma $latex dx^alpha otimes dx^beta otimes dx^gamma otimes dx^delta$).

El tensor de curvatura de Riemann tiene las siguientes propiedades:

  1. Antisimetrías: $latex R_{alpha beta gamma delta} = – R_{alpha beta delta gamma} = -R_{beta alpha gamma delta}$.
  2. Simetrías: $latex R_{alpha beta gamma delta} = R_{gamma delta alpha beta}$.
  3. Primera identidad de Bianchi: $latex R_{alpha[betagammadelta]} = R_{alphabetagammadelta} + R_{alphagammadeltabeta} + R_{alphadeltabetagamma} = 0$.
  4. Segunda identidad de Bianchi:$latex R_{alphabeta[gammadelta;epsilon]} = R_{alphabetagammadelta;epsilon} + R_{alphabetadeltaepsilon;gamma} + R_{alphabetaepsilongamma;delta} = 0$.

El tensor de energia-momento describe la densidad de energia, la densidad de momento y el flujo de momento de un campo de materia:

$latex T^{00}$ = densidad de energía.

$latex T^{0i} = $ densidad de momento.

$latex T^{ij} = $ flujo de momento $latex i$ a través de la superficie $latex j$.

Las identidades de Bianchi son muy importantes porque nos llevan a:

$latex G^{mu nu},_{;nu} = 0 Rightarrow T^{mu nu},_{;nu} = 0$

que son las cuatro ecuaciones que representan la conservación local de la energía y del momento (la perdida de energía y momento en una región se compensa con el flujo de energía y momento fuera de esa región) donde $latex ;$ indica la derivada covariante.

Tags: , , , , , , , , , , , , ,

Sean $latex V_1, cdots, V_r$ espacios vectoriales de dimensión finita sobre $latex mathbb{R}$ y sean $latex V_1^*, cdots, V_r^*$ sus espacios duales.

Definimos el producto tensorial como el espacio vectorial de aplicaciones multilineales de $latex V_1^* times ldots times V_r^*$ en $latex mathbb{R}$, es decir:

$latex V_1 otimes ldots otimes V_r := mathcal{L}(V_1^* times ldots times V_r^*, mathbb{R})$

Si $latex v_1 in V_1, ldots , v_r in V_r$ y $latex sigma_1 in V_1^*, ldots, sigma_r in V_r^*$, entonces definimos $latex v_1 otimes , ldots , otimes v_r in V_1 otimes , ldots , otimes V_r$ como:

$latex v_1 otimes ldots otimes v_r (sigma_1, ldots, sigma_r)= sigma_1(v_1) ldots sigma_r(v_r)$

Si $latex dim V_j = n_j$ y sea $latex { e_i^j}_{i=1}^{n_j}$ una base de $latex V_j$ con $latex j=1,ldots,r$, entonces:

$latex {e_{i_1}^1 otimes ldots otimes e_{i_r}^r }_{1 leq i_j leq n_j, 1 leq j leq r }$

es una base de $latex V_1 otimes ldots otimes V_r$, de manera que $latex dim V_1 otimes ldots otimes V_r = n_1ldots n_r$.

Sea $latex V$ un espacio vectorial de dimensión $latex dim V = n$ y $latex V^*$ su dual. Construimos el espacio vectorial

$latex V^{(r,s)}:=(otimes^r V) otimes (otimes^s V^*)$

donde $latex otimes^k E:= E otimes overset{k)}{ldots} otimes E$ es la $latex k$-ésima potencia tensorial de $latex E$. A los elementos de $latex V^{(r,s)}$ se les llama tensores $latex r$ veces contravariantes y $latex s$ veces covariantes sobre $latex V$. Si $latex { e_1, ldots, e_n}$ es una base de $latex V$ y $latex { e^1, ldots, e^n}$ su base dual (los elementos $latex e^i$ son $latex 1$-formes: $latex e^i: V longrightarrow mathbb{R} in V^*$), entonces todo elemento de $latex V^{(r,s)}$ lo podemos escribir como:

$latex t = t^{i_1,ldots, i_r}_{j_1,ldots, j_s}e_{i_1} otimes ldots otimes e_{i_r} otimes e^{j_1} otimes ldots otimes e^{j_s}$

No es dificil demostrar $latex mathcal{L}(V,V) cong V otimes V^*$, $latex mathcal{L}(V times V, mathbb{R}) cong V^* otimes V^*$ y, en general:

$latex mathcal{L}(V times overset{k)}{ldots} times V, V) cong V otimes (otimes^k V^*)$.

Sea $latex M$ una variedad diferenciable y $latex m in M$. Entonces:

$latex T_m^{(r,s)} = (otimes^r T_mM) otimes (otimes^s T_m^*M)$

es un tensor $latex r$ veces contravariante y $latex s$ veces covariante de $latex M$ en $latex m$ y

$latex T^{(r,s)}M = bigsqcup_{m in M} T_m^{(r,s)}M$

es la variedad de tensores de tipo $latex (r,s)$ de $latex M$. Denotamos por $latex pi : T^{(r,s)}M longrightarrow M$ a la proyección que a cada tensor en $latex m$ le hace corresponder el punto $latex m$.

Un campo tensorial $latex r$ veces contravariante y $latex s$ veces covariante en $latex M$, de tipo $latex (r,s)$, es una aplicación diferenciable $latex K : M longrightarrow T^{(r,s)}M$ tal que $latex pi circ K = id$, es decir, que para cada $latex min M$ tenemos que $latex K_m := K(m) in T_m^{(r,s)}M$ (tenemos un campo tensorial definido en cada punto de la variedad).

Si $latex (U, varphi)$ es una carta, entonces:

$latex K|_U = K^{i_1,ldots,i_r}_{j_1,ldots,j_s} frac{partial}{partial varphi^{i_1}} ldots otimes frac{partial}{partial varphi^{i_r}} otimes dvarphi^{j_1} otimes ldots otimes dvarphi^{j_s}$

Los campos tensoriales son una generalización de:

  1. funciones: una función diferenciable $latex h:M longrightarrow mathbb{R}$ determina un campo tensorial de tipo $latex (0,0)$.
  2. campos vectoriales: un campo vectorial $latex X: M longrightarrow TM$ es un campo tensorial de tipo $latex (1,0)$, pues $latex TM = T^{(1,0)}$.
  3. $latex 1$-formas: una $latex 1$-forma $latex w: M longrightarrow T^*M$ es un campo tensorial de tipo $latex (0,1)$, ya que $latex T^*M = T^{(0,1)}$.

Tags: , ,

Dada una variedad diferenciable $latex M$, una curva diferenciable en $latex M$ es una aplicación diferenciable $latex alpha:]a,b[ longrightarrow M$ donde $latex a,b in mathbb{R}$ y $latex a<b$.

Dado un punto $latex m$:

  1. denotaremos por $latex mathcal{C}(M,m)$ al conjunto de curvas diferenciables $latex alpha:]a,b[ longrightarrow M$ tales que $latex a<0<b$ y $latex alpha(0) = m$.
  2. denotamos por $latex mathcal{F}(M,m)$ al conjunto de funciones diferenciables $latex f:U longrightarrow mathbb{R}$ tales que $latex U$ es un entorno abierto de $latex m$ en $latex M$.

Definimos en $latex mathcal{C}(M,m)$ la relación de equivalencia $latex sim$ de la siguiente manera:

$latex alpha sim beta Leftrightarrow (f circ alpha)'(0) = (f circ beta)'(0),,forall alpha,beta in mathcal{C}(M,m),,forall f in mathcal{F}(M,m)$

Llamamos espacio tangente a $latex M$ en $latex m$ a

$latex T_mM := frac{mathcal{C}(M,m)}{sim}$

A los elementos de $latex T_mM$ se les llama vectores tangentes a $latex M$ en $latex m$, y son clases de equivalencia de curvas diferenciables respecto de la relación anterior.

Denotamos por $latex p:mathcal{C}(M,m) longrightarrow T_mM$ a la aplicación que a cada curva $latex alpha$ le asocia su clase de equivalencia. Además, $latex v(f) :=(f circ alpha)'(0)$ donde $latex f in mathcal{F}(M,m)$ y $latex v = p(alpha)$ con $latex alpha in mathcal{C}(M,m)$.

Si la dimensión de la variedad es $latex n$ y $latex (U,varphi)$ es una carta cualquiera con $latex m in U$, entonces:

  1. denotamos por $latex { e_i}_{i=1}^{n}$ a la base canónica de $latex mathbb{R}^n$.
  2. definimos las curvas coordenadas $latex tau_i in mathcal{C}(M,m)$ como $latex tau_i(t) := varphi^{-1}(varphi(m)+te_i)$.
  3. definimos los vectores tangentes coordenados respecto de la carta $latex (U,varphi)$ como $latex frac{partial}{partial varphi^i}|_m = p(tau_i) in T_mM$.

Los vectores

$latex { frac{partial}{partial varphi^1}|_m, cdots, frac{partial}{partial varphi^n}|_m }$

forman una base de $latex T_mM$, por lo que para todo $latex v in T_mM$ tenemos que:

$latex v = sum_{i=1}^n v(varphi^i) frac{partial}{partial varphi^i}|_m$

y $latex dim T_mM = dim M = n$.

Llamamos espacio cotangente a $latex M$ en $latex m$ al espacio vectorial dual de $latex T_mM$ y que denotamos por $latex T_m^*M$, o sea:

$latex T_m^*M := (T_mM)^* = mathcal{L}(T_mM,mathbb{R})$

A los elementos de $latex T_m^*M$ se les denomina 1-formas de $latex M$ en $latex m$. Las 1-formas:

$latex { dvarphi_m^1, cdots, dvarphi_m^n }$

forman una base de $latex T_m^*M$ que es, a su vez, la base dual de:

$latex { frac{partial}{partial varphi^1}|_m, cdots, frac{partial}{partial varphi^n}|_m }$

Es evidente, por tanto, que $latex dim T_m^*M = dim T_m M = n$.

Llamamos fibrado tangente $latex TM$ a la variedad diferenciable que resulta de la unión disjunta de todos los espacios tangentes $latex T_mM$ con $latex m in M$, es decir:

$latex TM := bigsqcup_{m in M} T_mM$

Podemos construir un atlas para el fibrado a partir de los atlas de cada una de las variedades tangentes. A partir de esta construcción deducimos que $latex dim TM = 2n$

De la misma manera, llamamos fibrado cotangente $latex T^*M$ a la variedad diferenciable que resulta de la unión disjunta de todos los espacios cotangentes $latex T_m^*M$ con $latex m in M$, es decir:

$latex T^*M := bigsqcup_{m in M} T_m^*M$

Ademas, $latex dim T^*M = dim TM = 2n$

Tags: , , , , ,

Como ya comentamos, existen diferentes soluciones analíticas de las ecuaciones de Einstein correspondientes a los diferentes tipos de BH en equilíbrio. En la tesis “Evolution formalism of Einstein equations: numerical and geometrical issues” de I. Cordero podemos encontrarlas.

Para empezar, la consideración de variedades de Lorentz con simetría esférica y tensor de Ricci nulo, nos lleva a la métrica de Schwarzschild ($latex J=0, Q=0$) que podemos escribir en diferentes sistemas de coordenadas:

Coordenadas de Schwarzschild $latex (r,theta,varphi,tau)$ con $latex r > 2M$ y siendo $latex tau$ el tiempo propio:

$latex g = frac{1}{1-frac{2M}{r}}dr otimes dr + r^2 (dtheta otimes dtheta + sin^2 theta dvarphi otimes dvarphi)-(1-frac{2M}{r})dtau otimes dtau$

que en forma matricial queda:

$latex g=begin{bmatrix} frac{1}{1-frac{2M}{r}} & 0 & 0 & 0 \ 0 & r^2 & 0 & 0 \ 0 & 0 & r^2 sin^2 theta & 0 \ 0 & 0 & 0 & -(1-frac{2M}{r}) end{bmatrix}$

y en física se suele escribir:

$latex ds^2 = frac{1}{1-frac{2M}{r}}dr^2+ r^2 (dtheta^2 + sin^2 theta dvarphi^2)-(1-frac{2M}{r})dtau^2$

Además, como $latex dOmega^2 = dtheta^2 + sin^2theta dvarphi^2$ es la métrica de $latex S^2$ ($latex S^2(theta,varphi) = (sin theta cos varphi, sin theta sin varphi, cos theta)$ en $latex ]0,pi[ times ]0,2pi[$ de manera que $latex g_{11} = S^2_theta cdot S^2_theta = 1$, $latex g_{12} = g_{21} = S^2_theta cdot S^2_varphi = 0$ y $latex g_{22} = S^2_varphi cdot S^2_varphi = sin^2 theta$, con lo que $latex g = dtheta otimes dtheta + sin^2 theta dvarphi otimes dvarphi$), tenemos:

$latex ds^2 = frac{1}{1-frac{2M}{r}}dr^2+ r^2 dOmega^2 – (1-frac{2M}{r})dtau^2$

Coordenadas isotrópicas $latex (bar{r},theta,varphi,tau)$ con $latex r = bar{r} (1 + frac{M}{2bar{r}})^2$ respecto de las de Schwarzschild:

$latex ds^2 = (1+frac{M}{2bar{r}})^4(dbar{r}^2+ bar{r}^2 dOmega^2 )- big (frac{1-frac{M}{2bar{r}}}{1+frac{M}{2bar{r}}} big) dtau^2$

Las métricas inducidas por estas dos métricas en las hipersuperficies de tiempo propio constante son conformemente planas y singulares en el horizonte.

Coordenadas de Painlevé-Gullstrand-Lemaître $latex (r,theta,varphi, T)$ con $latex dT = dtau + frac{sqrt{frac{2M}{r}}}{1-frac{2M}{r}}dr$ respecto de las de Schwarzschild:

$latex ds^2 = dr^2 + r^2 dOmega^2 + 2 sqrt{frac{2M}{r}}dTdr – (1-frac{2M}{r})dT^2$

que en forma matricial queda:

$latex g=begin{bmatrix} 1 & 0 & 0 & sqrt{frac{2M}{r}} \ 0 & r^2 & 0 & 0 \ 0 & 0 & r^2 sin^2 theta & 0 \ sqrt{frac{2M}{r}} & 0 & 0 & -(1-frac{2M}{r}) end{bmatrix}$

Coordenadas de Eddington-Finkelstein $latex (t, r, theta, varphi)$ con $latex t = tau + 2M ln |frac{r}{2M} – 1|$ respecto de las de Schwarzschild:

$latex ds^2 = frac{1}{1+frac{2M}{r}} dr^2 + r^2 dOmega^2 + frac{4M}{r} dtdr – (1-frac{2M}{r})dt^2$

que en forma matricial queda:

$latex g=begin{bmatrix} 1+frac{2M}{r} & 0 & 0 & frac{2M}{r} \ 0 & r^2 & 0 & 0 \ 0 & 0 & r^2 sin^2 theta & 0 \ frac{2M}{r} & 0 & 0 & -(1-frac{2M}{r}) end{bmatrix}$

Las métricas inducidas por estas dos métricas en las hipersuperficies de tiempo constante son planas y regulares en el horizonte.

En el llibre “Geometria diferencial i relativitat” de J. Girbau tambe comenta les coordenadas de Kruskal-Szekeres $latex (u,v,theta,varphi)$:

$latex ds^2 = frac{32M^3}{r} e^{-frac{r}{2M}} (du^2 – dv^2) + r^2 dOmega^2$

donde

$latex u=sqrt{frac{r}{2M}-1} e^{frac{r}{4M}} cosh frac{tau}{4M}$

y

$latex v=sqrt{frac{r}{2M}-1} e^{frac{r}{4M}} sinh frac{tau}{4M}$

No hay singularidad física en $latex r=2M$, pero hay dos en $latex r=0$.

En segundo lugar, tenemos la métrica de Kerr ($latex J neq 0, Q = 0$) que también podemos escribir en diferentes sistemas de coordenadas:

Coordenadas de Boyer-Lindquist $latex (r,theta,varphi,t)$:

$latex ds^2 = frac{rho^2}{Delta} dr^2 + rho^2 dtheta^2 + tilde{w}^2(dvarphi – wdt)^2 – (frac{rho sqrt{Delta}}{Sigma})^2dt^2$

donde

$latex Delta = r^2 -2Mr + a^2$

$latex rho^2 = r^2 + a^2 cos^2 theta$

$latex Sigma^2 = (r^2 + a^2)^2 – a^2 Delta sin^2 theta$

$latex w = frac{2aMr}{Sigma^2}$

$latex tilde{w} = frac{Sigma sin theta}{rho}$

y siendo $latex a$ el momento angular del BH. Fijando $latex a=0$ obtenemos el BH de Schwarzchild en coordenadas de Schwarzchild.

Coordenadas de Kerr-Schild $latex (r,theta,bar{varphi},bar{t})$:

$latex ds^2 = frac{Z^{2k}-1}{Z-1} dr^2 + rho^2 dtheta^2+ sin^2 theta rho^2 [1+Y(1+Z)] dbar{varphi}^2 – (1-Z) dbar{t}^2+$

$latex +2aepsilon sin^2 theta frac{Z^{k+1}-1}{Z-1}drdbar{varphi} -2 epsilon Z^k dr dbar{t} -2 a sin^2 theta Z dbar{varphi}dbar{t}$

donde

$latex Y = frac{a^2 sin^2 theta}{rho^2}$, $latex Z = frac{2Mr}{rho^2}$

y $latex epsilon = +1(-1)$ regulariza el horizonte futuro (pasado) del BH. La relación con las anteriores coordenadas viene dada por

$latex dbar{varphi} = dvarphi – epsilon frac{a}{Delta} dr$

$latex dbar{t} = dt – epsilon [ frac{1+Y}{1+Y-Z} – frac{1-Z^k}{1-Z} ] dr$

donde $latex Delta$ es la función horizonte, que es cero en el horizonte y hace que la componente $latex g_{tt}$ de la métrica se anule.

Tags: , , , , , , , , , , , ,

Desde un punto de vista astrofísico, un BH es el resultado final de algunos tipos de colapso gravitatorio de objetos estelares concretos. Desde un punto de vista geométrico, un BH es una solución de las ecuaciones de Einstein que contiene una singularidad, curvatura infinita, en la métrica del espacio-tiempo. Llamamos horizonte a la frontera del BH.

Podemos clasificar los agujeros negros en función de su momento angular $latex J$ en estáticos ($latex J=0$) y en rotatorios ($latex J neq 0)$, y en función de su carga eléctrica $latex Q$ en neutros ($latex Q =0$) y cargados ($latex Q neq 0$).

En relatividad numérica hay dos maneras de tratar los BHs y manejar sus singularidades (que pueden ser singularidades de coordenadas, causadas por una inapropiada elección de las mismas, o singularidades físicas reales en el cenro del BH), pues el problema de los métodos numéricos es que no pueden tratar con sus infinitos asociados:

  1. Una manera de evitar el problema es la excisión de la singularidad en el espacio-tiempo. El método onsiste en introducir una frontera interior en el dominio computacional eliminando la región singular. Como el horizonte de un BH representa una frontera causal, está claro que el interior del mismo no puede afectar al exterior. Sin embargo, numéricamente esto no tiene porque ser cierto. Además, tenemos que estar seguros de que la frontera interior está localizada siempre dentro del horizonte, por lo que hay que implementar un buscador del horizonte. Mas concretamente, debe estar dentro del horizonte de sucesos, que depende del espacio-tiempo global. También es complicado aproximar una superficies esférica de escición en una malla cartesiana.
  2. La alternativa es el método de punción: si el comportamiento de una singularidad es analíticamente conocido, podemos factorizarlo de la solución. Una punción en los datos iniciales representa un puente de Einsten-Rosen conectando dos universos asintóticamente planos. La salida del agujero de gusano describe exactamente el campo exterior gravitatorio de un BH. El universo es compactificado mediante una transformación de coordenadas y localizado en el horizonte del BH. La punción representa la infinidad espacial del universo compactificado. La idea intuitiva consiste en tener dos o mas copias de $latex mathbb{R}^3$ con varias esferas eliminadas e identificando los contornos esféricos interiores. De esta manera, varias regiones asintóticamente planas se conectan por puentes. La siguiente figura muestra la idea en una dimensión menos: utilizando dos copias de $latex mathbb{R}^2$ y círculos en lugar de esferas:

En relatividad general, la gravedad actua completamente a través de los potenciales métricos, es decir, a través del espacio-tiempo por el que se mueve el fluido. Por tanto, para manejar la interacción gravitatoria del BH con un fluido basta representar el BH en la estructura del espacio-tiempo.

La tesis doctoral “Relativistic simulations of compact object mergers for nucleonic matter and strange quark matter” de A. Bauswein comenta un nuevo método para trabajar con BHs en relatividad numérica, pués solo mediante una reformulación de las ecuaciones de la métrica podemos manejar la divergencia que tiene lugar en el espacio-tiempo de un BH.

Básicamente, su aproximación viene de combinar diferentes métodos bien establecidos: el método de punción, la descomposición conformal transverse traceless (CTT) de la curvatura extrínseca dando lugar a una nueva formulación de las ecuaciones CFC, y las solución de Bowen-York para las ecuaciones de ligadura del momento.

Esta combinación no es trivial pues hay que describir el significado de la masa de punción durante la evolución, el movimiento del BH, la acreción, la backreaction de las ondas gravitacionales (GW) y la construcción de los datos iniciales.

En el libro “Foundations of Differentiable Manifolds and Lie Groups” de F. Warner, encontramos que un grupo de Lie es una variedad diferenciable $latex G$ en la que tenemos definida una operación $latex G times G longrightarrow G$ tal que:

$latex rho: G times G longrightarrow G ,/, (sigma,tau) mapsto rho(sigma,tau) = sigma tau^{-1}$

es una aplicación diferenciable.

Si fijamos $latex sigma in G$ entonces podemos definir la translación a izquierda por $latex sigma$ y la traslación a derecha por $latex sigma$ de la siguiente manera:

$latex l_sigma(tau) = sigma tau$

$latex r_sigma(tau) = tau sigma$

para cualquier $latex tau in G$. En ambos casos tenemos un difeomorfismo. Si $latex H subset G$ entonces denotamos $latex r_sigma(H)$ por $latex Hsigma$ y $latex l_sigma(H)$ por $latex sigma H$.

El grupo general lineal $latex Gl(n,mathbb{K})$ es el ejemplo mas importante de grupo de Lie. Está formado por los automorfismos de $latex mathbb{K}^n$. Por ser un espacio vectorial, como ya comentamos, es una variedad diferenciable y tiene estructura de grupo con la composición. También puede pensarse como $latex mathcal{M}_n(mathbb{K})$.

Junto al diagrama de clases, que forma parte de los diagramas de modelado estructurado (visión estática del modelo) que ya hemos comentado, existen otros diagramas de modelado de comportamiento (visión dinámica del modelo) ampliamente utilizados como son los diagramas de casos de uso y los diagramas de secuencia.

El diagrama de casos de uso muestra cada una de funcionalidades del sistema modelado en una visión de caja negra de las mismas. Representa al sistema como un rectángulo que contiene tantas elipses como casos de uso diferentes junto con el actor o actores con los que interactuará.

Un ejemplo diagrama de casos de uso de nuestro sistema SPH con un solo caso de uso Simulation podria ser:

El diagrama de secuencia permite modelar cada caso de uso y describe la interacción entre los diferentes objetos de un sistema a través del tiempo para la consecución del mismo. Es por tanto una visión de caja blanca. Se representan los diferentes objetos como rectangulos o conjuntos de rectangulos en función de si son simples o multiples junto con flechas entre los objetos que representan los mensajes entre estos.

Un posible diagrama de secuencia del caso de uso Simulation podria ser:

En el libro “Foundations of differenciable manifolds and Lie groups” de Warner comenta como construir nuevas variedades diferenciables a partir de variedades conocidas.

Dada una variedad $latex (M,mathcal{A})$ donde

$latex mathcal{A} = { varphi_i: U_i longrightarrow mathbb{R}^n }_{i in I}$

y si consideramos un abierto $latex N subset M$ entonces $latex (N,mathcal{A}_N)$ con:

$latex mathcal{A}_N = { varphi_i|_{U_i cap N}: U_i cap N longrightarrow mathbb{R}^n }_{i in I}$

es una variedad diferenciable.

Si $latex (M, mathcal{A}_1)$ y $latex (N,mathcal{A}_2)$ son dos variedades diferenciables de dimensión $latex p$ y $latex q$ respectivamente donde

$latex mathcal{A}_1 = { varphi_i: U_i longrightarrow mathbb{R}^p}_{i in I}$

y

$latex mathcal{A}_2 = { psi_j: V_j longrightarrow mathbb{R}^q}_{j in J}$

entonces $latex (M times N, mathcal{A}_1 times mathcal{A}_2)$ es una variedad diferenciable de dimensión $latex p+q$ donde:

$latex mathcal{A}_1 times mathcal{A}_2 = { varphi_i times psi_j: U_i times V_j longrightarrow mathbb{R}^p times mathbb{R}^q}_{i in I, j in J}$

Esta propiedad nos permite demostrar que el toro es una variedad diferenciable. Si consideramos las variedades $latex S(r)$ y $latex S(R)$ con $latex r<R$ entonces $latex T(R,r) = S(R) times S(r)$.

Finalmente, si $latex (M,mathcal{A})$ es una variedad diferenciable de dimensión $latex n$ con $latex mathcal{A} = { varphi_i: U_i longrightarrow mathbb{R}^n}_{i in I}$ y sea $latex f: Mlongrightarrow N$ una aplicación biyectiva. Si

$latex mathcal{A}_f = { varphi_i circ f^{-1}: f(U_i) longrightarrow mathbb{R}^n}_{i in I}$

entonces $latex (N, mathcal{A}_f)$ es una variedad diferenciable de dimensión $latex n$.

Esta propiedad nos permite decir que todo espacio vectorial de dimensión finita es una variedad diferenciable. Efectivamente, sea $latex V$ un espacio vectorial de dimensión $latex n$ sobre $latex mathbb{R}$ y $latex {u_1,ldots,u_n}$ una base. Entonces existe un único isomorfismo lineal $latex f: mathbb{R}^n longrightarrow V$ tal que $latex f(e_i) = u_i$ donde $latex { e_1, ldots, e_n}$ es la base canónica de $latex mathbb{R}^n$.  Entonces $latex (V,mathcal{A_f})$ es una variedad diferencial de dimensión $latex n$.

VisIt es una herramienta open source que nos permite visualizar y analizar conjuntos de datos extremadamente grandes, del orden de tera y peta, en multiples plataformas. Podemos visualizar rápidamente nuestros datos, animarlos, manipularlos y almacenar los resultados obtenidos.

Algunas características interesantes son:

  • Tipos estandar de gráficos: Curve plot, Mesh plot, Contour plot, Surface plot, Vector plot,  Tensor plot, Volume plot, etc.
  • Podemos trabajar en 1D, 2D, 3D y variando en tiempo.
  • Permite definir diferentes tipos de mallas: rectilineas, curvilineas, desestructuradas, puntuales y AMR (Adaptive Mesh Refinement), etc.
  • Manipulación de datos (slicing, clipping, project, etc.) e interrogaciones (analisis comparativo, debugging, etc.).
  • Opciones para anotaciones, iluminación y rendering.
  • Podemos trabajar con datos escalares, vectoriales y tensoriales.
  • Permite paralelizaciones.

Para visualizar los datos podemos utilizar VisIt como aplicación mediante fichero ($latex approx$ 100 formatos diferentes), base de datos (se pueden desarrollar nuevos plug-ins) o como librería mediante código. También podemos crear animaciones mediante flipbook, keyframing o scripting. Podemos trabajar en local, en remoto o utilizando la arquitectura cliente/servidor.

VisIt se compone de cuatro componentes:

  • Graphical User Interface (GUI): se ejecuta en local y permite, entre otras cosas, seleccionar los ficheros, crear graficos, fijar atributos, etc.
  • Viewer: se lanza en local y es donde se muestran las visualizaciones con las que podemos interactuar con el ratón.
  • Database server: puede lanzarse en remoto y permite el acceso a los datos.
  • Compute engine / Parallel compute server: puede lanzarse en remoto y es el encargado de generar los gráficos.

A continuación una pequeña animación flipbook que hemos creado con datos que ofrece VisIt en su documentación.

[youtube http://www.youtube.com/watch?v=EaCWZGLwwZs?rel=0]

Los ficheros .visit, en este caso wave.visit, son ficheros de texto que contienen los nombres de todos los ficheros que guardan, cada uno, el estado de la simulación en un instante de tiempo determinado: wave0010.silo, wave0010.silo, ..., wave0700.silo. En este caso son .silo que es un formato propio de VisIt.

Tags: , , , ,

En el libro “Geometria Diferencial i Relativitat” de J. Girbau aparecen algunos ejemplos concretos de variedades diferenciables.

Para empezar, $latex { (mathbb{R}^n, id_{mathbb{R}^n}) }$ es un atlas que dota al espacio euclídeo $latex mathbb{R}^n$ de estructura de variedad diferenciable de dimensión $latex n$ y clase $latex C^infty$.

También podemos considerar la esfera $latex S^n(r)$ ($latex n$-esfera, hiperesfera) de radio $latex r$ centrada en el origen de $latex mathbb{R}^{n+1}$. Si llamamos $latex N$ al polo norte, de coordenadas $latex (0,0,ldots,0,r)$, y $latex S$ al polo sur, de coordenadas $latex (0,0,ldots,0,-r)$, sean $latex U_1 = S^n(r)-{N}$ y $latex U_2=S^n(r)-{S}$ y sean $latex varphi_1:U_1 longrightarrow mathbb{R}^n$ la proyección estereográfica desde $latex N$ sobre el plano $latex x^{n+1}=0$ y $latex varphi_2:U_2 longrightarrow mathbb{R}^n$ la proyección estereográfica desde $latex S$ sobre el plano del ecuador. Se puede comprobar que $latex { (U_1, varphi_1), (U_2, varphi_2) }$ es un atlas $latex C^infty$ de dimensión $latex n$ que dota a $latex S^{n}(r)$ de una estructura de variedad diferenciable.

Otro ejemplo concreto de variedad diferenciable es el plano proyectivo real $latex mathbb{R}P^n$. Consideramos en $latex mathbb{R}^{n+1} – { 0 }$ la relación de equivalencia:

$latex x sim y Leftrightarrow exists lambda in mathbb{R}: x = lambda y$.

De esta manera,

$latex mathbb{R}P^n := frac{mathbb{R}^{n+1} – { 0 }}{sim}$

Es fácil comprobar que $latex { (U_i, varphi_i)}, i=0,ldots,n$ donde

$latex U_i = { [ x^0 :ldots : x^n ], x_i neq 0 }, i=0,ldots,n$,

siendo $latex [ x^0: ldots :x^n ]$ coordenadas homogéneas, y

$latex varphi_i: U_i longrightarrow mathbb{R}^n ,/, [x^0:ldots:x^n] mapsto (frac{x^0}{x^i},ldots,hat{frac{x^i}{x^i}},ldots,frac{x^n}{x^i})$

donde $latex hat{frac{x^i}{x^i}}$ significa que falta la componente $latex i$-ésima, dota a $latex mathbb{R}P^n$ de una estructura de variedad diferenciable. Este ejemplo es de naturaleza diferente del anterior ya que $latex mathbb{R}P^n$ no se presenta dentro de ningún $latex mathbb{R}^n$.

En la tesis doctoral “Relativistic simulations of compact object mergers for nucleonic matter and strange quark matter” de A. Bauswein comenta la implementación numérica del modelo físico para simular colisiones de objetos compactas.

La implementación numérica consta de dos partes: una que resuelve las ecuaciones de Euler tridimensionales relativistas que nos dan la evolución hidrodinámica del sistema y la otra que nos da una solución de las ecuaciones de campo de Einstein que nos proporcionan el campo gravitatorio donde se mueve el fluido.

Como ya comentamos en este post, utilizando SPH las ecuaciones de la hidrodinámica se tranforman en un conjunto de ODEs que podemos resolver mediante el método Runge-Kutta de cuarto orden. Las derivadas de los potenciales de la métrica se evaluan sobre una malla superpuesta y se mapean en las partículas. Se utiliza un paso de tiempo adaptativo para satisfacer la condición de Courant-Friedrichs-Levi. Todas las cantidades se representan por sus valores en las partículas, que corresponde con una visión de la hidrodinámica Lagrangiana, y la evolución temporal de las cantidades hidrodinámicas se calcula en la partículas.

Se añade un esquema de viscosidad artificial dependiente del tiempo para manejar choques hidrodinámicos. Este modelo resuelve un problema de Riemann local dado por los dos estados de partículas vecinas y añadiendo la correspondiente contribución a las ecuaciones de conservación del momento y de la energía.

A diferencia del SPH Newtoniano, la gravedad no se puede implementar dentro del esquema de partículas. Debemos resolver las ecuaciones de campo de Einstein, para lo que se utiliza el formalismo 3+1 o ADM que consiste en foliar el espacio-tiempo en hipersuperfícies espaciales con coordenada de tiempo constante. En esta aproximación, las ecuaciones de campo de Einstein se dividen en un conjunto de ecuaciones de ligadura y un conjunto de ecuaciones de evolución formulando de esta manera un problema de Cauchy.

Dos variedades diferenciales $latex M$ y $latex N$ son difeomorfas si existe un difeomorfismo entre ellas. Sin considerar ninguna estructura adicional, dos variedades difeomorfas son equivalentes.

Sea $latex f:M longrightarrow N$ una aplicación entre variedades de clase $latex C^k$.

Diremos que $latex f$ es un difeomorfismo de clase $latex C^k$ si es bijectiva, diferenciable y $latex f^{-1}$ también es diferenciable de clase $latex C^k$.

Diremos que $latex f$ es diferenciable de clase $latex C^k$ si lo es en cada punto de $latex M$. Usualmente, como en variedades, utilizamos diferenciable para diferenciable de clase $latex C^infty$.

Diremos que $latex f$ es diferenciable en $latex m in M$ de clase $latex C^k$  si existen cartas adaptadas $latex (U,Phi)$ de $latex M$ y $latex (V,Psi)$ de $latex N$ tales que $latex m in U$ y su representación local $latex bar{f}$ es diferenciable en $latex Phi(m)$ de clase $latex C^k$.

Dos cartas $latex (U,Phi)$ de $latex M$ y $latex (V,Psi)$ de $latex N$ estan adaptadas a $latex f$ si $latex f(U) subset V$. Llamamos representación local de $latex f$ en las cartas $latex (U,Phi)$ y $latex (V,Psi)$ a $latex bar{f} = Psi circ f circ Phi^{-1}: A longrightarrow B$.

A continuación se muestra una imagen con todos los fichero generados de forma automática:

y el código que contienen algunos de estos ficheros:

Cada instanciación de una clase define un objeto: Persona es una clase, Jose es un objeto, una instancia de Persona. La clase es una plantilla que permite crear objetos.

Una clase almacena un estado: un conjunto de atributos que guardan información del objeto; y un comportamiento: un conjunto de métodos que nos permiten trabajar con el objeto.

A continuación, nuestro diagrama de clases con una muestra de atributos, métodos y relaciones:

En el libro “Geometria Diferencial i Relativitat” de J. Girbau, encontramos que una variedad diferenciable de dimensión $latex n$ y de clase $latex C^k$ es un par $latex (M,mathcal{A})$ formado por un conjunto $latex M$ y una estructura diferenciable $latex mathcal{A}$ de dimensión $latex n$ y clase $latex C^k$ (de aquí en adelante supondremos dimensión $latex n$ y clase $latex C^k$) y tal que $latex M$ con la topología asociada es Hausdorff y 2AN.

Si no hay posibilidad de confusión se denota la variedad simplemente por $latex M$ y es usual utilizar la expresión variedad diferenciable cuando es de clase $latex C^infty$.

Llamamos estructura diferenciable a cada una de las clases de equivalencia por la relación $latex C^kmbox{-compatibles}$: dos atlas $latex A_1$ y $latex A_2$ de clase $latex C^k$ son $latex C^kmbox{-compatibles}$ si $latex A_1 cup A_2$ también es un atlas de clase $latex C^k$. Como dos atlas $latex C^kmbox{-compatibles}$ definen la misma topología en $latex M$, podemos hablar de topología asociada a una estructura diferenciable.

Un atlas de clase $latex C^k$ es una familia de aplicaciones biyectivas $latex mathcal{A}={ varphi_i: U_i rightarrow A_i}_{i in I}$ llamadas cartas locales, donde $latex I$ es un conjunto de índices, $latex U_i subset M$ abierto y $latex A_i subset mathbb{R}^n$ abierto cumpliendo:

  1. $latex cup_{i in I} U_i = M$
  2. $latex forall i,j in I$ tal que $latex U_i cap U_j neq emptyset$, entonces $latex varphi_i(U_i cap U_j)$ y $latex varphi_j(U_i cap U_j)$ son abiertos y las aplicaciones $latex varphi_j circ varphi_i^{-1}:varphi_i(U_i cap U_j) longrightarrow varphi_j(U_i cap U_j)$ son $latex C^k$

Si $latex mathcal{A}$ es un atlas sobre un conjunto $latex M$ entonces

$latex mathcal{B}_{mathcal{A}} = { varphi_{i}^{-1}(W) : i in I, W subset mathbb{R}^n mbox{ abierto}}$

es una base de una topología en $latex M$.

Una variedad de Riemann es una variedad diferenciable en la que tenemos definida una métrica: un campo tensorial diferenciable $latex g$ dos veces covariante, un tensor $latex (0,2)$, simétrico y definido positivo. A $latex g$ se le llama métrica de Riemann.

Si $latex (M,g)$ es una variedad de Riemann y $latex x(t)$ es una curva diferenciable sobre $latex M$, entonces la longitud de $latex x(t)$ entre dos puntos $latex x(a)$ y $latex x(b)$ se define como:

$latex int_a^b sqrt{g(dot{x(t)},dot{x(t)})} dt$

Una variedad pseudoriemanniana es aquella en la que la métrica no cumple la propiedad de ser definida positiva y basta con que sea no degenerada.

Una variedad de Lorentz es aquella cuya métrica tiene signatura $latex (n-1, 1)$. Son especialmente importantes por sus aplicaciones en la física: la teoria de la relatividad general modela el espacio-tiempo mediante una variedad de Lorentz de dimensión $latex n=4$ y signatura $latex (3,1)$.

Cuando modelamos sistemas, ciertos componentes pueden estar relacionados con otros y estas interrelaciones deben ser modeladas. La relación de asociación es el tipo de relación básico entre dos clases. Se representa mediante una linea donde podemos incluir roles, cardinalidades, dirección y restricciones. Suele implementarse como variables de instancia de una clase en otra.

Con respecto a la implementacion, podriamos tener algo como:

//asociación: obj1 y obj2 tienen una relación.
public class Obj1 { void func(Obj2 obj2) {} };

Una relación de agregación es un tipo particular de asociación aparecido después y que permite describir cuando una clase está formada por otras clases. En las agregaciones, el ciclo de vida de una parte es independiente del ciclo de vida del todo, es decir, que si el padre se elimina, no se eliminan sus hijos. En el diagrama aparecen con una punta de flecha en forma de diamante negro en la clase padre. El código podría ser:

//agregación: obj1 posee a obj2 que se lo han prestado: cuando obj1 muere obj2 puede vivir.
public class Obj1 { private Obj2 obj2; Obj1(Obj2 obj2) { this.obj2 = obj2; } }

Una relación de composición permite describir el mismo tipo de relación que las agregaciones pero, por el contrario, la dependencia de sus ciclos de vida es mas fuerte: si el padre se elimina también se eliminan todos sus hijos, aunque una parte puede ser eliminada sin eliminar el todo. En el diagrama aparecen con una punta de flecha en forma de diamante blanco en la clase padre. En código:

//composición: obj1 posee a obj2 y es responsable de su tiempo de vida: cuando obj1 muere obj2 también muere.
public class Obj1 { private Obj2 obj2 = new Obj2(); }

La lista de clases obtenida junto con sus relaciones de generalización:

Imagen

(Diagrama realizado con Visual Paradigm for UML)

El diagrama de clases muestra los bloques con los que construiremos un sistema orientado a objetos. Describen la vida estática del modelo, mostrando que atributos y comportamientos ocurren sin entrar a detallar de que manera se realizan. También muestran las relaciones existentes entre las clases: generalizaciones que reflejan herencias, agregraciones para reflejar composición o uso y asociaciones para mostrar conexiones.

Para representar cada clase utilizamos un rectángulo con tres compartimentos. En el primero va el nombre de la clases, en el segundo sus atributos que pueden ir acompañados de su valor inicial y en el último los métodos junto a sus parámetros. La visibilidad se indica mediante los simbolos – para indicar que un atributo o método es privado, y + para indicar que un atributo o un método es público.

Una relación de generalización indica herencia. Se representa mediante una punta de flecha triangular en la clase padre que generaliza una clase hijo. Un objeto instanciado de la clase hijo heredará de manera implícita los atributos y métodos de su clase padre. En nombre de las clases abstractas va en cursiva.

Volviendo a nuestro diagrama particular, las clases abstractas las hemos indicado temporalmente mediante un borde discontinuo y los colores indican posibles paquetes: verde para clases específicas de SPH, amarillo para numericalMethods, azul para astrophysics y rojo para dataStructures.

Vamos a hacer una lista con todas las clases que creemos que vamos a necesitar para la implementación del método SPH y aplicarlo a simulaciones de problemas astrofísicos. A continuación pensaremos en sus atributos y métodos. Finalmente, interrelacionaremos las clases entre si. Para todo ello, nos apoyaremos en el Visual Paradigm. Empecemos:

  1. Dispondremos de un stage donde ocurre la acción de nuestra simulación. Este contendrá una clase abstracta spatialDecomposition que podremos implementar posteriormente en las clases derivadas bsptree, octree, kdtree, linkedList y que utilizaremos para selección de partículas vecinas incluso con $latex h$ variable. También necesitamos una clase mesh/grid para ir resolviendo las ecuaciones de campo de Einstein, que constituyen un problema PDE elíptico, e ir obteniendo el campo gravitacional en el que se mueve nuestro fluido. Por tanto, también necesitaremos un ellipticSolver que derivaremos en spectralSolver o finiteDifferenceSolver.
  2. Necesitamos una clases fluid que almacenará toda la información referente a un fluido .
  3. Los fluidos nos permitiran modelar compactObjects como blackHole, neutronStar, strangeStar… y en donde dispondremos de una operacion merge.
  4. Y también la clase particle que permitirá almacenar la información de cada una de las partículas de nuestros fluidos.
  5. La clase kernel nos permitirá evaluar la función kernel. Como ya comentamos en un post, para ofrecer la posibilidad de utilizar diferentes funciones kernel, la clase kernel será abstracta ofreciendo métodos virtuales que implementaremos en las clases derivadas gaussianKernel, quadraticKernel, cubicKernel, quinticKernel, tensorialKernel.
  6. Necesitamos también una clase odeSolver que utilizaremos para resolver las distintas ODE que aparecen al “discretizar” las ecuaciones de la hidrodinámica. Como tenemos diferentes posibilidades, la implementaremos en una clase abstracta cuyos métodos virtuales los implementaran rungeKutta2OdeSolver, rungeKutta4OdeSolver, velocityStormerVerletOdeSolver, leapfrogOdeSolver, modifiedEulerOdeSolver y que nos permitirá utilizar cualquiera de ellos indistintamente.
  7. La clase abstracta eos nos permitirá trabajar con diferentes ecuaciones de estado.
  8. Desde el punto de vista matemático, necesitamos definir una clase tensor y sus casos particulares scalar, vector y matrix.

El problema de Riemann es un caso especial de problema de valor inicial (IVP) en el que la PDE es:

$latex u_t + a u_x = 0$

y la condicion inicial (IC):

$latex u(x,0) = u_0(x) = begin{cases} u_Lmbox{ si } x < 0 \ u_R mbox{ si } x > 0 end{cases}$

donde $latex u_L$ y $latex u_R$ son dos valores constantes, de manera que tenemos una discontinuidad en $latex x=0$. En este caso, la solución es sencilla:

$latex u(x,t) = u_0(x-at) = begin{cases} u_L mbox{ si } x-at < 0\ u_R mbox{ si } x-at > 0 end{cases}$.

Podemos extender el problema a un conjunto de $latex m$ PDEs hiperbólicas:

$latex U_t + AU_x = 0$ con $latex -infty < x < infty, t>0$

donde la matriz $latex A$ es de coeficientes constantes. Al asumir la hiperbolicidad de $latex A$, tenemos $latex m$ valores propios reales $latex lambda_i$ y $latex m$ vectores propios independientes $latex K^{(i)}$. En este caso, la IC se escribe como:

$latex U(x,0) = U^{(0)}(x) = begin{cases} U_L, x<0\ U_R, x>0 end{cases}$

Los Riemann solvers son métodos numéricos que permiten resolver el problema de Riemann.

La siguiente operación en física:

$latex a+b+c = d+e$

solamente la podremos realizar si todas las variables están en las mismas unidades. En física solemos medir longitudes ($latex L$), masas ($latex M$), tiempos ($latex T$), temperaturas ($latex K$), etc. La idea es disponer de unas unidades de manera que ciertas constantes universales tomen el valor $latex 1$ simplificando así algunas ecuaciones físicas.

Dada la ecuación $latex v=frac{d}{dt}x$ podemos escribir la correspondiente ecuación de dimensión $latex [v] = frac{L}{T} = L T^{-1}$. De la misma manera, nos quedaria $latex [a] = L T^{-2}$ para las aceleraciones.

Sabemos que la velocidad de la luz en el vacío $latex c = 3 times 10^{10} frac{cm}{s}$ en el sistema cegesimal (CGS: cm, g, s) y su ecuación dimensional es $latex [c] = LT^{-1}$. Si queremos conseguir $latex c=1$, ¿qué unidades de longitud y tiempo necesitamos? Suponiendo que la longitud la continuamos midiendo en $latex cm$, aunque lo denotamos ahora con $latex uL$, vamos a definir una unidad de tiempo $latex uT$ de manera que $latex c=1$:

$latex 3 times 10^{10} frac{uL}{s} times frac{1}{x} frac{s}{uT} = 1 frac{uL}{uT} Leftrightarrow 1 uT = frac{1}{3 times 10^{10}} s$

Por lo que en $latex uL$ y $latex uT$ hemos conseguido que $latex c=1$.

Tenemos que la constante de gravitación universal $latex G = 6.67 times 10^{-8} frac{cm^3}{g s^2}$ en el sistema CGS, con $latex [G] = L^3 M^{-1} T^{-2}$. Si hacemos:

$latex 6.67 times 10^{-8} frac{uL^3}{g s^2} times frac{1}{(3 times 10^{10})^2} frac{s^2}{uT^2} times frac{1}{x} frac{g}{uM}=1 frac{uL^3}{uM uT^2} $

tenemos que $latex x = frac{2.2}{3} times 10^{-28}$ y nos queda que $latex 1 g = frac{2.2}{3} times 10^{-28} uM$. Así pués, $latex G=1$ en estas nuevas unidades.

Si prodecemos de esta manera hasta conseguir que  $latex c = G = k_B = 1$, donde $latex k_B$ es la constante de Boltzmann, obtenemos las unidades geometrizadas. Si lo hacemos para $latex c =k_B=h=1$, donde $latex h$ es la constante de Planck, obtenemos las unidades naturales.

En C++ podemos definir clases abstractas. Son clases pensadas para definir un comportamiento, es decir, un conjunto de métodos sin implementación. Sus clases derivadas están obligadas a implementar estos métodos, siempre i cuando queramos que sean clases concretas instanciables. Los métodos sin implementación de las clases abstractas reciben el nombre de métodos virtuales.

En el siguiente diagrama de clases UML podemos ver dos pequeños ejemplos de como podríamos utilizar las clases abstractas. En el primer caso, tenemos la clase abstracta kernel para implemetar la función kernel del método SPH. Uno de sus métodos debe devolver un valor a partir de los parámetros $latex d$ y $latex h$. La clase kernel define la signatura del método virtual calcular. Las clases derivadas quadraticKernel, cubicKernel y quinticKernel lo implementan de manera diferente. La clase que lo utiliza tiene como atributo un kernel abstracto, de manera que le podremos asignar cualquiera de las derivaciones concretas. Lo mismo ocurre con la clase abstracta Eos que tiene un método que devuelve la presión. Tenemos dos clases derivadas que lo implementan, cada cual de una manera. La clase cliente con un atributo Eos podrá instanciar cualquiera de las dos clases derivadas.

El siguiente código, generado automáticamente por Visual Paradigm, muestra como queda implementado el UML2 de kernel i cubicKernel en C++:

kernel.h

[sourcecode languaje=”cpp”]
namespace example {

class kernel {

public:
virtual double calculate(double d, double h) = 0;

};

}
[/sourcecode]

cubicKernel.h

[sourcecode languaje=”cpp”]
namespace example {

class cubicKernel : example::kernel {

public:
double calculate(double d, double h);

};

}
[/sourcecode]

cubicKernel.cpp

[sourcecode languaje=”cpp”]

#include "cubicKernel.h"

double example::cubicKernel::calculate(double d, double h) {

//aquí va la implementación concreta del kernel cúbico
throw "Not yet implemented";

}
[/sourcecode]

LORENE (Langage Objet pour la RElativité NumériquE), tal y como aparece en la página web de la sección Meudon del Observatorio de Paris, es un conjunto de clases C++ pensada para resolver varios problemas que aparecen en relatividad numérica, y de manera mas amplia en astrofísica computacional. Proporciona una serie de herramientas para resolver ecuaciones en derivadas parciales mediante métodos espectrales multidominio.

Las clases de LORENE implementan estructuras básicas, como vectores y matrices, objetos matemáticos mas abstractos, como tensores, y objetos astrofísicos, como estrellas o agujeros negros.

Es un software libre distribuido bajo licencia GNU General Public License.

Utilizando software de ingenieria inversa, podemos obtener el diagrama de clases de LORENE. De esta manera, podemos ver las relaciones existentes entre estas y comprender mejor el software desarrollado.

El siguiente diagrama, obtenido mediante Visual Paradigm for UML, muestra dos subdiagramas de clases, uno relacionado con agujeros negros y otro con estrellas (algunas clases se muestran expandidas con atributos y metodos mientras que otras se muestran sin expandir):

En [Monaghan 2005] se revisan la teoría y las aplicaciones del SPH desde su aparición hasta la actualidad.

El artículo empieza comentando cinco de los puntos fuertes del método SPH:

  1. La advección se trata de manera exacta: si a una partícula se le asigna un color, y se le especifica una velocidad, el transporte del color por la partícula del sistema es exacto.
  2. Cuando trabajamos con mas de un material, cada uno se describe mediante su propio conjunto de partículas de manera que los problemas de interfaz resultan triviales.
  3. Cierra la brecha entre el contínuo y la fragmentación de una forma natural.
  4. La resolución puede depender tanto de la posición como del tiempo.
  5. SPH tiene la ventaja computacional de que la computación ocurre solo donde está realmente la materia, con la consecuente reducción de calculo y almacenamiento.

En [Monaghan 1992] comenta algunos detalles de implementación del método SPH.

La inicialización de un cálculo SPH comienza especificando la masa, posición, velocidad y energía térmica de cada partícula. Si calculamos la densidad mediante la forma alternativa:

$latex frac{d}{dt}rho_a = sum_b m_b v_{ab} nabla_a W_{ab}$

donde $latex v_{ab} = v_a – v_b$ y que tiene algunas ventajas (evita la oscilación de los bordes del fluido y es computacionalmente mas eficiente el calculo de la tasa de cambio de todas las variables físicas), necesitamos asignar la densidad inicial de cada partícula. También la mezcla de elementos de cada partícula puede ser necesaria.

Para inicializar las partículas, es conveniente situarlas en una malla regular, que puede ser una malla Cartesiana regular con igual espaciado o una malla cúbica centrada en el cuerpo en tres dimensiones (mejor representación de las integrales por sumas, mejor detección de vecinos y mejor distribución de partículas en caso de compresión planar). Además, si el tamaño de celda asociado con una partícula $latex a$ es $latex Delta V_a$ entonces $latex m_a$ se toma como $latex rho Delta V_a$.

Si utilizamos la misma $latex h$ para cada partícula, entonces basta trabajar con listas enlazadas. Si el kernel esta basado en splines, las celdas enlazadas por la lista deben ser de $latex 2h$ de ancho de manera que solo las celdas vecinas pueden contribuir en las partículas en una celda dada.

Si cada partícula tiene su propia $latex h$ entonces el calculo de las sumas del SPH pueden formar parte de un calculo en arbol que es el método natural para calcular las fuerzas autogravitantes de un conjunto de partículas. Este codigo en arbol puede ser vectorizado.

Trabajando en GRMHD nos encontramos el siguiente problema de Cauchy:

  1. Ecuaciones de Einstein: $latex (gamma,K)_{t=0} longrightarrow (gamma,K)_t$
  2. Hidrodinámica: $latex (rho,q,tau)_{t=0} longrightarrow (rho,q,tau)_t$
  3. Electromagetismo: $latex (E,B)_{t=0} longrightarrow (E,B)_t$

donde $latex gamma$ es la métrica y $latex K$ la curvatura extrinseca; $latex rho$ es la densidad de masa, $latex q$ la densidad de momento, $latex tau$ la densidad de energia; $latex E$ es el campo de fuerza eléctrica y $latex B$ la inducción magnética.

En [Rosswog 2009] comenta algunos métodos para ODEs, ya que una de las características relevantes de SPH es que transforma las PDEs de la hidrodinámica en un sistema de ODEs y hacer una simulación en SPH no es mas que resolver dicho sistema.

Las ODEs en astrofísica suelen ser de segundo orden pero podemos transformarlas en un sistema de ODEs de primer orden introduciendo como nuevas variables a las velocidades. De esta manera obtenemos $latex n$ ecuaciones acopladas de la forma:

$latex frac{d}{dt}y_i = f_i(t,y_1,ldots,y_n)$ con $latex i = 1,ldots,n$

Tenemos, en este caso, dos tipos de problemas: los problemas de contorno y los problemas de valor inicial. Este último caso, que es el que nos interesa, se trata de determinar el valor de $latex y_i^n$ en $latex t^n$ a partir de la condición inicial $latex y_i^0$ en $latex t^0$.

El método de Euler es el método de integración pedagógico. La solución avanza de $latex t^n$ a $latex t^{n+1} equiv t^n + Delta t$ según:

$latex vec{y}^{n+1} = vec{y}^n + Delta t vec{f}(t^n,vec{y}^n)$

donde $latex vec{y}^{n+1} = vec{y}(t^{n+1})$ y $latex Delta t$ es el paso de tiempo escogido. Es un método explícito de primer orden y es asimétrico en tiempo, pues la derivada es correcta al principio del paso de tiempo pero no al final. De ahí el nombre de Euler hacia adelante. En el caso de posiciones, velocidades y aceleraciones, tenemos:

$latex vec{x}^{n+1} = vec{x}^n + Delta t vec{v}^n$

$latex vec{v}^{n+1} = vec{v}^n + Delta t vec{a}^n$

La posibilidad contraria es:

$latex vec{y}^{n+1} = vec{y}^n + Delta t vec{f}(t^{n+1},vec{y}^{n+1})$

llamado el método de Euler hacia atrás, que es un método implícito. Todos los siguientes métodos que comentaremos son explícitos.

El siguiente método es el Euler modificado que utiliza el valor de la derivada al principio y al final del paso de tiempo:

$latex vec{y}^{n+1} = vec{y}^n + Delta t frac{vec{f}^n + vec{f}^{n+1,p}}{2}$

que es un método de segundo orden. En el caso de posiciones, velocidades y aceleraciones tenemos:

$latex vec{x}^{n+1} = vec{x}^n + Delta t frac{vec{v}^n + vec{v}^{n+1,p}}{2}$

$latex vec{v}^{n+1} = vec{v}^n + Delta t frac{vec{a}^n + vec{a}^{n+1,p}}{2}$

En el método de Runge-Kutta de segundo orden tomamos la derivada en el punto medio de los intervalos de tiempo tomados:

  1. Tomamos como paso $latex frac{Delta t}{2}$, calculamos $latex vec{y}^{mid} = vec{y}^n + frac{1}{2} Delta t vec{f}(t^n,vec{y}^n)$ y evaluamos la derivada en ese punto $latex vec{f}^{mid} = vec{f}(t^{frac{n+1}{2}},vec{y}^{mid})$
  2. Tomamos el paso entero con esta derivada: $latex vec{y}^{n+1} = vec{y}^n + Delta t vec{f}^{mid}$

En el caso de posiciones, velocidades y aceleraciones tenemos:

$latex vec{x}^{n+1} = vec{x}^n + Delta t vec{v}^{mid}$

$latex vec{v}^{n+1} = vec{v}^n + Delta t vec{a}^{mid}$

Runge-Kutta de cuarto orden:

$latex vec{k_1} = Delta t vec{f}(t^n,vec{y}^n)$

$latex vec{k_2} = Delta t vec{f}(t^n + frac{Delta t}{2},vec{y}^n + frac{vec{k}_1}{2})$

$latex vec{k_3} = Delta t vec{f}(t^n + frac{Delta t}{2},vec{y}^n + frac{vec{k}_2}{2})$

$latex vec{k_4} = Delta t vec{f}(t^n + Delta t,vec{y}^n + vec{k}_3)$

$latex vec{y}^{n+1} = vec{y}^n + frac{1}{6}(vec{k}_1+ 2 vec{k}_2 + 2 vec{k}_3 + vec{k}_4)$

Velocity Stormer-Verlet:

En el caso de posiciones, velocidades y aceleraciones tenemos:

$latex vec{x}^{n+1} = vec{x}^n + vec{v}^n Delta t + frac{1}{2} vec{a}^n Delta t^2$

$latex vec{v}^{n+1} = vec{v}^n + Delta t frac{vec{a}^n + vec{a}^{n+1}}{2}$

Dado que la conservación exacta es en nuestro puede ser mas importante  que el orden alto, aconseja el último método.

Para un observador inicial, si denotamos con $latex vec{E}$ al campo eléctrico, $latex vec{B}$ al campo magnético, $latex rho$ a la densidad de carga y $latex vec{J}$ a la densidad de corriente, entonces tenemos las ecuaciones de Maxwell:

$latex nabla cdot vec{E} = rho$

$latex nabla times vec{E} + vec{B}_t = 0$

$latex nabla cdot vec{B} = 0$

$latex nabla times vec{B} – vec{E}_t = vec{J}$

y la ecuación de continuidad o de conservación de carga:

$latex rho_t + nabla cdot vec{J} = 0$

Si el observador inercial se encuentra en el espacio-tiempo de Minkowski, las ecuaciones de Maxwell se expresan como dos ecuaciones de ligadura:

$latex nabla cdot vec{E} = rho$

$latex nabla cdot vec{B} = 0$

y seis ecuaciones de evolución:

$latex vec{E}_t = nabla times vec{B} – vec{J}$

$latex vec{B}_t = -nabla times vec{E}$

Si en un instante $latex t=t_0$ se cumplen las ecuaciones de ligadura y si la carga eléctrica se conserva en un entorno de $latex t=t_0$,

$latex rho_t + nabla cdot vec{J} = 0$,

entonces las ecuaciones de ligadura se cumplen en ese entorno (como consecuencia de las ecuaciones de evolución).

Tags: ,

En [Rosswog 2009] tenemos las ecuaciones de la hidrodinámica en forma Lagrangiana discretizadas y en su forma mas básica. Las partículas avanzaran en el tiempo siguiendo las siguientes ecuaciones:

Para empezar, no hay necesidad de resolver la ecuación de continuidad ya que la masa de las partículas permanece fija. Podemos obtener las densidades mediante:

$latex rho_a = sum_b m_b W_{ab}$

La ecuación del momento queda:

$latex frac{d}{dt}vec{v}_a = – sum_b m_b (frac{P_a}{rho_a^2} + frac{P_b}{rho_b^2} + Pi_{ab} ) nabla_a W_{ab} $

La ecuación de evolución para la energía interna específica puede escribirse como:

$latex frac{d}{dt} u_a = sum_b m_b (frac{P_a}{rho_a^2} + frac{1}{2}Pi_{ab}) vec{v}_{ab} nabla_b W_{ab} $

Rosswog las llama “vanilla ice” SPH.

Tags: , , , , , ,

En las Jornadas sobre los problemas del milenio celebradas en Barcelona del 1 al 3 de junio de 2011, Diego Córdoba dió una charla sobre el problema Clay de las ecuaciones de Navier-Stokes para la que escribió estas notas.

Un fluido es ideal si es incompresible, homogéneo  y perfecto.

La idea del problema consiste en determinar si:

…un fluido incompresible con energía finita puede desarrollar singularidades en tiempo finitos.

Mas formalmente, si consideramos un fluido viscoso ($latex nu > 0$), homogéneo ($latex rho = 1$) e incompresible ($latex nabla cdot u = 0$), tenemos:

$latex u_t + u cdot nabla u = – nabla p + nu Delta u + f$ con $latex nu >0, x in mathbb{R}^3, t geq 0$

$latex nabla cdot u = 0$

$latex u(x,0) = u_0$

donde cada partícula en el tiempo $latex t$ está en la posición $latex x = (x_1,x_2,x_3)$ del dominio que ocupa el fluido $latex Omega subset mathbb{R}^3$, $latex u(x,t) = (u_1(x,t), u_2(x,t), u_3(x,t))$ es el campo de velocidades, $latex p = p(x,t)$ son las presiones en el seno del fluido y $latex rho = rho(x,t)$ es la densidad. Además, $latex nu = cte geq 0$ es la viscosidad y $latex f$ la fuerza externa.

La fuerza exterior  debe verificar:

$latex |partial_x^alphapartial_t^m f| leq C_{alpha,k,m}(1+|x|+t)^{-k}$ para todo $latex alpha, k, m > 0$

y el dato inicial las siguientes condiciones de regularidad:

$latex |partial_x^alpha u_{0i}| leq C_{alpha,k}(1+|x|)^{-k}$ para todo $latex alpha, k > 0$

Las soluciones $latex (u,p)$ admisibles para $latex x in mathbb{R}^3$ estan o en $latex C^{infty}(mathbb{R}^3 times [0,infty))$ con decaimiento en el infinito de la presión y de energía finita, es decir, $latex int_{mathbb{R}^3} |u|^2 dx < infty$ para todo $latex t$, o estan en $latex C^{infty}(mathbb{T}^3 times [0,infty))$ periódicas y la presión de media cero.

Sea $latex u_0$ satisfaciendo las condiciones de regularidad. El problema del Instituto Clay consiste, entonces, en determinar si siempre existen soluciones admisibles para $latex u_0$ o si existe algún caso en que no.

Tags: , ,

Eclipse es un IDE multiplataforma Open Source desarrollado por la Eclipse Foundation. Utiliza módulos, plug-in, para proporcionar diferentes funcionalidades, ya que se basa en la filosofía RCP (Rich Client Platform), de manera que tenemos una plataforma ligera en la que el usuario va añadiendo lo que realmente necesita.

El plug-in CDT (C/C++ Development Toolkit) permite desarrollar en Eclipse proyectos basados en estos lenguajes de programación.

Otro plug-in que encontramos es el MDT (Model Development Tools), el cual permite trabajar, entre otras cosas, con UML2 con el que especificar y diseñar software basado en el paradigma de la orientación a objetos.

Tags: , , , ,

En [Monaghan 1992] comenta el caso del método SPH en relatividad especial.

Para empezar asumimos que el fluido está constituido por bariones, por lo que el tensor de energia-momento es:

$latex T^{mu nu} = (n m_0 c^2 + n tau + P) U^mu U^nu + P g^{mu nu}$

donde los indices griegos van de $latex 0$ a $latex 3$ y los coeficientes de la métrica se definen

$latex g_{00} = -1$ y $latex g_{ij} = 1$

En estas ecuaciones, $latex n$ representa la densidad de bariones, $latex P$ es la presión, $latex tau$ la energía térmica, $latex c$ la velocidad del sonido, $latex U^nu$ la 4-velocidad con $latex U_nu U^nu = -1$ y $latex m_0$ la masa.

Las ecuaciones del momento se siguen de:

$latex frac{partial}{partial x^nu} T^{i nu} = 0$

que es:

$latex frac{d}{dt} q = – frac{1}{N} nabla P$

y en forma SPH queda:

$latex frac{d}{dt}q_a = -sum_b nu_b (frac{P_a}{n_a^2} + frac{P_b}{N_b^2}) nabla_a W_{ab}$

donde $latex nu_b$ es el número de bariones asociados a la partícula $latex b$.

Y la de la energía se sigue de:

$latex frac{partial}{partial x^j} T^{0j} = 0$

que es:

$latex frac{d}{dt} epsilon = – frac{1}{N} nabla cdot (Pv)$

y en foma SPH queda:

$latex frac{d}{dt} epsilon_a = -sum_b m_b (frac{P_a v_a}{N_a^2} + frac{P_b v_b}{N_b^2}) nabla W_{ab} $

Tags: , , , , , ,

Según Monaghan en [Monaghan, 1992], uno de los padres del método, buscaban un método con el que fuera fácil trabajar y diera buenos resultados y encontraron eso y mucho mas.

Además de no  necesitar malla para calcular las derivadas, puesto que se calculan de manera analítica a partir de una fórmula de interpolación, las ecuaciones de conservación del momento y la energía pasan a ser un conjunto de ecuaciones diferenciales ordinarias fáciles de entender.

En esencia, el método SPH es un método de interpolación que permite expresar una función a partir de los valores en un conjunto desordenado de puntos a los que llamamos particulas.

La idea es que aproximamos cualquier función $latex A(r)$ de la siguiente manera:

$latex A_I(r) = int A(r’) W(r-r’,h) dr’$

donde $latex r$ es el vector posición, $latex W$ es una función de pesos o kernel (con las propiedades ya comentadas en otros posts) y $latex h$ es la longitud de suavizado. Este tipo de aproximaciones reciben el nombre de interpolación integral.

Para trabajar de manera numérica, la interpolador integral se aproxima por:

$latex A_S(r) = sum_b m_b frac{A_b}{rho_b} W(r-r_b,h)$

donde el índice del sumatorio $latex b$ identifica a cada partícula y el sumatorio es sobre todas las partículas. La partícula $latex b$ tiene masa $latex m_b$, posición $latex r_b$, densidad $latex rho_b$ y velocidad $latex v_b$. El valor de la función $latex A$ en $latex r_b$ es $latex A_b$.

Tags: , ,

En el artículo [Rosswog 2009], Stephan Rosswog hace un repaso detallado del método SPH centrandose especialmente en sus aplicaciones en astrofísica. Repasamos el apartado que hace referencia a las ecuaciones de la hidrodinámica en forma Lagrangiana.

A diferencia de los metodos basados en malla, que son Eulerianos, es decir, métodos donde  describimos el fluido desde un punto fijo del espacio, el Smoothed Particle Hydrodynamics es totalmente Lagrangiano, por lo que describimos el fluido desde un sistema de coordenadas fijado en una particula del fluido en movimiento.

La derivada Lagrangiana o sustancial respecto del tiempo, $latex frac{d}{dt}$, se relaciona con la derivada Euleriana respecto al tiempo, $latex frac{partial}{partial t}$ de la siguiente manera:

$latex frac{d}{dt} = frac{dx^i}{dt} frac{partial}{partial x^i} + frac{partial}{partial t} = vec{v} cdot nabla + frac{partial}{partial t}$

Aplicando esta relación a las ecuaciones en forma Euleriana, las ecuaciones de la hidrodinámica en forma Lagrangiana quedan:

  1. Ecuación de continuidad: $latex frac{d}{dt} rho = – rho nabla cdot vec{v}$
  2. Ecuacion del momento: $latex frac{d}{dt} vec{v} = -frac{nabla P}{rho} + vec{f}$
  3. Ecuación de la energía: $latex frac{d}{dt}u = frac{P}{rho^2} frac{d}{dt} rho = – frac{P}{rho} nabla cdot vec{v}$
  4. Ecuación de estado, que describe la termodinámica del fluido estelar: $latex P = (gamma -1) cdot rho cdot epsilon$ (ecuación del gas ideal)

Tags: , , , , , , , , , , ,

La hipótesis de Rieman dice que todos los ceros no triviales de la función zeta de Riemann:

$latex zeta(z) = sum_{n=1}^{infty} frac{1}{n^z}$

tienen parte real $latex frac{1}{2}$.

En el entretenido libro “La música de los números primos” de Marcus de Sautoy, éste comenta la posibilidad de que la teoría de números y la física estén mas estrechamente relacionadas de lo que pensamos:

Los físicos creen que la razón por la que los ceros de Riemann deben situarse todos sobre la recta es que terminarán por ser las frecuencias de un tambor matemático. A un cero que se situara fuera de la recta le correspondería una frecuencia imaginaria prohibida por la teoría.

Y para afianzar esta idea, comenta un problema clásico de hidrodinámica resuelto por Bernhard Riemann argumentando de la misma manera:

El problema se refiere a una esfera de fluido en rotación que se mantiene unida gracias a interacciones gravitacionales recíprocas entre las partículas que la componen. Una estrella, por ejemplo, es una enorme bola de gas giratorio que se mantiene unido por su propia gravedad. La cuestión es: ¿qué sucederá con la bola si se le da una patada?¿Se limitará a temblar ligeramente o se desitegrará? Para responder a estas preguntas es necesario determinar si ciertos números imaginarios determinados están o no alineados. Si lo están, la esfera de fluido en rotación quedará intacta.

Los Nachlass son un conjunto de notas inéditas de Riemann que están en la biblioteca de Gotinga. Cuenta el libro que cuando el físico Jon Keating pidió dos partes de los Nachlass, una correspondiente a sus intentos de demostración de la hipotesis de Riemann y otra a sus estudios en hidrodinámica, el bibliotecario le entrego un único grupo de documentos. De ahí que Marcus escriba:

Una vez más, los Nachlass revelaban hasta qué punto Riemann se adelantó a su tiempo. Es imposible que no fuera consciente del significado que implicaba su solución al problema de dinámica de fluidos. Su método había demostrado por qué cietos números imaginarios que aparecían en su análisis de la esfera de fluido se colocaban en línea recta; y al mismo tiempo -y en los mismos folios- estaba intentando demostrar por qué los ceros de su paisaje zeta se situaban todos sobre la misma línea.

 

Tags: , , , ,

Las heramientas CASE (Computer Aided Software Engineering) pretenden aplicar los métodos de la ingenieria a la producción de software y surgen en respuesta a la crisis del software de finales de los sesenta. Cubren todas las etapas del ciclo de vida del software y básicamente pretende generar software de calidad minimizando tiempo y costes tanto en el desarrollo como en el mantenimiento del mismo.

El lenguaje UML (Unified Modeling Language) es un estandar que permite especificar y diseñar software basado en el paradigma de la orientación a objetos.

En este enlace podemos ver  diferentes herramientas CASE que soportan UML y que permiten  generar código automáticamente en diferentes lenguajes de programación.

Tags: ,

Existen diferentes posibilidades a la hora de definir una función kernel:

  1. Gaussiana [Gingold & Monaghan, 1977]:$latex W(r,h) = alpha_D cdot e^{-q^2}$ con $latex 0 leq q leq 2$ donde $latex q=frac{r}{h}$, $latex r$ es la distancia entre dos partículas determinadas y $latex alpha_D$, el factor dimensional, que es $latex frac{1}{pi h^2}$ en dos dimensiones y $latex frac{1}{pi^{frac{3}{2}} h^3}$ en tres.
  2. Cuadrática [Gingold & Monaghan, 1977]: $latex alpha_D (frac{3}{16} q^2 – frac{3}{4} q + frac{3}{4})$ con $latex 0 leq q leq 2$ y donde $latex alpha_D$ es $latex frac{2}{pi h^2}$ en 2D y $latex frac{5}{4 pi h^3}$ en 3D.
  3. B-spline cúbico [Monaghan & Lattancio, 1985]: $latex W(r,h) = alpha_D begin{cases} 1 – frac{3}{2} q^2 + frac{3}{4} q^3 ,& mbox{if } 0 leq q leq 1 \ frac{1}{4}(2-q)^3 ,& mbox{if } 1 < q leq 2 \ 0 ,& mbox{if } q > 2end{cases}$ con $latex alpha_D = frac{10}{7 pi h^2}$ en dos dimensiones y $latex frac{1}{pi h^3}$ en 3 dimensiones.
  4. Quíntica: $latex W(r,h) = alpha_D (1-frac{q}{2}^4)(2q + 1)$ con $latex 0 leq q leq 2$ y $latex alpha_D = frac{7}{4 pi h^2}$ en 2D y $latex frac{7}{8 pi h^3}$ en 3D.

Pensando en posibles implementaciones, parece lógico utilizar una classe abstracta Kernel de la que heredaran las anteriores implementando sus métodos, de manera que, las posibles clases que llamen a la misma puedan hacerlo siempre de la misma manera independientement de cual estemos utilizando.

Tags: , , , ,

El siguiente código muestra la definición de una clase en el lenguaje C++. En el podemos ver los elementos típicos de las mismas.

En el fichero Complex.h, tendriamos la representación del tipo:

[sourcecode languaje=”cpp”]

#ifndef ComplexH
#define ComplexH

class Complex {

private:

//atributos
double re;
double im;

public:

//constructores
Complex();
Complex(double re, double im);

//metodos
double get_re();
double get_im();
void set_re(double re);
void set_im(double im);
double get_mod();
double get_arg();

//destructor
~Complex();

};
[/sourcecode]

Mientras que en el fichero Complex.cpp tendriamos la implementación del mismo:

[sourcecode languaje=”cpp”]
#include "Complex.h"

//constructores
Complex::Complex() {
re = 0;
im = 0;
}

Complex::Complex(double re, double im) {
this.re = re;
this.im = im;
}

//destructor
Complex::~Complex() {
}

//metodos
double Complex::get_re() {
return this.re;
}

double Complex::get_im() {
return this.im;
}

void Complex::set_re(double re) {
this.re = re;
}

void Complex::set_im(double im) {
this.im = im;
}

double Complex::get_mod() {
return …;
}

double Complex::get_arg(){
return …;
}

};
[/sourcecode]

« Older entries § Newer entries »

FireStats icon Powered by FireStats