You are currently browsing the tag archive for the ‘PDE elípticas’ tag.

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

Anuncios

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

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

\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 \mathcal{D} conmuta con \Delta (métrica plana), tenemos:

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

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

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

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

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

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

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

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

donde inicialmente:

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

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

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

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

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

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

¡que vuelve a ser lineal!

Continuamos con:

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

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

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

y, finalmente:

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

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

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

donde calculamos al principio:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

y la ecuación:

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

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

\approx \frac{(\alpha\psi)_{i-1,j,k} - 2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i+1,j,k}}{h_x^2} +

+ \frac{(\alpha\psi)_{i,j-1,k}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j+1,k}}{h_y^2} +

+ \frac{(\alpha\psi)_{i,j,k-1}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j,k+1}}{h_z^2} =

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

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

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

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

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:

\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 X^i:

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

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

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

\Theta_\beta := \mathcal{D}_i \beta^i,

que discretizada queda:

\frac{(\Theta_\beta)_{i-1,j,k}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i+1,j,k}}{h_x^2} +

\frac{(\Theta_\beta)_{i,j-1,k}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i,j+1,k}}{h_y^2} +

\frac{(\Theta_\beta)_{i,j,k-1}-2(\Theta_\beta)_{i,j,k}+(\Theta_\beta)_{i,j,k+1}}{h_z^2} =

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

\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}{3} \partial_x \Theta_\beta \approx

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

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

- \frac{1}{3} (\partial_x \Theta_\beta)_{i,j,k}.

De la misma manera:

\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

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

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

- \frac{1}{3} (\partial_y \Theta_\beta)_{i,j,k}.

Y, por último:

\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

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

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

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

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

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

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

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

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

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

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

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

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

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

\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

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

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

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

+ \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 F(X^x_{i,j,k})=0 y entonces tenemos:

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

\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

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

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

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

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

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

\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

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

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

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

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

con:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\Delta (\alpha\psi) = 2 \pi (\alpha\psi)^{-1} ( D + \tau + 2 \rho h (w^2-1) + 6 p) +

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

\approx \frac{(\alpha\psi)_{i-1,j,k} - 2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i+1,j,k}}{h_x^2} + \frac{(\alpha\psi)_{i,j-1,k}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j+1,k}}{h_y^2} + \frac{(\alpha\psi)_{i,j,k-1}-2(\alpha\psi)_{i,j,k}+(\alpha\psi)_{i,j,k+1}}{h_z^2} =

=2 \pi (\alpha\psi)_{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}) +

+ \frac{7}{8}(\alpha\psi)_{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:

\partial_{\psi\alpha_{i,j,k}} F(\psi\alpha_{i,j,k}) = -2 ( \frac{1}{h_x^2} + \frac{1}{h_y^2} + \frac{1}{h_z^2} ) +

+ 2 \pi (\psi\alpha)_{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}) -

