You are currently browsing the category archive for the ‘Métodos Numéricos’ category.

Después de muchísimo tiempo sin escribir, vuelvo con este post que resume todo el trabajo que hemos realizado durante este último año.

En julio del año pasado nos enteramos de un artículo publicado por X. Yang and R. Mittal, de la Johns Hopkins University, en el que aceleraban de manera espectacular el algoritmo de Jacobi y lo utilizaban para resolver ecuaciones en derivadas parciales de tipo elíptico.

A pesar del relativo eco mediático que tuvo, y aunque aceleraba muchísimo Jacobi, seguía sin ser competitivo con los métodos utilizados actualmente para solucionar este tipo de ecuaciones. Sin embargo, como prometía en cuanto a su sencillez desde el punto de vista tanto de implementación como de paralelización, decidimos trabajar un tiempo sobre el mismo.

Finalmente, en junio presentamos unos proceedings en el CEDYA 2015 (pag. 733) y hace tres semanas enviamos un paper a JCP, donde presentamos una serie de mejoras realizadas que permiten, tal y como allí comentamos, que el SRJ sea prácticamente competitivos en ejecución secuencial con los algoritmos utilizados actualmente.

Por un lado, su inmediata implementación (todos tenemos un Jacobi implementado de nuestro curso de métodos numéricos, y es trivial convertirlo en un SRJ 🙂 ) hace pensar que mucha gente que no tenga y necesite un resolvedor elíptico eficiente sin invertir mucho esfuerzo en su implementación quiera utilizarlo. Por otro, su trivial paralelización, por ejemplo en entornos GPU, nos hace pensar en su extraordinario potencial en el ámbito de la supercomputación.

En la página de nuestro grupo de investigación tenemos disponibles todos los esquemas SRJ presentados en el paper.

¡Disfrutenlo!

Anuncios

Tenemos:

  1. \bar{r} := \frac{r}{r+a}
  2. \Delta := \frac{(1-\bar{r})^4}{a^2} \partial_{\bar{r}\bar{r}} + \frac{(1-\bar{r})^4}{a^2}\frac{2}{\bar{r}} \partial_{\bar{r}}
  3. \Delta \Theta_X = 6 \pi [ \frac{(1-\bar{r})^2}{a} \partial_{\bar{r}} S^*_{\bar{r}} + \frac{1 - \bar{r}}{a} \frac{2}{\bar{r}} S^*_{\bar{r}} + \frac{1-\bar{r}}{a} \frac{\cot \theta}{\bar{r}} S^*_{\theta}]
  4. \Delta X^{\bar{r}} = 8 \pi S^*_{\bar{r}} - \frac{1}{3} \frac{(1-\bar{r})^2}{a} \partial_{\bar{r}} \Theta_X
  5. \hat{A}^{\bar{r}\bar{r}} = \frac{4}{3}\frac{(1-\bar{r})^2}{a} \partial_{\bar{r}} X^{\bar{r}} - \frac{2}{3}2\frac{1-\bar{r}}{a}\frac{1}{\bar{r}} X^{\bar{r}}
  6. \Delta \Theta_\beta = \frac{3}{2}[\frac{(1-\bar{r})^4}{a^2} \partial_{\bar{r}\bar{r}}u + \frac{(1-\bar{r})^3(2-\bar{r})}{2a^2} \frac{4}{\bar{r}}u + \frac{(1-\bar{r})^2}{a^2} \frac{2}{\bar{r}^2}u]

con:

u:=\alpha \psi^{-6} \hat{A}^{\bar{r}\bar{r}}

y:

\{\frac{2 (\bar{r}_i - 1)^4 (\bar{r_i} - (\bar{r}_{i+1} - \bar{r}_i))}{a^2 (\bar{r}_i - \bar{r}_{i-1}) \bar{r_i} (\bar{r}_{i+1} - \bar{r}_{i-1} )},

\frac{(\bar{r}_i - 1)^2 (\frac{-2}{h_\theta^2 \bar{r}_i^2} + \frac{(\bar{r}_i - 1)^2 ((r_{i+1}-r_i)-(r_i-r_{i-1})-2)}{(r_{i+1}-r_i)(r_i - r_{i-1})} + \frac{-2 \csc^2 \theta_i}{h_\varphi^2 \bar{r}_i^2})}{a^2},

\frac{2 (\bar{r}_i - 1)^4 ((\bar{r}_{i} - \bar{r}_{i-1}) + \bar{r_i})}{a^2 (\bar{r}_{i+1} - \bar{r}_i) \bar{r}_i (\bar{r}_{i+1} - \bar{r}_{i-1} )} \}

Para aproximar la primera y segunda derivada de una función f(x) mediante tres puntos estamos habituados a las fórmulas:

f'(x) \approx \frac{f_{i+1}-f_{i-1}}{2h} = \frac{-1}{2h} f_{i-1} + \frac{1}{2h} f_{i+1},

f''(x) \approx \frac{f_{i-1} - 2f_i + f_{i+1}}{h^2} = \frac{1}{h^2} f_{i-1} + \frac{-2}{h^2} f_i + \frac{1}{h^2} f_{i+1}.

En estas expresiones estamos asumiendo que los puntos están equiespaciados una distancia h. ¿Cómo quedan las formulas en el caso de que la distancia entre los dos primeros puntos lx sea diferente a la distancia entre los dos últimos rx? Existen varias maneras de calcularlo, por ejemplo mediante interpolación de Lagrange como ya hicimos en este post, y quedan:

f'(x) \approx \frac{-rx}{lx(lx+rx)} f_{i-1} + \frac{rx - lx}{lx rx} f_i + \frac{lx}{(lx+rx)rx} f_{i+1},

f''(x) \approx \frac{2}{lx(lx+rx)} f_{i-1} + \frac{-2}{lx rx} f_i + \frac{2}{(lx+rx)rx} f_{i+1}.

Supongamos que tenemos el Laplaciano expresado en un sistema de coordenadas curvilineas cualesquiera:

\bigg [ c_1 \frac{\partial^2}{\partial x_1^2} + c_2 \frac{\partial}{\partial x_1} +c_3 \frac{\partial^2}{\partial x_2^2} + c_4 \frac{\partial}{\partial x_2} + c_5 \frac{\partial^2}{\partial x_3^2} + c_6 \frac{\partial}{\partial x_3} \bigg ] u(x_1,x_2,x_3).

La expresión correspondiente en diferencias finitas, utilizando las aproximaciones para las primeras y segundas derivadas de este post, nos queda:

c_{1_{i,j,k}}\frac{1}{h_{x_1}^2} (u_{i-1,j,k} - 2 u_{i,j,k} + u_{i+1,j,k}) + c_{2_{i,j,k}} \frac{1}{2h_{x_1}}(u_{i+1,j,k}-u_{i-1,j,k}) +

+ c_{3_{i,j,k}}\frac{1}{h_{x_2}^2} (u_{i,j-1,k} - 2 u_{i,j,k} + u_{i,j+1,k}) + c_{4_{i,j,k}} \frac{1}{2h_{x_2}}(u_{i,j+1,k}-u_{i,j-1,k}) +

+ c_{5_{i,j,k}}\frac{1}{h_{x_3}^2} (u_{i,j,k-1} - 2 u_{i,j,k} + u_{i,j,k+1}) + c_{6_{i,j,k}} \frac{1}{2h_{x_3}}(u_{i,j,k+1}-u_{i,j,k-1}),

que, agrupando por términos, podemos reescribir como:

(c_{1_{i,j,k}} \frac{1}{h_{x_{1}}^2} - c_{2_{i,j,k}} \frac{1}{2h_{x_{1}}}) u_{i-1,j,k} + (c_{1_{i,j,k}} \frac{1}{h_{x_{1}}^2} + c_{2_{i,j,k}} \frac{1}{2h_{x_{1}}}) u_{i+1,j,k} +

+ (c_{3_{i,j,k}} \frac{1}{h_{x_{2}}^2} - c_{4_{i,j,k}} \frac{1}{2h_{x_{2}}}) u_{i,j-1,k} + (c_{3_{i,j,k}} \frac{1}{h_{x_{2}}^2} + c_{4_{i,j,k}} \frac{1}{2h_{x_{2}}}) u_{i,j+1,k} +

