PDE elípticas

You are currently browsing articles tagged PDE elípticas.

El software Chombo, del Berkeley Lab, combina los métodos en diferencias finitas con las mallas adaptativas (AMR) para resolver, entre otras, PDEs elípticas.

Las siguientes imágenes, en 2D y 3D, se han obtenido a partir de su AMRPoisson:

schSolSouschSol

Tags: , , ,

En la discretización que hicimos teníamos dos sistemas acoplados, uno para las $latex X^i$ y otro para las $latex beta^i$. Procedemos ahora a desacoplarlos.

Para empezar, tomamos la divergencia (plana) del sistema:

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

y, teniendo en cuenta que $latex mathcal{D}$ conmuta con $latex Delta$ (métrica plana), tenemos:

$latex Delta (mathcal{D}_i X^i) = 8 pi mathcal{D}^j S^*_j – frac{1}{3} Delta (mathcal{D}_j X^j)$,

por lo que:

$latex Delta (mathcal{D}_i X^i) = frac{3}{4} 8 pi mathcal{D}^j S^*_j$.

De esta manera, si definimos $latex Theta_X := mathcal{D}_i X^i$, nos queda:

$latex Delta Theta_X = frac{3}{4} 8 pi mathcal{D}^j S^*_j = 6 pi (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )$,

que discretizado queda:

$latex frac{(Theta_X)_{i-1,j,k}-2(Theta_X)_{i,j,k}+(Theta_X)_{i+1,j,k}}{h_x^2} + $

$latex frac{(Theta_X)_{i,j-1,k}-2(Theta_X)_{i,j,k}+(Theta_X)_{i,j+1,k}}{h_y^2} + $

$latex frac{(Theta_X)_{i,j,k-1}-2(Theta_X)_{i,j,k}+(Theta_X)_{i,j,k+1}}{h_z^2} = $

$latex = 6 pi (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )_{i,j,k} $,

donde inicialmente:

$latex (S^*_a)_{i,j,k} = (psi^6)_{i,j,k}rho_{i,j,k}h_{i,j,k}w^2_{i,j,k}(v_a)_{i,j,k}$,

$latex (partial_x S^*_x + partial_y S^*_y +partial_z S^*_z )_{i,j,k} = $

$latex frac{(S^*_x)_{i+1,j,k}-(S^*_x)_{i-1,j,k}}{2h_x} + frac{(S^*_x)_{i,j+1,k}-(S^*_x)_{i,j-1,k}}{2h_y} + frac{(S^*_x)_{i,j,k+1}-(S^*_x)_{i,j,k-1}}{2h_z}$

y que es lineal.

El primer sistema acoplado de ecuaciones quedaría ahora:

$latex partial_{xx} X^x + partial_{yy} X^x + partial_{zz} X^x = 8 pi S^*_x – frac{1}{3} partial_x Theta_X approx $

$latex approx frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} + frac{X^x_{i,j-1,k}-2X^x_{i,j,k}+X^x_{i,j+1,k}}{h_y^2} + frac{X^x_{i,j,k-1}-2X^x_{i,j,k}+X^x_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_x)_{i,j,k} – frac{1}{3} (partial_x Theta_X)_{i,j,k}$,

¡que vuelve a ser lineal!

Continuamos con:

$latex partial_{xx} X^y + partial_{yy} X^y + partial_{zz} X^y = 8 pi S^*_y – frac{1}{3} partial_y Theta_X approx $

$latex approx frac{X^y_{i-1,j,k}-2X^y_{i,j,k}+X^y_{i+1,j,k}}{h_x^2} + frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} + frac{X^y_{i,j,k-1}-2X^y_{i,j,k}+X^y_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_y)_{i,j,k} – frac{1}{3} (partial_y Theta_X)_{i,j,k}$

y, finalmente:

$latex partial_{xx} X^z + partial_{yy} X^z + partial_{zz} X^z = 8 pi S^*_z – frac{1}{3} partial_z Theta_X approx $

$latex approx frac{X^z_{i-1,j,k}-2X^z_{i,j,k}+X^z_{i+1,j,k}}{h_x^2} + frac{X^z_{i,j-1,k}-2X^z_{i,j,k}+X^z_{i,j+1,k}}{h_y^2} + frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi (S^*_z)_{i,j,k} – frac{1}{3} (partial_z Theta_X)_{i,j,k}$,

donde calculamos al principio:

$latex (partial_x Theta_X)_{i,j,k} = frac{(Theta_X)_{i+1,j,k}-(Theta_X)_{i-1,j,k}}{2h_x}$

$latex (partial_y Theta_X)_{i,j,k} = frac{(Theta_X)_{i,j+1,k}-(Theta_X)_{i,j-1,k}}{2h_y}$

$latex (partial_z Theta_X)_{i,j,k} = frac{(Theta_X)_{i,j,k+1} – (Theta_X)_{i,j,k-1}}{2h_z}$

A continuación, discretizamos las siguientes ecuaciones:

$latex hat{A}^{xx} = 2 partial_x X^x – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx$

$latex approx frac{2}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{xx}_{i,j,k}$,

