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

Anuncios