Block Cyclic Reduction

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 $latex x$, tres puntos para las $latex y$ y siete puntos para las $latex z$ (necesitamos que sean números impares).

Para la dimensión $latex z$ tenemos:

$latex 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)$

$latex -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 $latex y$ no queda:

$latex 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 $latex x$ obtenemos:

$latex 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 $latex n$ dimensiones $latex x_1, x_2, ldots, x_n$ con $latex t_1, t_2, ldots, t_n$ puntos en cada dimensión, ya que nos queda $latex A_i$ una matriz tridiagonal por bloques con $latex t_i times t_i$ bloques donde cada uno está formado por matrices $latex A_{i-1}$ en la diagonal, $latex -I_{i-1}$ arriba y debajo de la diagonal y $latex 0$ en el resto, todas de dimensión $latex t_{i-1} times t_{i-1}$ con $latex i=1,ldots, n$. En $latex A_n$ tenemos el stencil $latex (-1,2n,-1)$.

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

$latex -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 $latex ^1$ hace referencia a las $latex x$ y no corresponde a potencias… Para resolver las $latex k$ impares simplemente:

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

Y resulta que definiendo $latex B_2:=(A_2)^2 – 2 I_2$, nos queda una matriz $latex B_1$ de la misma estructura que $latex 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 $latex N=2^{k+1}-1$).

Llegados a este punto, si estamos con $latex v^n$, resolvemos como queramos, pero si estamos en otra variable $latex v^i$ con $latex ineq 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 ($latex A_i rightarrow B_i rightarrow ldots rightarrow F_i$, $latex F$ por poner algo :-)); el otro cuando, llegados a una sola ecuación aplicando la recursividad anterior, nos quedan mas dimensiones por resolver ($latex B_{j+1}, B_{j+2}, ldots $).

Tags: ,

Reply

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *


¡IMPORTANTE! Responde a la pregunta: ¿Cuál es el valor de 7 2 ?
 
FireStats icon Powered by FireStats