$latex hat{A}^{xy} = hat{A}^{yx}= partial_x X^y + partial_y X^x approx $

$latex approx frac{X^y_{i+1,j,k}-X^y_{i-1,j,k}}{2h_x} + frac{X^x_{i,j+1,k}-X^x_{i,j-1,k}}{2h_y} = hat{A}^{xy}_{i,j,k} = hat{A}^{yx}_{i,j,k}$,

$latex hat{A}^{xz} = hat{A}^{zx} = partial_x X^z + partial_z X^x approx $

$latex approx frac{X^z_{i+1,j,k}-X^z_{i-1,j,k}}{2h_x} + frac{X^x_{i,j,k+1}-X^x_{i,j,k-1}}{2h_z} = hat{A}^{xz}_{i,j,k} = hat{A}^{zx}_{i,j,k}$,

$latex hat{A}^{yy} = 2 partial_y X^y – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} +frac{2}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{yy}_{i,j,k}$,

$latex hat{A}^{yz} = hat{A}^{zy} = partial_y X^z + partial_z X^y approx $

$latex approx frac{X^z_{i,j+1,k}-X^z_{i,j-1,k}}{2h_y} + frac{X^y_{i,j,k+1}-X^y_{i,j,k-1}}{2h_z} = hat{A}^{yz}_{i,j,k} = hat{A}^{zy}_{i,j,k}$,

$latex hat{A}^{zz} = 2 partial_z X^z – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} + frac{2}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{zz}_{i,j,k}$.

Por tanto, la siguiente ecuación:

$latex Delta psi = -2 pi psi^{-1} E^* – psi^{-7} frac{(hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2}{8}$

queda:

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

$latex =-2 pi psi^{-1}_{i,j,k} E^*_{i,j,k} – $

$latex – frac{psi^{-7}_{i,j,k}}{8} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 ) $,

con:

$latex partial_{psi_{i,j,k}} F(psi_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) -2 pi psi_{i,j,k}^{-2} E^*_{i,j,k} – $

$latex – frac{7}{8} psi^{-8}_{i,j,k} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 )$,

donde:

$latex E^*_{i,j,k} = psi^{6}_{i,j,k} (D_{i,j,k}+tau_{i,j,k})$

y la ecuación:

$latex Delta (alphapsi) = (alpha psi) (2 pi psi^{-2} (E^*+2S^*) + $

$latex + frac{7}{8} psi^{-8} ((hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2) )$

como:

$latex approx frac{(alphapsi)_{i-1,j,k} – 2(alphapsi)_{i,j,k}+(alphapsi)_{i+1,j,k}}{h_x^2} + $

$latex + frac{(alphapsi)_{i,j-1,k}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j+1,k}}{h_y^2} + $

$latex + frac{(alphapsi)_{i,j,k-1}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j,k+1}}{h_z^2} = $

$latex = (alpha psi)_{i,j,k} (2 pi psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) + $

$latex + frac{7}{8} psi^{-8}_{i,j,k} ((hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2) )$,

donde:

$latex partial_{(alpha psi)_{i,j,k}} F((alpha psi)_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) – 2 pi psi^{-2}_{i,j,k} (E^*_{i,j,k}+2S^*_{i,j,k}) + $

$latex – frac{7}{8} psi^{-8}_{i,j,k} ((hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2) )$

con:

$latex S^*_{i,j,k} = psi^6_{i,j,k}(rho_{i,j,k}h_{i,j,k}(w^2_{i,j,k}-1) + 3 p_{i,j,k})$.

Finalmente, tenemos el otro sistema acoplado:

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

con el que procedemos de igual manera que con las $latex X^i$:

$latex Delta(mathcal{D}_i beta^i) = mathcal{D}_i (mathcal{D}_j (2 alpha psi^{-6} hat{A}^{ij})) – frac{1}{3} Delta (mathcal{D}_i beta^i)$,

de manera que:

$latex Delta Theta_beta = frac{3}{4} mathcal{D}^i (mathcal{D}_j (2 alpha psi^{-6} hat{A}^{ij})) =$

$latex frac{3}{2}(partial_{xx}(alpha psi^{-6} hat{A}^{xx}) + partial_{yy}(alpha psi^{-6} hat{A}^{yy}) + partial_{zz}(alpha psi^{-6} hat{A}^{zz})$,

con:

$latex Theta_beta := mathcal{D}_i beta^i$,

que discretizada queda:

$latex frac{(Theta_beta)_{i-1,j,k}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i+1,j,k}}{h_x^2} + $

$latex frac{(Theta_beta)_{i,j-1,k}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i,j+1,k}}{h_y^2} + $

$latex frac{(Theta_beta)_{i,j,k-1}-2(Theta_beta)_{i,j,k}+(Theta_beta)_{i,j,k+1}}{h_z^2} = $