+ (c_{5_{i,j,k}} \frac{1}{h_{x_{3}}^2} - c_{6_{i,j,k}} \frac{1}{2h_{x_{3}}}) u_{i,j,k-1} + (c_{5_{i,j,k}} \frac{1}{h_{x_{3}}^2} + c_{6_{i,j,k}} \frac{1}{2h_{x_{3}}}) u_{i,j,k+1} -

-2 (c_{1_{i,j,k}} \frac{1}{h_{x_{1}}^2} + c_{3_{i,j,k}} \frac{1}{h_{x_{2}}^2} + c_{5_{i,j,k}} \frac{1}{h_{x_{3}}^2}) u_{i,j,k}.

He escrito un modulo en Mathematica que, a partir de los coeficientes del Laplaciano y de las distancias para la discretización nos generan automáticamente estos coeficientes:

finDif

Los anteriores ejemplos corresponden al caso de cilíndricas, esféricas y esféricas compactificadas.

psiAlpha

phiAlpha1

psiAlpha2

psiAlpha3

El caso mas sencillo es cuando tenemos cinco puntos:

(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3,y_3), (x_4,y_4),

de manera que:

x_1-x_0 = 2(x_2-x_1) = 2(x_3-x_2) = x_4-x_3.

Tal y como escribimos en el anterior post, el polinomio general para cinco puntos quedaria:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x) + y_3 l_3(x) + y_4 l_4(x),

donde:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4},

l_1(x) = \frac{x-x_0}{x_1-x_0} \frac{x-x_2}{x_1-x_2} \frac{x-x_3}{x_1-x_3} \frac{x-x_4}{x_1-x_4},

l_2(x) = \frac{x-x_0}{x_2-x_0} \frac{x-x_1}{x_2-x_1} \frac{x-x_3}{x_2-x_3} \frac{x-x_4}{x_2-x_4},

l_3(x) = \frac{x-x_0}{x_3-x_0} \frac{x-x_1}{x_3-x_1} \frac{x-x_2}{x_3-x_2} \frac{x-x_4}{x_3-x_4},

l_4(x) = \frac{x-x_0}{x_4-x_0} \frac{x-x_1}{x_4-x_1} \frac{x-x_2}{x_4-x_2} \frac{x-x_3}{x_4-x_3}.

Dado que tratamos de aproximar la segunda derivada, tenemos, en el caso de equidistancia quedaría:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12 h^2}.

Vamos a ver que queda ahora en nuestro caso. Reescribimos los l_i(x) en función de x_0 de manera que, por ejemplo, tenemos:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4} =

=\frac{x-x_0-2h}{x_0-x_0-2h} \frac{x-x_0-3h}{x_0-x_0-3h} \frac{x-x_0-4h}{x_0-x_0-4h} \frac{x-x_0-6h}{x_0-x_0-6h} =

=\frac{(x-x_0-2h)(x-x_0-3h)(x-x_0-4h)(x-x_0-6h)}{144h^4},

que para x=x_2=x_0+3h (2h+h), la segunda derivada queda:

\frac{d^2}{dx^2} l_0(x) |_{x=x_2} = -\frac{1}{72h^2}.

Haciendo lo mismo para todas las derivadas segundas de todos los l_i(x), obtenemos la siguiente formula en diferencias finitas:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-3} + 81 u_{u-1} - 160 u_i + 81 u_{i+1} - u_{i+3}}{72h^2},

donde los índices hacen referencia a una malla inicial equiespaciada.

Para el mismo denominador, antes teniamos:

\frac{d^2}{dx^2}u_i = \frac{-6u_{i-2}+96u_{i-1}-180u_i+96u_{i+1}-6u_{i+2}}{72 h^2}.

Podemos hacer lo mismo para cada uno de los puntos x=x_0, x=x_1=x_0+2h, x=x_3=x_0+4h y x=x_4=x_0+6h:

\frac{d^2}{dx^2}u_{i-3} = \frac{40 u_{i-3} -243 u_{i-1} + 352 u_i -162 u_{i+1} + 13 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i-2} = \frac{7 u_{i-3} - 32 u_i + 27 u_{i+1} - 2 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+1} =\frac{-2 u_{i-3} + 27 u_{i-1} - 32 u_i -+7 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+3} =\frac{13 u_{i-3} -162 u_{i-1} + 352 u_i -243 u_{i+1} + 40 u_{i+3}}{36 h^2}

 En el caso de x=x_0+4h (3h+h), tenemos:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-4} + 256_{u-1} - 510 u_i + 256 u_{i+1} - u_{i+4}}{140h^2}.

Dados n+1 puntos:

(x_0,y_0), \ldots, (x_n,y_n),

definimos el polinomio interpolador de Lagrange:

L(x) := \Sigma_{i=0}^n y_i l_i(x)

donde:

l_i(x) := \Pi_{0 \leq j \leq n,j \neq i} \frac{x-x_j}{x_i-x_j}.

Con esto, tenemos:

\frac{d^n}{dx^n}L(x) = \Sigma_{i=0}^n y_i \frac{d^n}{dx^n}l_i(x)

Supongamos que tenemos n=2:

(x_0,y_0), (x_1,y_1).

En este caso, tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x)

con:

l_0(x) = \frac{x-x_1}{x_0-x_1}

l_1(x) = \frac{x-x_0}{x_1-x_0}

Como tenemos dos puntos, únicamente podemos aproximar la primera derivada:

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x)

con:

\frac{d}{dx} l_0(x) = \frac{1}{x_0-x_1}

\frac{d}{dx} l_1(x) = \frac{1}{x_1-x_0},

de manera que:

\frac{d}{dx}L(x) = \frac{y_1 - y_0}{x_1 - x_0}.

Con esto, tenemos las aproximaciones de primer orden:

u_x(x_0) = u_x(x_1) = \frac{y_1 - y_0}{x_1-x_0},

que, en índices, queda:

\frac{d}{dx}u_i \approx \frac{u_{i+1}-u_i}{h_x},

\frac{d}{dx}u_{i+1} \approx \frac{u_{i+1}-u_i}{h_x},

donde h_x = x_1 - x_0 (y que es lógico, ya que la derivada será una constante).

Supongamos ahora que tenemos tres puntos:

(x_0,y_0),(x_1,y_1),(x_2,y_2).

En este caso tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x),

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x) + y_2 \frac{d}{dx} l_2(x),

\frac{d^2}{dx^2} L(x) = y_0 \frac{d^2}{dx^2} l_0(x) + y_1 \frac{d^2}{dx^2} l_1(x) + y_2 \frac{d^2}{dx^2} l_2(x).

Ahora tenemos:

l_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{x^2 + (-x_1-x_2)x + x_1x_2}{(x_0-x_1)(x_0-x_2)},

l_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{x^2 + (-x_0-x_2)x + x_0x_2}{(x_1-x_0)(x_1-x_2)},

l_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{x^2 + (-x_0-x_1)x + x_0x_1}{(x_2-x_0)(x_2-x_1)}.

\frac{d}{dx} l_0(x) = \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)},

\frac{d}{dx} l_1(x) = \frac{2x+(-x_0-x_2)}{(x_1-x_0)(x_1-x_2)},

\frac{d}{dx} l_2(x) = \frac{2x+(-x_0-x_1)}{(x_2-x_0)(x_2-x_1)}.

\frac{d^2}{dx^2} l_0(x) = \frac{2}{(x_0-x_1)(x_0-x_2)},

\frac{d^2}{dx^2} l_1(x) = \frac{2}{(x_1-x_0)(x_1-x_2)},

\frac{d^2}{dx^2} l_2(x) = \frac{2}{(x_2-x_0)(x_2-x_1)}.

De esta manera, por ejemplo, tenemos que:

\frac{d^2}{dx^2}L(x) = y_0 \frac{2}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2}{(x_2-x_0)(x_2-x_1)},