+ \frac{49}{8} (\psi\alpha)_{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:

\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}{3} \partial_x (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

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

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

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

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

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

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

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

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

\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 (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

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

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

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

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

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

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

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

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

\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 (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z) \approx

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

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

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

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

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

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

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

con:

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

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:

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

En primer lugar, teniamos:

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

donde:

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

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:

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

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

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

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

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

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

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

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

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

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

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

donde:

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

E:= D + \tau

es:

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

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

con:

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

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

queda:

\Delta (\alpha\psi) = 2 \pi (\alpha\psi)^{-1} ( D + \tau + 2 \rho h (w^2-1) + 6 p) +

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

Y la última:

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

que escribimos como:

\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}{3} \partial_x (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

\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 (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

\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 (\partial_x \beta^x + \partial_y \beta^y + \partial_z \beta^z)

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 h^{ij} sean cero, en terminos de las variables que éste utiliza. Empezamos con una auxilar:

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

donde:

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

S_j := \rho h w^2 v_j.

La primera es:

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

donde:

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

E:= D + \tau

La siguiente:

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

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

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

Y la última:

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

\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 L es el operador de Killing conforme actuando sobre la parte longitudinal X^i sin traza y 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 \gamma_{\mu \nu} := g_{\mu \nu} + n_{\mu} n_{\nu} (o \boldsymbol{\gamma} := \boldsymbol{g} + \boldsymbol{n} \otimes \boldsymbol{n} ) con \boldsymbol{n} = \frac{dt}{|dt|}.
  • la curvatura extrínseca \boldsymbol{K:=-\frac{1}{2}\mathcal{L}_{\boldsymbol{n}} \boldsymbol{\gamma}} (o, con índices, K_{\mu \nu} = -\frac{1}{2} \mathcal{L}_{\boldsymbol{n}} \gamma_{\mu \nu}).

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

u(x,y,z)

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

s(x,y,z).

Ya tenemos \Delta u = s, es decir,

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

[a,b] \times [c,d] \times [e,f]

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

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

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

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

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

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

\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 \Delta u = 6. Las condiciones de contorno, en \Omega = [0,1]^3 quedan:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

+ (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

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

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

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

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

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

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

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

Empezamos suponiendo que conocemos:

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

\frac{\partial}{\partial y}|_{0,0}u, \frac{\partial}{\partial y}|_{1,0}u

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

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

Discretizamos:

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

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

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

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

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

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

\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

\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

\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

\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

\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

\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

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

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

\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 6 ecuaciones con 6 incognitas y la matriz tiene rango 6, por lo que la solución es única.

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

\frac{\partial}{\partial x}|_{0,0}u, \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}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}|_{0,2}u, \frac{\partial}{\partial y}|_{1,2}u, \frac{\partial}{\partial y}|_{2,2}u

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

Si discretizamos:

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

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

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

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

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

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

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

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

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

\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

\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

\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

\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

\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

\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

\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

\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

\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

\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

\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

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

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

\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 9 ecuaciones con 9 incognitas pero la matriz tiene rango 8, por lo que tenemos infinitas soluciones. Hay que conservar.

Suponemos \Delta u = f en 2D, es decir,

\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 x en una frontera, condición Neumann respecto y en una frontera y condición Neumann respecto x e y en dos fronteras.

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

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, 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 9 por lo que existe solución y es única. Desde el punto de vista algebraico, es el punto (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 9 hiperplanos

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

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

\ldots

en el espacio \mathbb{R}^9.

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

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

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 12 incognitas (u_{i,j} con i=0..3 y j=1..3) y el rango de A_2 es 12, por lo que la solución, nuevamente, es única.

Para el caso en el que conocemos \frac{\partial}{\partial y}|_{1,0}u, \frac{\partial}{\partial y}|_{2,0}u, \frac{\partial}{\partial y}|_{3,0}u en lugar de 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:

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:

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 12 ecuaciones con 12 incognitas (u_{i,j} con i=1..3 y j=0..3) y el rango de A_3 es 12, por lo que la solución es única.

Finalmente, suponemos conocidos \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 7 ecuaciones mas a las 9 que ya teniamos por lo que nos queda una matrix A_4 \in \mathcal{M}(16 \times 16):

\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 u_{0,0} por 4, y las correspondientes a u_{0,1}, u_{0,2}, u_{0,3}, u_{1,0},u_{2,0}, u_{3,0}  por 2, quedando:

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

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

\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 i,j=1,1,

\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 i,j=1,2,

\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 i,j=1,3,

\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 i,j=2,1,

\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 i,j=2,2,

\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 i,j=2,3,

\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 i,j=3,1,

\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 i,j=3,2,

\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 i,j=3,3,

de donde:

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

\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 u_{0,1}, u_{0,2}, u_{0,3} conocemos \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:

\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 i,j=0,1

\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 i,j=0,2

\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 i,j=0,3

y

\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

\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

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

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

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

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

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

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 1D supongamos que ahora tenemos \frac{\partial^2}{\partial x^2}u = f en [a,b] con u(a)=u_a pero \frac{\partial}{\partial x} = du_b. Suponiendo de nuevo n=8, las ecuaciones nos quedan:

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

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

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

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

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

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

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

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

\frac{u_1 - u_{-1}}{2h} = \frac{d}{dx}u_{0} = du_0,

que, despejando, nos da:

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:

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

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

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

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

\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 -2 y -1 en la diagonal en función de si es Dirichlet o Neumann.

diciembre 2017
L M X J V S D
« Ago    
 123
45678910
11121314151617
18192021222324
25262728293031