Introducción informal a Matlab y Octave - Matrices y algebra lineal (II)
Tutorial creado por Guillem Borrell i Nogueras. Extraido de: http://torroja.dmt.upm.es/%7Eguillem/matlab/
05 de Noviembre de 2006
Encriptación
9 - Matrices y algebra lineal (II)
3.3 Sistemas de ecuaciones lineales.
Un sistema de ecuaciones lineales es:| a11x1+a12x2+a13x3+⋯+a1NxN | = | b1 |
| a21x1+a22x2+a23x3+⋯+a2NxN | = | b2 |
| a21x1+a22x2+a23x3+⋯+a2NxN | = | b3 |
| ⋮ | = | ⋮ |
| aM1x1+aM2x2+aM3x3+⋯+aMNxN | = | bM |
siempre puede expresarse de la forma
ax=b
Pero una formulación tan sencilla suele esconder grandes complicaciones. Todos coincidiremos en que aparentemente la solución del problema anterior es:
x=a-1b
Pero el problema numérico va mucho más allá. La inversa de una matriz sólo existe cuando ésta es regular y no siempre tendremos tanta suerte. Buscamos una solución al problema general de la forma
a(MנN)x(N)=b(M)
En el que no siempre existe una inversa para a. Este análisis mucho más general requiere conocimientos que se escapan en parte de los cursos de cálculo numérico y plantea una dificultad adicional: el manejo intuitivo de matrices y vectores. Cuando uno se sumerge en el mundo del cálculo matricial es bueno que encuentre modos de equivocarse lo menos posible. Es muy común que multiplicaciones de matrices incompatibles acaben con un algoritmo brillante. No siempre trabajaremos con vectores columna post multiplicando matrices cuadradas. Un truco bastante utilizado en el cálculo numérico es utilizar una notación parecida a la de Einstein en la que aparecen las filas y columnas como coordenadas co variantes y contra variantes:
ajixi=bj
La regla mnemotécnica es que los índices que aparezcan repetidos como subíndices y superíndices de izquierda a derecha representan un sumatorio y desaparecen en el resultado final. La regla se mantiene en la solución, obviamente con permutación de los índices:
xi=(a-1)ijbj
Esta notación es especialmente útil porque ayuda a no equivocarse en los cálculos, subíndice son filas y superíndice son columnas. Por ejemplo para el producto de dos vectores el producto interno (producto escalar) sería:
xjyj=k
donde k es un escalar. En cambio si los dos vectores aparecen permutados obtenemos el producto externo de dos vectores en el que se amplían los índices:
yjxi=aji
Para entender perfectamente la siguiente operación utilizaremos Matlab: >> y=ones(3,1); >> x=2*ones(1,3); >> x*y ans = 6 >> y*x ans = 2 2 2 2 2 2 2 2 2Todo ello es útil cuando tenemos operaciones un poco más complejas como la siguiente:
ykzkakixibj
¿Cuáles serán las dimensiones de la matriz (o vector o escalar) resultante? Tenemos claro desde un principio que a es una matriz, x un vector columna e y y z son vectores fila. Sin hacer ningún cálculo sabemos que la solución tiene la forma:
ykzkakixibj=cj
donde la repetición de índices en una misma posición significa operación escalar y no matricial. Vamos a exponer un ejemplo de la operación para conseguir una idea más gráfica: octave:27> y y = 4 4 4 octave:28> z z = 5 5 5 octave:29> a a = 2 2 2 2 2 2 2 2 2 octave:30> x x = 3 3 3 octave:31> b b = 1 1 1 octave:32> (((y.*z)*a)*x)*b ans = 1080 1080 1080Pero centrémonos más en el problema importante:
aijxj=bi
Es el problema más importante del análisis numérico. Casi todos algoritmos se reducen al planteamiento de un sistema de ecuaciones lineales. Los más interesantes son los que vienen de una ecuación en derivadas parciales. Las diferencias finitas, volúmenes finitos, elementos finitos y métodos espectrales terminan en la formulación de un problema de este tipo. El análisis preciso de estos problemas es una parte esencial de cualquier curso de análisis numérico y tiene muchos claros y oscuros dependiendo siempre de la forma de la matriz a. El siguiente árbol clasifica la mayoría de los problemas con su tratamiento:
- a es cuadrada y regular
- a no tiene ninguna estructura reconocible
- La mayoría de los elementos de a son no nulos. Métodos directos o iterativos dependiendo del tamaño de la matriz.
- La mayoría de los elementos de a son nulos. Matrices sparse. Resolución por métodos iterativos.
- a tiene una estructura determinada
- a es tridiagonal. Resolución por un método directo. No disponible en Matlab.
- a es hermitiana (a=a⊤′). Resolución por un método directo. No disponible en Matlab.
- a no tiene ninguna estructura reconocible
- a es subcondicionada o cuadrada singular. Descomposición en valores singulares o SVD (pseudoinversa)
- a es sobrecondicionada, es decir, rectangular con más filas que columnas. Es un problema de mínimos cuadrados o de SVD (pseudoinversa)
3.3.1 Matrices cuadradas regulares.
Este es el único caso en el que es rigurosamente cierto que:x=a-1b
Matlab proporciona con el operador resolución de ecuaciones lineales un método universal para resolver estos problemas. Aunque Matlab sea capaz de darnos herramientas sencillas para resolver estos sistemas con matrices de gran tamaño, en la referencia [3] se mencionan algunas de las dificultades planteadas.Cuando hemos clasificado los métodos de resolución según las características de la matriz del sistema hemos diferenciado tres maneras de resolver el problema: los métodos directos, los iterativos y el uso de las matrices sparse. Los tres métodos no son exclusivos ya que las matrices sparse no son más que un método de almacenamiento alternativo para ahorrar memoria; hablaremos de ellas más adelante. Cuándo utilizar un método directo o un método iterativo es una sutileza necesaria puesto que ahorra quebraderos de cabeza importantes.
3.3.2 Métodos directos.
Utilizaremos un método directo de resolución, es decir, el equivalente a invertir la matriz a, cuando la matriz esté densamente llena de números (no estemos utilizando almacenamientos alternativos como matrices sparse) y el tamaño sea moderado. En el caso de las matrices sparse los algoritmos de resolución algebraicos (descomposición LU o Cholesky) no pueden competir en la mayoría de los casos con la rapidez de los métodos iterativos. Cuando las matrices son demasiado grandes la resolución puede ser muy costosa. Esto no significa que Matlab sea incapaz de resolver sistemas de ecuaciones de gran tamaño. He aquí un ejemplo que demuestra lo contrario.>> a=rand(1000); >> b=rand(1000,1); >> tic;a\b;toc ans = 1.0539Acaba de resolver un sistema de ecuaciones con mil grados de libertad en un segundo. ¿Y la precisión?
octave:4> max((a*(a\b))-b) ans = 2.3159e-13Pues del orden de la precisión de los propios cálculos. Los algoritmos de resolución de Matlab son muy sofisticados1 y poseen métodos para mantener un error aceptable en cualquier condición. El uso de un método directo o uno iterativo se debe principalmente a motivos de velocidad. Utilizaremos un método iterativo sólo en el caso de que exista una ganancia de velocidad apreciable o una ventaja en la aplicación del algoritmo de resolución como se ve en el ejercicio resuelto 8.5.
3.3.3 Métodos iterativos.
Utilizaremos un método iterativo cuando no podemos invertir la matriz del sistema en un tiempo razonable. Su objetivo es llegar a la solución realizando menos operaciones a costa de perder precisión (los métodos directos son exactos y su precisión depende únicamente de la precisión de la máquina)Los métodos iterativos llegan a la solución mediante la evaluación directa de la matriz del sistema. Si conseguimos llegar a una solución evaluando la matriz cien veces el número de operaciones será del orden de 100n, a tener en cuenta si el número de elementos es varios órdenes de magnitud mayor al número de iteraciones necesarias. Por supuesto en el caso que casi todos los elementos de la matriz sean distintos de cero. En el siguiente apartado veremos cómo tratar las matrices ``casi vacías''.
Quizás el método iterativo más utilizado es el del Gradiente Conjugado Precondicionado o PGC que analizaremos también en la siguiente sección. Una de las posibilidades de las funciones que implementan métodos iterativos es que no necesitan una matriz del sistema, les basta con una función que calcule ax. Esto es ideal para formular los problemas de un modo alternativo o para utilizar las matrices sparse como veremos a continuación.
3.3.4 Matrices sparse
Se dice que una matriz es sparse2 cuando el número de elementos no nulos es del orden de n siendo n el número de elementos de la matriz. Hay multitud de casos de sistemas cuyas matrices son sparse, un ejemplo clásico son las matrices generadas por esquemas de elementos finitos tanto en el caso de mallas estructuradas como el de no estructuradas.La finalidad del uso de las matrices sparse es triple. Primero ahorrar la enorme cantidad de memoria desperdiciada almacenando una matriz llena de ceros; segundo, definir métodos para realizar todas las operaciones disponibles en una matriz pero en el caso de almacenamiento sparse y por último que los errores de resolución no sean los equivalentes a resolver una matriz con n elementos sino una con n. Resumiendo, aprovechar la forma de la matriz para obtener todas las propiedades beneficiosas que se pueda.
3.3.4.1 Análisis de matrices
Es muy importante saber qué forma tiene nuestra matriz. Si es diagonal superior, tridiagonal o sparse influirá en gran medida el algoritmo de resolución del problema. Para analizar la naturaleza de las matrices contamos con funciones de gran utilidad:- issparse
- Si la matriz argumento es sparse retornará verdadero.
- issquare
- Si la matriz es cuadrada retorna el número de elementos por lado, si no lo es retorna falso.
- spy
- Es un método muy utilizado cuando ya sabemos que nuestra matriz va a ser sparse. Pinta un diagrama cuyos puntos son los elementos de la matriz que son distintos de cero. Un ejemplo de su uso se encuentra en la sección 3.3.4.3.
3.3.4.2 Almacenamiento de matrices sparse
Existen bastantes métodos de almacenamiento de matrices sparse. Los más conocidos son Compressed Row Storage, Compressed Column Storage, Block Compressed Row Storage, Compressed Diagonal Storage, Jagged Diagonal Storage y Skyline Storage.Probablemente el más sencillo sea el primero. El algoritmo descompone la matriz en tres vectores, uno con los elementos de la matriz distintos de cero, otro con el índice de columna donde se encuentra y un tercero el índice de los vectores anteriores que inicia una nueva fila. Si se descompone la matriz :
| A= | 漂R> 缂R> 缂R> 缂R> 缂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?> ?> ?> ?> ?NT> |
val=(10,-2,3,9,3,7,8,7,3,…,13,4,2,-1)
col_ ind=(1,5,1,2,6,2,3,4,1,…,5,6,2,5,6)
row_ ptr=(1,3,6,9,13,17,20)
El tipo de almacenamiento es oculto en Matlab principalmente por motivos de sencillez.Pero lo más importante del uso de matrices sparse es que ni sepamos que hemos cambiado el tipo de almacenamiento. Por ejemplo, si a partir de una matriz normal pasamos a almacenamiento sparse gracias a la función sparse,
- sparse(M)
- Almacena la matriz M como matriz sparse
>> M=[1 0 2 3;0 0 0 3;3 2 0 0] M = 1 0 2 3 0 0 0 3 3 2 0 0 >> SpM=sparse(M) SpM = Compressed Column Sparse (rows=3, cols=4, nnz=6) (1 , 1) -> 1 (3 , 1) -> 3 (3 , 2) -> 2 (1 , 3) -> 2 (1 , 4) -> 3 (2 , 4) -> 3Nos devolverá la información de cómo ha almacenado la matriz3, nunca los vectores que está utilizando para almacenarla. Lo más importante es que esta nueva matriz sparse puede operar con matrices de todo tipo ya sean sparse o convencionales transparentemente.
>> b=[1;2;3;4]
b =
1
2
3
4
>> SpM*b
ans =
19
12
7
Valora este capítulo:
Autor y licencia de 'Introducción informal a Matlab y Octave - Matrices y algebra lineal (II)'
|
Opiniona sobre 'Introducción informal a Matlab y Octave - Matrices y algebra lineal (II)' (13)
Tu nombre debe tener tres caracteres como mínimo.
Es necesario que te des de alta con una cuenta de correo válida.
Es necesario que te des de alta con una cuenta de correo válida.
El contenido del título de tu opinión debe tener tres caracteres como mínimo.
Es obligatorio que selecciones una valoración del recurso.
El contenido del comentario de tu opinión debe tener tres caracteres como mínimo.
Opina sobre este tutorial |
Wikis relacionados con 'Introducción informal a Matlab y Octave - Matrices y algebra lineal (II)'
Las fotografias de flores (flora en general) quizas sean las que mejor se dejan enmarcar....
Más »
Este trabajo ha tenido en cuenta los supuestos teóricos analizados en el artículo “Competencias: Un...
Más »
En la edición anterior, se explicó las bases de Netfilter/IPTables. En esta segunda entrega, se...
Más »
Género gramatical y sexo no son, como muchos ingenuos o espontáneos usuarios de la lengua...
Más »
En la primera parte se introdujo un estudio transtextual de la obra, además de discutir...
Más »
