You are currently browsing the monthly archive for noviembre 2012.

Ya comentamos en el post anterior que formato tienen los sistemas 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 n dimensiones solo significa colocar nuestra expresión en n bucles anidados. Así de sencillo…

En el caso de 1D, tendriamos:

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] );

}

En 1D tenemos:

u_{xx} = f(x),

con 0 < x < 1 y u(0)= u(1)=0. Particionamos el dominio en n subintervalos introduciendo los puntos x_i = ih con 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 x_i):

\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 \xi_i \in (x_{i-1}, x_{i+1}) por lo que:

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

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

\frac{1}{h^2} \begin{bmatrix} -2 & 1 & & & & & \\ 1 & -2 & 1 & & & & \\ & 1 & -2 & 1 & & & \\ & & \cdot & \cdot & \cdot & & \\ & & & \cdot & \cdot & \cdot & \\ & & & & 1 & -2 & 1 \\ & & & & & 1 & -2\end{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:

\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 \xi_i \in (x_{i-2},x_{i+2}).

En 2D tenemos:

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

con 0 < x < 1, 0 < y < 1 y u=0 en el cuadrado unidad. Particionamos el dominio en n_x \times n_y subintervalos introduciendo los puntos (x_i,y_j) = (ih_x, jh_y) con h_x=\frac{1}{n_x} y 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:

\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 v_{0\,j} = v_{n_x\,j} = v_{i\,0} = u_{i\,n_y} = 0. Tenemos un error local de truncamiento del orden de O(h_x^2+h_y^2).

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

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

En 3D tenemos:

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

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

\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:

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

si h_x = h_y = h_z = h.

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:

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

para encontrar

\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:

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

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

(3) Resolver la ecuación:

\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 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:

\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 \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 \beta anterior, por ejemplo, en coordenadas esféricas, tendríamos:

(1) Resolver ecuación:

\Delta \mu = \mu_S

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

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

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

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

(i) \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) \Delta \chi = r S(\boldsymbol{\beta})^r - \frac{r}{3} \frac{\partial \Theta}{\partial r} + 2 \Theta , donde \chi = r \beta^r

(4) Deducir \eta de una de las siguientes ecuaciones:

(i) \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 -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 \beta^r que puede tener problemas con el orden.

(ii) \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.

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

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

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

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

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

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

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

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

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

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

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

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

para la primera y:

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

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

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

para la segunda.

De esta manera, sustituyendo, tenemos:

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

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

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

\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:

\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:

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}),

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

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

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

Vamos a suponer que estamos en dimensión n y miraremos como quedan los diferentes kernels y sus derivadas. Denotaremos con \boldsymbol{q} a (q_1,\ldots,q_n) y 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 n es una n-tupla de enteros no negativos \alpha = (\alpha_1,\ldots,\alpha_n) con \alpha_j \in \mathbb{Z}_+ para j= 1\ldots n, de manera que \alpha \in \mathbb{Z}_+^n. Llamamos  módulo de \alpha a:

|\alpha| = \sum_{i=1}^n \alpha_i

Además, podemos definir

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

y

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

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

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

kernel Gaussiano:

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

\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

\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

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:

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 q = \sqrt{ \sum_{i=1}^n q_i^2 }, por lo que \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:

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

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

\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 0 \leq |\alpha| \leq 3 (si |\alpha|>3 entonces W \equiv 0) podemos escribir:

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

=\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:

 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 |\alpha| = 0 \ldots 4:

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

=\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:

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 |\alpha| = 0 \ldots 5:

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

=\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}

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

\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 a situada en r_a = (x_a,y_a) tenemos allí el kernel W_a(r-r_a,h). Dada otra partícula vecina b de posición r_b = (x_b, y_b), el punto

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

queda sobre el kernel:

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 q de todos los objetos del sistema, sino también sus momentos p. Matemáticamente, el espacio de configuraciones puede definirse como una variedad M de manera que, para cada posicion q \in M, los momentos p toman valores en el espacio  cotangente T_q^*M. De esta manera, el espacio de fases puede verse de manera natural como el fibrado cotangente:

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

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

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

\frac{d\boldsymbol{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:

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

\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),

P_a and P_b

are obtained from the EoS and \frac{d\boldsymbol{r}_a}{dt} = \boldsymbol{v}_a.

In SHD we have:

d

Finally, the GHD equations are:

d

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 x^3 + bx^2 + cx+d = 0, la fórmula de las soluciones es:

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

donde:

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

Construimos una malla de 100\times 100 puntos de nuestros kernels en 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:

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

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

En 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):

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:

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

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?

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

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:

\gamma_{ij} = \psi^4 \delta_{ij}

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

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

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

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

\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 E = \rho h W^2 - P, S = \rho h (W^2 -1) + 3P y

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:

\Delta B^i = S_\beta

\Delta \chi = \partial_i B^i

si definimos \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.

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.

noviembre 2012
L M X J V S D
« Oct   Dic »
 1234
567891011
12131415161718
19202122232425
2627282930