que, en el caso de tener los puntos equiespaciados (h_x:=x_{i+1}-x_{i} con 0 \leq i \leq 1), queda la aproximación de segundo orden:

\frac{d^2}{dx^2}u(x) = \frac{y_0 - 2 y_1 + y_2}{h_x^2},

que además, como era de esperar, es independiente de x:

\frac{d^2}{dx^2}u_{i-1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_{i+1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}.

¿Qué pasa con la primera derivada? En este caso, el resultado si que depende de x (por lo que tendremos un resultado diferente en función de si la x vale x_0, x_1 o x_2) y lo que obtenemos es:

\frac{d}{dx} L(x) = y_0 \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

por lo que:

\frac{d}{dx} L(x_0) = y_0 \frac{2x_0 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_0 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_0 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_1) = y_0 \frac{2x_1 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_1 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_1 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_2) = y_0 \frac{2x_2 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_2 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_2 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

que, con equiespaciado, quedan:

L'(x_0) = \frac{3 y_0 -4 y_1 + y_2}{2 h_x},

L'(x_1) = \frac{y_0 -2 y_1 + y_2}{2 h_x} y

L'(x_2) = \frac{y_0 -4 y_1 + 3 y_2}{2 h_x}.

Laplaciano en cartesianas:

\Delta u = \Sigma_i \frac{\partial^2}{\partial x_i^2}u

1d

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

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

2d

\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2} + \frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2} = f_{i,j}

i fijo:

\frac{1}{h_y^2}u_{i,j-1} + \frac{1}{h_y^2}u_{i,j+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_x^2}u_{i-1,j} + \frac{-1}{h_x^2}u_{i+1,j})

j fijo:

\frac{1}{h_x^2}u_{i-1,j} + \frac{1}{h_x^2}u_{i+1,j} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_y^2}u_{i,j-1} + \frac{-1}{h_y^2}u_{i,j+1})

3d

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

i,j fijos:

\frac{1}{h_z^2}u_{i,j,k-1} + \frac{1}{h_z^2}u_{i,j,k+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k})

i,k fijos:

\frac{1}{h_y^2}u_{i,j-1,k} + \frac{1}{h_y^2}u_{i,j+1,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

j,k fijos:

\frac{1}{h_x^2}u_{i-1,j,k} + \frac{1}{h_x^2}u_{i+1,j,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

Ya hace casi dos meses que no escribo nada en el blog :-(. Es hora de añadir algo.

El Block Cyclic Reduction es un algoritmo rápido para resolver sistemas de ecuaciones de manera directa cuya matriz tiene estructura tridiagonal por bloques.

Supongamos que tenemos una discretización de un volumen mediante cinco puntos para las x, tres puntos para las y y siete puntos para las z (necesitamos que sean números impares).

Para la dimensión z tenemos:

A_3 = \left(  \begin{array}{ccccccc}  6 & -1 & 0 & 0 & 0 & 0 & 0 \\  -1 & 6 & -1 & 0 & 0 & 0 & 0 \\  0 & -1 & 6 & -1 & 0 & 0 & 0 \\  0 & 0 & -1 & 6 & -1 & 0 & 0\\  0 & 0 & 0 & -1 & 6 & -1 & 0 \\  0 & 0 & 0 & 0 & -1 & 6 & -1\\  0 & 0 & 0 & 0 & 0 & -1 & 6  \end{array}  \right)

-I_3 =\left(  \begin{array}{ccccccc}  -1 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & -1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & -1 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & -1 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & -1 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & -1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & -1  \end{array}  \right)

Para la dimensión y no queda:

A_2 = \left(  \begin{array}{ccc}  A_3 & -I_3 & 0 \\  -I_3 & A_3 & -I_3 \\  0 & -I_3 & A_3  \end{array}  \right); \,\,\, -I_2 = \left(  \begin{array}{ccc}  -I_3 & 0 & 0 \\  0 & -I_3 & 0 \\  0 & 0 & -I_3  \end{array}  \right)

Y, finalmente, para la dimensión x obtenemos:

A_1 = \left(  \begin{array}{ccccc}  A_2 & -I_2 & 0 & 0 & 0 \\  -I_2 & A_2 & -I_2 & 0 & 0 \\  0 & -I_2 & A_2 & -I_2 & 0 \\  0 & 0 & -I_2 & A_2 & -I_2 \\  0 & 0 & 0 & -I_2 & A_2  \end{array}  \right).

Ahora todo es fácilmente generalizable a n dimensiones x_1, x_2, \ldots, x_n con t_1, t_2, \ldots, t_n puntos en cada dimensión, ya que nos queda A_i una matriz tridiagonal por bloques con t_i \times t_i bloques donde cada uno está formado por matrices A_{i-1} en la diagonal, -I_{i-1} arriba y debajo de la diagonal y 0 en el resto, todas de dimensión t_{i-1} \times t_{i-1} con i=1,\ldots, n. En A_n tenemos el stencil (-1,2n,-1).

En el caso que estamos tratando, para cada fila par de A_1 podemos multiplicarla por A_2 y sumarla a las filas impares anterior y posterior, de manera que eliminamos las variables pares:

-v^1_{j-2} + ((A_2)^2 - 2I_2)v^1_j - v^1_{j+2} = b^1_{j-1} + A_2 b^1_j + b^1_{j+1},

(el ^1 hace referencia a las x y no corresponde a potencias… Para resolver las k impares simplemente:

A_2 v^1_k = b^1_{k} + v^1_{k-1} + v^1_{k+1}. )

Y resulta que definiendo B_2:=(A_2)^2 - 2 I_2, nos queda una matriz B_1 de la misma estructura que A_1, por lo que podemos aplicar este metodo de manera recursiva hasta llegar a tener una sola ecuación (para esto, no basta con que tengamos un número impar de nodos, necesitamos algo como N=2^{k+1}-1).

Llegados a este punto, si estamos con v^n, resolvemos como queramos, pero si estamos en otra variable v^i con i\neq n, la matriz vuelve a tener la misma estructura, pero con una dimensión menos, y volvemos a aplicar el mismo algoritmo.

Tenemos dos lugares distintos donde aplicar recursividad: el primero por ir agrupando de tres en tres las ecuaciones de una dimensión de trabajo (A_i \rightarrow B_i \rightarrow \ldots \rightarrow F_i, F por poner algo :-)); el otro cuando, llegados a una sola ecuación aplicando la recursividad anterior, nos quedan mas dimensiones por resolver (B_{j+1}, B_{j+2}, \ldots ).

Para empezar, vamos a calcular numéricamente, y en mecánica clásica,  el campo creado por una partícula con masa en el espacio.

Si colocamos las fronteras de nuestro dominio \Omega \in \mathbb{R}^3 lo suficientemente lejos, el cálculo del potencial gravitatorio se reduce a resolver la ecuación de Poisson

\Delta u(x,y,z) = 4 \pi G \rho(x,y,z) si (x,y,z) \in \Omega

con

u(\bar{x},\bar{y},\bar{z}) = 0 si (\bar{x},\bar{y},\bar{z}) \in \partial \Omega,

y donde \rho(x,y,z) es la densidad de masa.

Aunque en estos casos tenemos solución analítica utilizando la función de Green, vamos a calcularla numéricamente utilizando Chombo, y mediante superposición, extenderlo a un sistema de masas puntuales.

Para empezar con la aproximación, utilizaremos gaussianas para simular las funciones \delta correspondientes a la distribución de masa puntual (mas similares a una \delta cuanto mas estrechas sean).

A continuación, colocaremos una partícula en el centro de un cubo y utilizaremos un malla adaptativa, definida manualmente (aunque Chombo también capaz de hacerlo automáticamente), que se hará mas fina a medida que se acerque a la misma. Esto permitirá tener resolución y tiempo de cálculo donde realmente se necesita:

masPun1masPun2

Colocamos fronteras Dirichlet homogeneas, y ejecutamos. El resultado es:

masPunSol1

que, al ser tridimensional, cuesta un poco de ver. Básicamente, es un campo con simétria esférica. A continuación cortamos con un plano por el centro del campo

masPunSol2masPunSol2b masPunSol2c

y lo elevamos, pues es campo gravitatorio se puede pensar como una curvatura 🙂

masPunSol3b masPunSol3

Suponemos ahora dos masas puntuales y procedemos de la misma manera. El tamaño de la partícula representa su masa.

2masPunSol1

Como se puede apreciar, hemos añadido un poco mas de resolución en la parte de la partícula mas masiva. Lo que obtenemos es (en 3D y en corte):

Finalmente, si elevamos los cortes, obtenemos:

2masPunSol3 2masPunSol3b

Finalmente, algunas gráficas correspondientes a dos nuevos casos en los que separo cada vez mas las partículas:

En el libro de Análisis numérico de Burden y Faires aparecen una serie de ejemplos y ejercicios resueltos de PDEs elípticas en 2D. Vamos a intentar resolverlos utilizando Chombo…

La primera ecuación corresponde con un ejemplo y es

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

en

\Omega = \{ (x,y): 0<x<0.5, 0<y<0.5\}

y cuyas condiciones en la frontera \partial \Omega son

u(x,0)=0, u(0,y)=0, u(x,0.5)=200x, u(0.5,y)=200y.

El resultado es:

ejeM1

Para el ejercicio 12.1.1

u_{xx} + u_{yy} = 4

en

\Omega = \{ (x,y): 0<x<1, 0<y<2\}

y cuyas condiciones en la frontera \partial \Omega son

u(x,0)=x^2, u(0,y)=y^2, u(x,2)=(x-2)^2, u(1,y)=(y-1)^2.

El resultado es:

ejeR1

A continuación, para el 12.1.3.a

u_{xx} + u_{yy} = 0

en

\Omega = \{ (x,y): 0<x<1, 0<y<1\}

y cuyas condiciones en la frontera \partial \Omega son

u(x,0)=0, u(0,y)=0, u(x,1)=x, u(1,y)=y.

El resultado es:

ejeR3a

En el 12.1.3.b encontramos

\Delta u = -(\cos(x+y)+\cos(x-y))

en

\Omega = \{ (x,y): 0<x<\pi, 0<y<\frac{\pi}{2}\}

y

u(x,0)=\cos x , u(0,y)=\cos y, u(x,\frac{\pi}{2})=0, u(\pi,y)=-\cos y,

obteniendo:

 ejeR3b ejeR3b2

Definimos a los kernels como funciones del tipo:

W_{ab}=W(\boldsymbol{r}_a - \boldsymbol{r}_b,h),

donde a es la partícula en la que está centrada la función y b es una partícula dentro del soporte compacto de la función kernel, controlado éste último por h, la smoothing length (longitud de suavizado).

En este post básicamente pretendo aclarar lo que significa \nabla_a W_{ab} cuando, por ejemplo, tenemos definido el kernel como:

W(q) = \alpha_D \exp (-q^2) con 0 \leq q \leq 2.

Para empezar, \alpha_D es una constante de dimensionalidad, por lo que la fórmula está escrita de manera compacta y sirve para cualquier dimensión. Además, tenemos que q = \frac{r}{h}, siendo r la distancia ente las partículas, por lo que:

r = |\boldsymbol{r}_a - \boldsymbol{r}_b| =^{(3D)} \sqrt{(x_a-x_b)^2 + (y_a-y_b)^2 + (z_a-z_b)^2}.

Si fijamos la posición de la partícula a, la función que nos da la distancia de esta a cualquier punto dentro del soporte compacto es:

r_a (\boldsymbol{r}) = |\boldsymbol{r}_a-\boldsymbol{r}| =^{(3D)} \sqrt{(x_a-x)^2 + (y_a-y)^2 + (z_a-z)^2},

siendo q_a lo mismo añadiendo el factor h.

Por lo tanto, en este caso tenemos, en tres dimensiones y donde b es una partícula en una posición arbitraria (x,y,z):

\nabla_a W_{ab}(q) =^{(3D)} (\partial_x W_{ab}(q_a), \partial_y W_{ab}(q_a), \partial_z W_{ab}(q_a)) =

= \alpha_D \exp(-q^2) (-2q) (\partial_x (q_a), \partial_y (q_a), \partial_z (q_a)) = \alpha_D \exp(-q^2) (-2q) \nabla_a q_a

donde:

\nabla_a q_a = \frac{-1}{h r_a} (x_a-x,y_a-y,z_a-z).

De la misma manera, si tenemos:

W(q) = \alpha_D \begin{cases} 1-\frac{3}{2} q^2 + \frac{3}{4} q^3, 0 \leq q < 1\\ \frac{1}{4} (2-q)^3, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

entonces:

\nabla_a W_{ab}(q) = \alpha_D \begin{cases} (-3q + \frac{9}{4}q^2) \nabla_a q_a, 0 \leq q < 1 \\ -\frac{3}{4} (2-q)^2 \nabla_a q_a, 1 \leq q < 2 \\ 0, q \geq 2 \end{cases}

Así pues:

\nabla W(q) = \frac{d}{dq} W(q) \nabla q.

 

A la hora de resolver las diferentes ecuaciones elípticas CFC tenemos dos posibilidades para fijar las condiciones en la frontera, cada una con sus mas y sus menos.

La primera consiste en hacer un desarrollo multipolar de los terminos fuente en armónicos esféricos, de manera que cuantos mas términos consideremos mas cerca podremos colocar la frontera.

La segunda consiste en compactificar el dominio, lo que nos permite reducir todo el universo a un cubo unidad y considerar Minkowski en su frontera, puesto que ésta corresponde a infinito.

Hace un tiempo estuve trabajando con autómatas celulares y tengo ganas de compartir algunas cosas interesantes de estos, aunque ahora ando un poco escaso de tiempo…

Para empezar a abrir boca, un enlace al artículo original de John von Neumann “Theory of Self-Reproducing Automata” y otros dos a la wikipedia: la definición de lo que son y el ya clásico juego de la vida.

Al resolver mediante diferencias finitas PDEs lo que estamos haciendo, de alguna manera, es poner en marcha un pseudoautómata celular que lo hace. Me interesa explorar sus generalizaciones: del estado, de la función de transición (como podemos encontrarlas, por ejemplo a partir de algoritmos genéticos), de la vecindad (trabajar con vecinos no contíguos),…

Por cierto, el record del Bounded gaps between primes podría estar ya en 6966, por lo que hemos rebajado el problema en 4 ordenes de magnitud (Iniciamos el viaje en 70,000,000 y habría que llegar a 2 para demostrar la conjetura de los números primos gemelos). Ya queda menos :-).

Algunas imagenes VisIt resultantes de interpolaciones entre mallas esféricas y mallas cartesianas de componentes de la métrica en el formalismo 3+1:

  • \beta^z (pseudocolor):
  • \alpha (contour):

alphaSphCar

Una nueva imagen curiosa del error en la resolución numérica de una Poisson 2d utilizando multigrid…

caraError2d

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:

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…

Llego ya un tiempo implementando los diferentes esquemas multigrid en C++ y, tras un periodo de larga oscuridad, por fin aparece un poco de luz: ayer , por primera vez, pareció funcionar todo (en una dimensión y con fronteras Dirichlet, claro… Aunque todo lo demás ya está también implementado).

En el primer nivel necesitamos un esquema de relajación. Esta implementado un weighted Jacobi, por su, a mi humilde entender, natural paralelización en CUDA (y que cuando w=1 tenemos un Jacobi), un SOR (sobrerelajación, que con w=1 deriva en un Gauss-Seidel), que suele ser la opción de siempre y un weighted red-black Gauss-Seidel, de cara a paralelizaciones tradicionales (OpenMP, MPI).

En el segundo nivel tenemos los algoritmos de corrección con esquemas V-Cycle, en iterativo, y \mu-Cycle, en recursivo.

Finalmente, tenemos el FMG en el último nivel, tambien en iterativo y en recursivo. En total, 12 combinaciones posibles (en realidad, doce para una dimensión, doce para dos dimensiones y doce para tres dimensiones)

A ver que tal se dan hoy las ejecuciones en dos dimensiones.

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