$latex frac{3}{2}((partial_{xx}(alpha psi^{-6} hat{A}^{xx}))_{i,j,k} + (partial_{yy}(alpha psi^{-6} hat{A}^{yy}))_{i,j,k} + (partial_{zz}(alpha psi^{-6} hat{A}^{zz})_{i,j,k})$,

De esta manera, tenemos:

$latex Delta beta^x = partial_x (2 alpha psi^{-6} hat{A}^{xx}) + partial_y (2 alpha psi^{-6} hat{A}^{xy}) + partial_z (2 alpha psi^{-6} hat{A}^{xz}) – frac{1}{2} partial_x Theta_beta approx $

$latex approx frac{beta^x_{i-1,j,k}-2beta^x_{i,j,k}+beta^x_{i+1,j,k}}{h_x^2} + frac{beta^x_{i,j-1,k}-2beta^x_{i,j,k}+beta^x_{i,j+1,k}}{h_y^2} + frac{beta^x_{i,j,k-1}-2beta^x_{i,j,k}+beta^x_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{xx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{xy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{xz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_x Theta_beta)_{i,j,k}$.

De la misma manera:

$latex Delta beta^y = partial_x (2 alpha psi^{-6} hat{A}^{yx}) + partial_y (2 alpha psi^{-6} hat{A}^{yy}) + partial_z (2 alpha psi^{-6} hat{A}^{yz}) – frac{1}{3} partial_y Theta_beta approx $

$latex approx frac{beta^y_{i-1,j,k}-2beta^y_{i,j,k}+beta^y_{i+1,j,k}}{h_x^2} + frac{beta^y_{i,j-1,k}-2beta^y_{i,j,k}+beta^y_{i,j+1,k}}{h_y^2} + frac{beta^y_{i,j,k-1}-2beta^y_{i,j,k}+beta^y_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{yx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{yy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{yz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_y Theta_beta)_{i,j,k}$.

Y, por último:

$latex Delta beta^z = partial_x (2 alpha psi^{-6} hat{A}^{zx}) + partial_y (2 alpha psi^{-6} hat{A}^{zy}) + partial_z (2 alpha psi^{-6} hat{A}^{zz}) – frac{1}{3} partial_z Theta_beta approx $

$latex approx frac{beta^z_{i-1,j,k}-2beta^z_{i,j,k}+beta^z_{i+1,j,k}}{h_x^2} + frac{beta^z_{i,j-1,k}-2beta^z_{i,j,k}+beta^z_{i,j+1,k}}{h_y^2} + frac{beta^z_{i,j,k-1}-2beta^z_{i,j,k}+beta^z_{i,j,k+1}}{h_z^2} = $

$latex = (partial_x (2 alpha psi^{-6} hat{A}^{zx}))_{i,j,k} + (partial_y (2 alpha psi^{-6} hat{A}^{zy}))_{i,j,k} + (partial_z (2 alpha psi^{-6} hat{A}^{zz}) )_{i,j,k} – $

$latex – frac{1}{3} (partial_z Theta_beta)_{i,j,k}$.

Parece que, del sistema no lineal acoplado inicial, hemos llegado a un sistema de diez ecuaciones desacopladas donde ocho de ellas son lineales y solo dos son no linales. No pinta mal. Ya escribiremos próximamente sobre las condiciones de contorno…

Tags: , , , ,

Vamos a discretizar las ecuaciones que comentamos en este post. Para ello, discretizaremos las derivadas de la siguiente manera:

$latex partial_x u = frac{u_{i+1,j,k}-u_{i-1,j,k}}{2h_x}$,

$latex partial_y u = frac{u_{i,j+1,k}-u_{i,j-1,k}}{2h_y}$,

$latex partial_z u = frac{u_{i,j,k+1}-u_{i,j,k-1}}{2h_z}$,

$latex partial_{xx} u = frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{h_x^2}$,

$latex partial_{yy} u = frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{h_y^2}$,

$latex partial_{zz} u = frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{h_z^2}$,

$latex partial_{xy} u = frac{u_{i-1,j-1,k}-u_{i+1,j-1,k}-u_{i-1,j+1,k}+u_{i+1,j+1,k}}{4h_xh_y}$,

$latex partial_{xz} u = frac{u_{i-1,j,k-1}-u_{i+1,j,k-1}-u_{i-1,j,k+1}+u_{i+1,j,k+1}}{4h_xh_z}$,

$latex partial_{yz} u = frac{u_{i,j-1,k-1}-u_{i,j+1,k-1}-u_{i,j-1,k+1}+u_{i,j+1,k+1}}{4h_yh_z}$.

El primer grupo de ecuaciones quedaría:

$latex partial_{xx} X^x + partial_{yy} X^x + partial_{zz} X^x = 8 pi psi^6 rho h w^2 v_x – frac{1}{3} partial_x (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} + frac{X^x_{i,j-1,k}-2X^x_{i,j,k}+X^x_{i,j+1,k}}{h_y^2} + frac{X^x_{i,j,k-1}-2X^x_{i,j,k}+X^x_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi psi^6_{i,j,k} rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{x_{i,j,k}} – frac{1}{3} ( frac{X^x_{i-1,j,k}-2X^x_{i,j,k}+X^x_{i+1,j,k}}{h_x^2} +$

$latex + frac{X^y_{i-1,j-1,k}-X^y_{i+1,j-1,k}-X^y_{i-1,j+1,k}+X^y_{i+1,j+1,k}}{4h_xh_y} +$

$latex + frac{X^z_{i-1,j,k-1}-X^z_{i+1,j,k-1}-X^z_{i-1,j,k+1}+X^z_{i+1,j,k+1}}{4h_xh_z} )$,

y además, para los esquemas de relajación no lineales, reescribimos la igualdad anterior como $latex F(X^x_{i,j,k})=0$ y entonces tenemos:

$latex partial_{X^x_{i,j,k}} F(X^x_{i,j,k}) = -2 ( frac{4}{3}frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2})$.

$latex partial_{xx} X^y + partial_{yy} X^y + partial_{zz} X^y = 8 pi psi^6 rho h w^2 v_y – frac{1}{3} partial_y (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx frac{X^y_{i-1,j,k}-2X^y_{i,j,k}+X^y_{i+1,j,k}}{h_x^2} + frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} + frac{X^y_{i,j,k-1}-2X^y_{i,j,k}+X^y_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi psi^6_{i,j,k} rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{y_{i,j,k}} – frac{1}{3} ( frac{X^x_{i-1,j-1,k}-X^x_{i+1,j-1,k}-X^x_{i-1,j+1,k}+X^x_{i+1,j+1,k}}{4h_xh_y} +$

$latex + frac{X^y_{i,j-1,k}-2X^y_{i,j,k}+X^y_{i,j+1,k}}{h_y^2} +$

$latex + frac{X^z_{i-1,j,k-1}-X^z_{i+1,j,k-1}-X^z_{i-1,j,k+1}+X^z_{i+1,j,k+1}}{4h_yh_z} )$,

con:

$latex partial_{X^y_{i,j,k}} F(X^y_{i,j,k}) = -2 ( frac{1}{h_x^2} +frac{4}{3} frac{1}{h_y^2} + frac{1}{h_z^2})$.

$latex partial_{xx} X^z + partial_{yy} X^z + partial_{zz} X^z = 8 pi psi^6 rho h w^2 v_z – frac{1}{3} partial_z (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx frac{X^z_{i-1,j,k}-2X^z_{i,j,k}+X^z_{i+1,j,k}}{h_x^2} + frac{X^z_{i,j-1,k}-2X^z_{i,j,k}+X^z_{i,j+1,k}}{h_y^2} + frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2} = $

$latex = 8 pi psi^6_{i,j,k} rho_{i,j,k} h_{i,j,k} w^2_{i,j,k} v_{z_{i,j,k}} – frac{1}{3} ( frac{X^x_{i-1,j,k-1}-X^x_{i+1,j,k-1}-X^x_{i-1,j,k+1}+X^x_{i+1,j,k+1}}{4h_xh_z} +$

$latex + frac{X^y_{i,j-1,k-1}-X^y_{i,j+1,k-1}-X^y_{i,j-1,k+1}+X^y_{i,j+1,k+1}}{4h_yh_z} )$

$latex + frac{X^z_{i,j,k-1}-2X^z_{i,j,k}+X^z_{i,j,k+1}}{h_z^2}$

con:

$latex partial_{X^z_{i,j,k}} = F(X^z_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{4}{3} frac{1}{h_z^2})$.

A continuación, discretizamos las siguientes ecuaciones:

$latex hat{A}^{xx} = 2 partial_x X^x – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx$

$latex approx frac{2}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{xx}_{i,j,k}$,

$latex hat{A}^{xy} = hat{A}^{yx}= partial_x X^y + partial_y X^x approx $

$latex approx frac{X^y_{i+1,j,k}-X^y_{i-1,j,k}}{2h_x} + frac{X^x_{i,j+1,k}-X^x_{i,j-1,k}}{2h_y} = hat{A}^{xy}_{i,j,k} = hat{A}^{yx}_{i,j,k}$,

$latex hat{A}^{xz} = hat{A}^{zx} = partial_x X^z + partial_z X^x approx $

$latex approx frac{X^z_{i+1,j,k}-X^z_{i-1,j,k}}{2h_x} + frac{X^x_{i,j,k+1}-X^x_{i,j,k-1}}{2h_z} = hat{A}^{xz}_{i,j,k} = hat{A}^{zx}_{i,j,k}$,

$latex hat{A}^{yy} = 2 partial_y X^y – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} +frac{2}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} – frac{1}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{yy}_{i,j,k}$,

$latex hat{A}^{yz} = hat{A}^{zy} = partial_y X^z + partial_z X^y approx $

$latex approx frac{X^z_{i,j+1,k}-X^z_{i,j-1,k}}{2h_y} + frac{X^y_{i,j,k+1}-X^y_{i,j,k-1}}{2h_z} = hat{A}^{yz}_{i,j,k} = hat{A}^{zy}_{i,j,k}$,

$latex hat{A}^{zz} = 2 partial_z X^z – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z) approx $

$latex approx -frac{1}{3}frac{X^x_{i+1,j,k}-X^x_{i-1,j,k}}{h_x} -frac{1}{3} frac{X^y_{i,j+1,k}-X^y_{i,j-1,k}}{h_y} + frac{2}{3} frac{X^z_{i,j,k+1}-X^z_{i,j,k-1}}{h_z}) = hat{A}^{zz}_{i,j,k}$.

Por tanto, la siguiente ecuación:

$latex Delta psi = -2 pi psi^{-1} (D + tau) – psi^{-7} frac{(hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2}{8}$

queda:

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

$latex =-2 pi psi^{-1}_{i,j,k} (D_{i,j,k}+tau_{i,j,k}) – $

$latex – frac{psi^{-7}_{i,j,k}}{8} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 ) $,

con:

$latex partial_{psi_{i,j,k}} F(psi_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) -2 pi psi_{i,j,k}^{-2} (D_{i,j,k}+tau_{i,j,k}) – $

$latex – frac{7}{8} psi^{-8}_{i,j,k} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 )$.

y la ecuación:

$latex Delta (alphapsi) = 2 pi (alphapsi)^{-1} ( D + tau + 2 rho h (w^2-1) + 6 p) + $

$latex + frac{7}{8} (alpha psi)^{-7} ((hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2)$

como:

$latex approx frac{(alphapsi)_{i-1,j,k} – 2(alphapsi)_{i,j,k}+(alphapsi)_{i+1,j,k}}{h_x^2} + frac{(alphapsi)_{i,j-1,k}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j+1,k}}{h_y^2} + frac{(alphapsi)_{i,j,k-1}-2(alphapsi)_{i,j,k}+(alphapsi)_{i,j,k+1}}{h_z^2} = $

$latex =2 pi (alphapsi)_{i,j,k}^{-1} (D_{i,j,k}+tau_{i,j,k} + 2 rho_{i,j,k} h_{i,j,k} (w^2_{i,j,k}-1)+6p_{i,j,k}) + $

$latex + frac{7}{8}(alphapsi)_{i,j,k}^{-7} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 ) $,

donde:

$latex partial_{psialpha_{i,j,k}} F(psialpha_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2} ) + $

$latex + 2 pi (psialpha)_{i,j,k}^{-2} (D_{i,j,k}+tau_{i,j,k} + 2 rho_{i,j,k} h_{i,j,k} (w^2_{i,j,k}-1)+6p_{i,j,k}) – $

$latex + frac{49}{8} (psialpha)_{i,j,k}^{-8} ( (hat{A}^{xx}_{i,j,k})^2+(hat{A}^{yy}_{i,j,k})^2+(hat{A}^{zz}_{i,j,k})^2+2(hat{A}^{xy}_{i,j,k})^2+2(hat{A}^{xz}_{i,j,k})^2+2(hat{A}^{yz}_{i,j,k})^2 )$.

Finalmente, tenemos:

$latex Delta beta^x = partial_x (2 alpha psi^{-6} hat{A}^{xx}) + partial_y (2 alpha psi^{-6} hat{A}^{xy}) + partial_z (2 alpha psi^{-6} hat{A}^{xz}) – $

$latex – frac{1}{3} partial_x (partial_x beta^x + partial_y beta^y + partial_z beta^z) approx $

$latex approx frac{beta^x_{i-1,j,k}-2beta^x_{i,j,k}+beta^x_{i+1,j,k}}{h_x^2} + frac{beta^x_{i,j-1,k}-2beta^x_{i,j,k}+beta^x_{i,j+1,k}}{h_y^2} + frac{beta^x_{i,j,k-1}-2beta^x_{i,j,k}+beta^x_{i,j,k+1}}{h_z^2} = $

$latex = frac{(alpha psi)_{i+1,j,k}^{-6} hat{A}_{i+1,j,k}^{xx} – (alpha psi)_{i-1,j,k}^{-6} hat{A}_{i-1,j,k}^{xx}}{h_x} + $

$latex + frac{(alpha psi)_{i,j+1,k}^{-6} hat{A}_{i,j+1,k}^{xy} – (alpha psi)_{i,j-1,k}^{-6} hat{A}_{i,j-1,k}^{xy}}{h_y} + $

$latex + frac{(alpha psi)_{i,j,k+1}^{-6} hat{A}_{i,j,k+1}^{xz} – (alpha psi)_{i,j,k-1}^{-6} hat{A}_{i,j,k-1}^{xz}}{h_z} – $

$latex – frac{1}{3} ( frac{beta^x_{i-1,j,k}-2beta^x_{i,j,k}+beta^x_{i+1,j,k}}{h_x^2} + $

$latex + frac{beta^y_{i-1,j-1,k}-beta^y_{i+1,j-1,k}-beta^y_{i-1,j+1,k}+beta^y_{i+1,j+1,k}}{4 h_x h_y} + $

$latex + frac{beta^z_{i-1,j,k-1}-beta^z_{i+1,j,k-1}-beta^z_{i-1,j,k+1}+beta^z_{i+1,j,k+1}}{4 h_x h_z} $,

con:

$latex partial_{beta^x_{i,j,k}} F(beta^x_{i,j,k}) = -2 ( frac{4}{3}frac{1}{h_x^2} + frac{1}{h_y^2} + frac{1}{h_z^2})$,

$latex Delta beta^y = partial_x (2 alpha psi^{-6} hat{A}^{yx}) + partial_y (2 alpha psi^{-6} hat{A}^{yy}) + partial_z (2 alpha psi^{-6} hat{A}^{yz}) – $

$latex – frac{1}{3} partial_y (partial_x beta^x + partial_y beta^y + partial_z beta^z) approx $

$latex approx frac{beta^y_{i-1,j,k}-2beta^y_{i,j,k}+beta^y_{i+1,j,k}}{h_x^2} + frac{beta^y_{i,j-1,k}-2beta^y_{i,j,k}+beta^y_{i,j+1,k}}{h_y^2} + frac{beta^y_{i,j,k-1}-2beta^y_{i,j,k}+beta^y_{i,j,k+1}}{h_z^2} = $

$latex = frac{(alpha psi)_{i+1,j,k}^{-6} hat{A}_{i+1,j,k}^{yx} – (alpha psi)_{i-1,j,k}^{-6} hat{A}_{i-1,j,k}^{yx}}{h_x} + $

$latex + frac{(alpha psi)_{i,j+1,k}^{-6} hat{A}_{i,j+1,k}^{yy} – (alpha psi)_{i,j-1,k}^{-6} hat{A}_{i,j-1,k}^{yy}}{h_y} + $

$latex + frac{(alpha psi)_{i,j,k+1}^{-6} hat{A}_{i,j,k+1}^{yz} – (alpha psi)_{i,j,k-1}^{-6} hat{A}_{i,j,k-1}^{yz}}{h_z} – $

$latex – frac{1}{3} ( frac{beta^x_{i-1,j-1,k}-beta^x_{i+1,j-1,k}-beta^x_{i-1,j+1,k}+beta^x_{i+1,j+1,k}}{4h_xh_y} +$

$latex + frac{beta^y_{i,j-1,k}-2beta^y_{i,j,k}+beta^y_{i,j+1,k}}{h_y^2} +$

$latex + frac{beta^z_{i-1,j,k-1}-beta^z_{i+1,j,k-1}-beta^z_{i-1,j,k+1}+beta^z_{i+1,j,k+1}}{4h_yh_z} )$,

con:

$latex partial_{beta^y_{i,j,k}} F(beta^y_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{4}{3} frac{1}{h_y^2} + frac{1}{h_z^2})$,

$latex Delta beta^z = partial_x (2 alpha psi^{-6} hat{A}^{zx}) + partial_y (2 alpha psi^{-6} hat{A}^{zy}) + partial_z (2 alpha psi^{-6} hat{A}^{zz}) – $

$latex – frac{1}{3} partial_z (partial_x beta^x + partial_y beta^y + partial_z beta^z) approx $

$latex approx frac{beta^z_{i-1,j,k}-2beta^z_{i,j,k}+beta^z_{i+1,j,k}}{h_x^2} + frac{beta^z_{i,j-1,k}-2beta^z_{i,j,k}+beta^z_{i,j+1,k}}{h_y^2} + frac{beta^z_{i,j,k-1}-2beta^z_{i,j,k}+beta^z_{i,j,k+1}}{h_z^2} = $

$latex = frac{(alpha psi)_{i+1,j,k}^{-6} hat{A}_{i+1,j,k}^{zx} – (alpha psi)_{i-1,j,k}^{-6} hat{A}_{i-1,j,k}^{zx}}{h_x} + $

$latex + frac{(alpha psi)_{i,j+1,k}^{-6} hat{A}_{i,j+1,k}^{zy} – (alpha psi)_{i,j-1,k}^{-6} hat{A}_{i,j-1,k}^{zy}}{h_y} + $

$latex + frac{(alpha psi)_{i,j,k+1}^{-6} hat{A}_{i,j,k+1}^{zz} – (alpha psi)_{i,j,k-1}^{-6} hat{A}_{i,j,k-1}^{zz}}{h_z} – $

$latex – frac{1}{3} ( frac{beta^x_{i-1,j,k-1}-beta^x_{i+1,j,k-1}-beta^x_{i-1,j,k+1}+beta^x_{i+1,j,k+1}}{4h_xh_z} +$

$latex + frac{beta^y_{i,j-1,k-1}-beta^y_{i,j+1,k-1}-beta^y_{i,j-1,k+1}+beta^y_{i,j+1,k+1}}{4h_yh_z} )$

$latex + frac{beta^z_{i,j,k-1}-2beta^z_{i,j,k}+beta^z_{i,j,k+1}}{h_z^2}$,

con:

$latex partial_{beta^z_{i,j,k}} F(beta^z_{i,j,k}) = -2 ( frac{1}{h_x^2} + frac{1}{h_y^2} + frac{4}{3} frac{1}{h_z^2} )$.

Tags: , , ,

Hace tiempo que no escribo nada pues estoy intentando terminar un programa que ya va tocando…

Para que no se diga, añado a continuación una imagen que he obtenido hace un momento, y que me ha hecho gracia, cuando pintaba el error entre la solución analítica de un problema de Poisson tridimensional y mi solución aproximada.

caraError

Se diría que estoy diseñando la cabeza de un nuevo personaje de animación, ¿no?

:mrgreen:

Tags: , , , , ,

Ya escribimos al respecto en este post. Aquí lo que haremos es reescribir las expresiones allí introducidas

En primer lugar, teniamos:

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

donde:

$latex S_j^* := sqrt{ frac{gamma}{f} } S = psi^6 S_j$,

$latex S_j := rho h w^2 v_j$.

En el caso de estar trabajando en cartesianas y teniendo en cuenta todo el trabajo realizado en el artículo, nos queda:

$latex partial_{xx} X^x + partial_{yy} X^x + partial_{zz} X^x = 8 pi psi^6 rho h w^2 v_x – frac{1}{3} partial_x (partial_x X^x + partial_y X^y + partial_z X^z)$,

$latex partial_{xx} X^y + partial_{yy} X^y + partial_{zz} X^y = 8 pi psi^6 rho h w^2 v_y – frac{1}{3} partial_y (partial_x X^x + partial_y X^y + partial_z X^z)$,

$latex partial_{xx} X^z + partial_{yy} X^z + partial_{zz} X^z = 8 pi psi^6 rho h w^2 v_z – frac{1}{3} partial_z (partial_x X^x + partial_y X^y + partial_z X^z)$.

A continuación, y para la siguiente ecuación, necesitamos:

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

que queda como:

$latex hat{A}^{xx} = 2 partial_x X^x – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z)$,

$latex hat{A}^{xy} = hat{A}^{yx}= partial_x X^y + partial_y X^x$,

$latex hat{A}^{xz} = hat{A}^{zx} = partial_x X^z + partial_z X^x$,

$latex hat{A}^{yy} = 2 partial_y X^y – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z)$,

