jueves, 29 de agosto de 2013

Resolución de sistemas compatibles determinados con wxMaxima (II): métodos de factorización

En la primera entrada de la serie expliqué algunos métodos directos para resolver, con wxMaxima, sistemas de ecuaciones lineales que son compatibles determinados.

Vamos a ver ahora, utilizando el mismo ejemplo, otro tipo de métodos para resolver este tipo de sistemas: los métodos de factorización. Podemos resumir dichos métodos diciendo que el objetivo es, dado un sistema Ax=b, encontrar una factorización de A en dos matrices cuyas inversas sean "más fáciles" de calcular (en realidad "más fácil" significa menos coste en términos computacionales). Si A = B · C, entonces el sistema Ax=b nos queda BCx=b y, por consiguiente, podríamos descomponer el sistema Ax=b en dos sistemas By=b, Cx=y (o aplicar el método de la inversa para despejar el vector de soluciones x = C^(-1) · (B^(-1) · b) ).

A continuación explico cómo ejecutar en wxMaxima 3 de estos métodos de factorización: LU, Cholesky y QR.

Recordad que estamos resolviendo el sistema:
10x+4y+4z+3t+5s=1
  4x+8y+2z+2t+5s=1
  4x+2y+6z+3t+4s=1
  3x+2y+3z+8t+5s=1
  5x+5y+4z+5t+6s=1
Y ya comprobamos en la entrada anterior que es un sistema compatible determinado.


Método de factorización LU:

lu_backsub(lu_factor(A),transpose(b));


Método de factorización de Cholesky:

CH:cholesky(A);
 
Generalmente, si A no cumple las condiciones necesarias para poder aplicar el método de Cholesky (que sea Hermítica y definida positiva), wxMaxima devuelve un mensaje de error. En cualquier caso, no está de más comprobar que la matriz devuelta por wxMaxima es la correcta:
is(CH.transpose(CH)=A);
true

invert(transpose(CH)).(invert(CH).b);
ratsimp(%);



Método de factorización QR:

En este caso es necesario cargar primero el paquete "lapack" para poder utilizar la función del procedimiento de factorización QR.
load(lapack);

[Q,R]:dgeqrf(A);
invert(R).(transpose(Q).b);

 

invert(R).(transpose(Q).b);



En la próxima entrada de esta serie trataré algunos de los métodos de resolución denominados iterativos: método de Jacobi, método de Gauss-Siedel y método de relajación.

No hay comentarios: