You are currently browsing the monthly archive for febrero 2013.

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

Textualmente del libro:

“As with any numerical code, debugging can be the most difficult part of creating a successful program. For multigrid, this situation is exacerbated in two ways. First, the interaction between the various multigrid components are very subtle, and it can be difficult to determine which part of a code is defective. Even more insidious is the fact that an incorrectly implemented multigrid code can perform quite well -somentimes better than other solution methods!”

Puedo confirmar, de primera mano, que ésto es así…

No tiene nada que ver con los temas habituales del blog, pero bueno, la vida es algo mas que ciencia y tecnología ¿no?. Simplemente unas frases que me gustaron de “El lado bueno de las cosas”:

“Estoy convencido de ello. Tienes que hacer todo lo que puedes y esforzarte al máximo y, si mantienes el optimismo, siempre te quedará el lado bueno de las cosas.”

“Te lo advierto. Tienes que estar atento a las señales. Cuando la vida te brinda un momento como éste, es un pecado no aprovecharlo. Es un pecado.”

Interesantes 🙂

Simplemente un reblog de nuestro crack Terry Tao sobre como capturar los conceptos esenciales de la derivación y la integración de manera algebraica para permitir su utilización sobre otros sistemas numéricos distintos de aquellos que soportan el concepto de límite, o sea, los reales y los complejos.

Posteriormente comenta como puede utilizar éstas cuando trabaja en teoría cuántica de campos para calcular integrales con variables bosónicas y fermiónicas (variables conmutativas y anticonmutativas en álgebra superconmutativa).

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)

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

En coordenadas cilíndricas podemos escribir:

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

que podemos discretizar como:

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

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

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

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

En coordenadas esféricas tenemos:

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

que podemos discretizar como:

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

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

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

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

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

En n dimensiones, el operador Laplaciano queda como:

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

en coordenadas cartesianas, y como:

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

que es equivalente a lo encontrado anteriormente.

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

\Delta u = f

como:

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

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

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

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

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

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

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

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

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

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

y al calcular el primer término:

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

que combinado con las fronteras, tenemos:

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

+ 2(f_{1,2} + f_{2,1} + f_{2,3} + f_{4,2}) + 4 f_{2,2} ] -

- \frac{1}{4}\frac{u_{0,1} + 2u_{0,2} + u_{0,3} + u_{1,0} + 2u_{2,0} + u_{3,0}}{(2h)^2}

y que es lo mismo que habíamos obtenido…

febrero 2013
L M X J V S D
« Ene   Mar »
 123
45678910
11121314151617
18192021222324
25262728