$latex hat{A}^{yz} = hat{A}^{zy} = partial_y X^z + partial_z X^y$,

$latex hat{A}^{zz} = 2 partial_z X^z – frac{2}{3} (partial_x X^x + partial_y X^y + partial_z X^z)$,

por lo que:

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

donde:

$latex E^*:= sqrt{ frac{gamma}{f} } E = psi^6 E$,

$latex E:= D + tau$

es:

$latex Delta psi = -2 pi psi^{-1} (D + tau) – psi^{-7} frac{(hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2}{8}$.

La siguiente:

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

con:

$latex S^*:= sqrt{ frac{gamma}{f} } S = psi^6 S$,

$latex S:= rho h (w^2-1) + 3 p$

queda:

$latex Delta (alphapsi) = 2 pi (alphapsi)^{-1} ( D + tau + 2 rho h (w^2-1) + 6 p) + $

$latex + frac{7}{8}(alphapsi)^{-7} ((hat{A}^{xx})^2+(hat{A}^{yy})^2+(hat{A}^{zz})^2+2(hat{A}^{xy})^2+2(hat{A}^{xz})^2+2(hat{A}^{yz})^2)$

Y la última:

$latex Delta beta^i = mathcal{D}_j (2 (alphapsi)^{-6} hat{A}^{ij}) – frac{1}{3} mathcal{D}^i (mathcal{D}_j beta^j)$,

que escribimos como:

$latex Delta beta^x = partial_x (2 (alpha psi)^{-6} hat{A}^{xx}) + partial_y (2 (alpha psi)^{-6} hat{A}^{xy}) + partial_z (2 (alpha psi)^{-6} hat{A}^{xz}) – $

$latex – frac{1}{3} partial_x (partial_x beta^x + partial_y beta^y + partial_z beta^z)$

$latex Delta beta^y = partial_x (2 (alpha psi)^{-6} hat{A}^{yx}) + partial_y (2 (alpha psi)^{-6} hat{A}^{yy}) + partial_z (2 (alpha psi)^{-6} hat{A}^{yz}) – $

$latex – frac{1}{3} partial_y (partial_x beta^x + partial_y beta^y + partial_z beta^z)$

$latex Delta beta^z = partial_x (2 (alpha psi)^{-6} hat{A}^{zx}) + partial_y (2 (alpha psi)^{-6} hat{A}^{zy}) + partial_z (2 (alpha psi)^{-6} hat{A}^{zz}) – $

$latex – frac{1}{3} partial_z (partial_x beta^x + partial_y beta^y + partial_z beta^z)$

Tags: , , ,

CoCoNuT es un código que permite realizar simulaciones de colapso estelar. Reescribimos las ecuaciones CFC, que son un caso particular de la aproximación FCF haciendo que las $latex h^{ij}$ sean cero, en terminos de las variables que éste utiliza. Empezamos con una auxilar:

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

donde:

$latex S_j^* := sqrt{ frac{gamma}{f} } S = psi^6 S_j$,

$latex S_j := rho h w^2 v_j$.

La primera es:

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

donde:

$latex E^*:= sqrt{ frac{gamma}{f} } E = psi^6 E$,

$latex E:= D + tau$

La siguiente:

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

con:

$latex S^*:= sqrt{ frac{gamma}{f} } S = psi^6 S$,

$latex S:= rho h (w^2-1) + 3 p$

Y la última:

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

Además, en CFC, tenemos:

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

donde $latex L$ es el operador de Killing conforme actuando sobre la parte longitudinal $latex X^i$ sin traza y $latex A^{ij}_{TT}$ es la parte transversal sin traza de la curvatura extrínseca , y de FCF tenemos:

  • la métrica inducida en cada hipersuperficie $latex gamma_{mu nu} := g_{mu nu} + n_{mu} n_{nu}$ (o $latex boldsymbol{gamma} := boldsymbol{g} + boldsymbol{n} otimes boldsymbol{n}$ ) con $latex boldsymbol{n} = frac{dt}{|dt|}$.
  • la curvatura extrínseca $latex boldsymbol{K:=-frac{1}{2}mathcal{L}_{boldsymbol{n}} boldsymbol{gamma}}$ (o, con índices, $latex K_{mu nu} = -frac{1}{2} mathcal{L}_{boldsymbol{n}} gamma_{mu nu}$).

Tags: , , ,

Una manera sencilla de tener una ecuación de Poisson en $latex 3D$ de la que conocer su solución analítica es la siguiente. Para empezar, consideramos una función:

$latex u(x,y,z)$

a la que le aplicamos el operador $latex Delta$ y obtendremos otra función:

$latex s(x,y,z)$.

Ya tenemos $latex Delta u = s$, es decir,

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

Para las condiciones de contorno es tan sencillo como considerar el domino:

$latex [a,b] times [c,d] times [e,f]$

y ver cuanto vale $latex u$ en cada uno de los extremos, de manera que obtenemos:

$latex u(a,y,z) = g_a(y,z), u(b,y,z) = g_b(y,z)$,

$latex u(x,c,z) = g_c(x,z), u(x,d,z) = g_d(x,z)$,

$latex u(x,y,e) = g_e(x,y), u(x,y,f) = g_f(x,y)$.

Por ejemplo, si consideramos $latex u(x,y,z)=x^2+y^2+z^2$, entonces:

$latex nabla cdot nabla u = (frac{partial}{partial x},frac{partial}{partial y},frac{partial}{partial z}) cdot (u_x,u_y,u_z)$

de manera que:

$latex Delta u = frac{partial}{partial x}2x + frac{partial}{partial y}2y + frac{partial}{partial z}2z = 6$

y tenemos la ecuación de Poisson $latex Delta u = 6$. Las condiciones de contorno, en $latex Omega = [0,1]^3$ quedan:

$latex u(0,y,z) = g_{xm}(y,z) = y^2+z^2$

$latex u(1,y,z) = g_{xM}(y,z) = y^2+z^2 + 1$

$latex u(x,0,z) = g_{ym}(x,z) = x^2+z^2$

$latex u(x,1,z) = g_{yM}(x,z) = x^2+z^2 + 1$

$latex u(x,y,0) = g_{zm}(x,y) = x^2+y^2$

$latex u(x,y,1) = g_{zM}(x,y) = x^2+y^2 + 1$

Resumiendo, la solución de $latex Delta u = 6$ siendo las funciones anteriores los valores de $latex u$ en $latex partial Omega$ es:

$latex u = x^2+y^2 + z^2$.

Otro ejemplo concreto para el caso de $latex u=0$ en $latex partial Omega$ siendo

$latex Omega = { (x,y,z): 0<x<1, 0<y<1,0<z<1}$ el cubo unidad.

Tomamos $latex u(x,y,z) = (x^4-x^2)(y^4-y^2)(z^4-z^2)$. Es sencillo comprobar que $latex u=0$ en $latex partial Omega$ (p.e. $latex u(1,y,z)=(1-1)(y^4-y^2)(z^4-z^2) = 0$). ¿Cuanto vale $latex Delta u$ en este caso?

$latex s(x,y,z) = Delta [(x^4-x^2)(y^4-y^2)(z^4-z^2)] = $

$latex = nabla cdot [(4x^3-2x)(y^4-y^2)(z^4-z^2),$

$latex ,(x^4-x^2)(4y^3-2y)(z^4-z^2),$

$latex (x^4-x^2)(y^4-y^2)(4z^3-2z)] = $

$latex = 2[(6x^2-1)(y^4-y^2)(z^4-z^2) +$

$latex + (x^4-x^2)(6y^2-1)(z^4-z^2) + $

$latex + (x^4-x^2)(y^4-y^2)(6z^2-1)]$.

De manera que la solución en el cubo unidad de la ecuación de Poisson

$latex u_{xx} + u_{yy} + u_{zz} = $

$latex = 2[(6x^2-1)(y^4-y^2)(z^4-z^2) +$

$latex + (x^4-x^2)(6y^2-1)(z^4-z^2) + $

$latex + (x^4-x^2)(y^4-y^2)(6z^2-1)]$

con condiciones homogeneas de tipo Dirichlet en la frontera tiene como solución:

$latex u(x,y,z) = (x^4-x^2)(y^4-y^2)(z^4-z^2)$

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

FireStats icon Powered by FireStats