noviembre 2012

You are currently browsing the monthly archive for noviembre 2012.

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

FireStats icon Powered by FireStats