Chapter3Matrices y lgebralinealSe ha dicho y repetido que Matlab es un programa orientado a clculo matricial. Todo nmero en Matlab, hasta los escalares, son en realidad extensiones de matrices o matrices degeneradas.

Esta caracterstica hace que todas las funciones puedan tomar matrices como argumento de forma natural. Aunque todas las funciones son en el fondo matriciales podemos clasificar unas cuantas cuyo propsito es precisamente manipular y crear matrices.

3.1Rutinas de creacin de matriceseyeMatriz identidadlinspace(base,limit,n)Devuelve un vector fila connelementos equiespaciados entrebaseylimitlogspace(base,limit,n)Devuelve un vector fila connelementos espaciados exponencialmente entre10basey 10limit.meshdom(x,y)Equivalente a meshgrid. Se usaba en versiones anteriores de Matlabmeshgrid(x,y)Crea una malla equiespaciada en dos dimensiones a partir de los vectoresxey. Retorna dos matrices, una con la coordenadaxy el otro con la coordenaday.onesDevuelve una matriz de las dimensiones solicitadas cuyos elementos son todos 1randDevuelve una matriz de las dimensiones solicitadas cuyos elementos son nmeros pseudoaleatorios entre 0 y 1.rand*donde*pueden ser varias letras. Devuelve una matriz de las dimensiones solicitadas cuyos elementos son nmeros pseudoaleatorios que siguen distintas distribucioneszerosDevuelve una matriz cuyos elementos son todos ceros.Las funciones que quizs requieran una pequea explicacin sonlinspace, logspaceymeshgrid. La primera es la extensin de la secuencia (2.5.2) cuando los elementos no son nmeros enteros. La usaremos muy frecuentemente y nos ser muy sencillo dominar su uso:>> linspace(-pi,e,20)ans = Columns 1 through 7: -3.141593 -2.833178 -2.524764 -2.216349 -1.907935 -1.599520 -1.291106 Columns 8 through 14: -0.982692 -0.674277 -0.365863 -0.057448 0.250966 0.559381 0.867795 Columns 15 through 20: 1.176210 1.484624 1.793038 2.101453 2.409867 2.718282logspacees exactamente la misma funcin pero con separacin logartmica:>> logspace(-pi,e,20)ans = Columns 1 through 6: 7.2178e-04 1.4683e-03 2.9870e-03 6.0765e-03 1.2361e-02 2.5147e-02 Columns 7 through 12: 5.1156e-02 1.0407e-01 2.1170e-01 4.3066e-01 8.7610e-01 1.7822e+00 Columns 13 through 18: 3.6256e+00 7.3756e+00 1.5004e+01 3.0523e+01 6.2092e+01 1.2631e+02 Columns 19 and 20: 2.5696e+02 5.2274e+02La funcin meshgrid es un poco ms compleja. Usaremos esta funcin cuando necesitemos una serie bidimensional de puntos equiespaciados que definiremos a partir de dos vectores, las diferencias segn el ejexy segn el ejey. Por ejemplo:>> [xx,yy]=meshgrid(1:2:10,1:3:9)xx = 1 3 5 7 9 1 3 5 7 9 1 3 5 7 9yy = 1 1 1 1 1 4 4 4 4 4 7 7 7 7 7>> 0.5*(xx+yy)ans = 1.0000 2.0000 3.0000 4.0000 5.0000 2.5000 3.5000 4.5000 5.5000 6.5000 4.0000 5.0000 6.0000 7.0000 8.0000diagExtrae e introduce diagonales en matrices.Probablemente esta ltima sea la de mayor sentido fsico puesto que las matrices por bandas se encuentran muy a menudo en problemas de EDP. El ejemplo que adjunta Matlab en la ayuda de la funcin es perfecto para entender su funcionamiento:>> m = 5;>> diag(-m:m) + diag(ones(2*m,1),1) + diag(ones(2*m,1),-1)ans = -5 1 0 0 0 0 0 0 0 0 0 1 -4 1 0 0 0 0 0 0 0 0 0 1 -3 1 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 1 -1 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 1 3 1 0 0 0 0 0 0 0 0 0 1 4 1 0 0 0 0 0 0 0 0 0 1 53.1.1Tipos de argumentos matriciales.Como vimos en el captulo dedicado a los tipos de argumentos podemos modificar el tipo del argumento matricial (su precisin) de modo explcito. Algunas de las funciones anteriores permiten inicializar una matriz incluyendo el atributo de su tipo, es decir, del tamao que ocupar cada elemento en memoria. Por ejemplo, si definimos la matriz unidad como entera de 8 bits:>> i8 = zeros(2,2,'int8')i8 = 0 0 0 0>> i8(1,1)=1023i8 = 127 0 0 0 Por desgracia los operadores an no disponen de la versatilidad de otros lenguajes ms polivalentes como Fortran o Python; la siguiente operacin nos dar error:>> i8*[123.534,324.123;123.431,324.123]??? Error using ==> mtimes Integers can only be combined with...3.2Rutinas de manipulacin de forma.Hay decenas de funciones de manipulacin de forma. No las intentaremos conocer todas porque es un empeo un poco intil. He aqu las ms utilizadasreshape(A,m,n,...)Reordena la matrizApara ajustarse a las dimensionesmn. Esta funcin slo manipula la forma, en el caso de que queramos ampliar su tamao en vez de aadir ceros dar un error.transposeTraspuesta. (Ver los operadores aritmticos para su abreviatura)ctransposeTraspuesta conjugada. (Ver los operadores aritmticos para su abreviatura)Una de las rutinas ms tiles es la siguiente:cat(opt,a,b,...)Concatena matrices dondeoptes la dimensin en la que se acoplarn las matrices y los argumentos siguientes matrices cuyas dimensiones permiten el 'pegado' entre ellas.Un ejemplo del uso decatsera:>> help cat>> A=ones(2,2);>> B=zeros(2,2);>> cat(1,A,B)ans = 1 1 1 1 0 0 0 0>> cat(2,A,B)ans = 1 1 0 0 1 1 0 0En el ejemplo vemos que con la opcinopt=1las matrices se concatenan en sentido longitudinal mientras que conopt=2se concatenan en sentido transversal. Esta funcin sirve como ejemplo de la gran potencia de la notacin matricial. Es muy importante no subestimar la potencia del uso de la creacin y la modificacin de matrices con la notacin tpica como vemos en el siguiente ejemplo. Haremos el clculo anterior sin utilizar ninguna funcin:>> A=ones(2,2);>> B=zeros(2,2);>> [A,B]ans = 1 1 0 0 1 1 0 0>> [A;B]ans = 1 1 1 1 0 0 0 0Que sirva como leccin; nunca hay que subestimar la potencia de la notacin matricial.

3.2.1Creacin directa de matrices de rango mayor que 2.La notacin y la salida por pantalla de Matlab est pensada para el clculo matricial, es decir, mantenernos siempre por debajo del rango 3. La notacin est pensada para separar filas y columnas pero nada ms; no hay ninguna notacin capaz de separar entre planos matriciales; no hay ninguna extensin tensorial.

El mtodo ms intuitivo sera crear la hiper matriz como matriz y luego utilizar la funcinreshapepara cambiarle la forma. Parece obvio que matlab cuente con mtodos un poco menos artesanales.repmat(A,[n,m,...])Crea una matriz por bloques de dimensiones[n,m,...]con copias de la matrizA.Por ejemplo, si queremos crear un cubo de nmeros de lado 5 donde latapay elfondosean matrices unidad y el resto sean ceros podemos hacer lo siguiente:>> cat(3,eye(5),repmat(zeros(5),[1,1,3]),eye(5))ans =ans(:,:,1) = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1ans(:,:,2) = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0ans(:,:,3) = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0ans(:,:,4) = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0ans(:,:,5) = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 3.3Sistemas 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 formaax=bPero una formulacin tan sencilla suele esconder grandes complicaciones. Todos coincidiremos en que aparentemente la solucin del problema anterior es:x=a-1bPero el problema numrico va mucho ms all. La inversa de una matriz slo existe cuando sta es regular y no siempre tendremos tanta suerte. Buscamos una solucin al problema general de la formaa(MN)x(N)=b(M)En el que no siempre existe una inversa paraa. Este anlisis mucho ms general requiere conocimientos que se escapan en parte de los cursos de clculo numrico y plantea una dificultad adicional: el manejo intuitivo de matrices y vectores. Cuando uno se sumerge en el mundo del clculo matricial es bueno que encuentre modos de equivocarse lo menos posible. Es muy comn 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 clculo numrico es utilizar una notacin parecida a la de Einstein en la que aparecen las filas y columnas como coordenadas co variantes y contra variantes:ajixi=bjLa regla mnemotcnica es que los ndices que aparezcan repetidos como subndices y superndicesde izquierda a derecharepresentan un sumatorio y desaparecen en el resultado final. La regla se mantiene en la solucin, obviamente con permutacin de los ndices:xi=(a-1)ijbj

Esta notacin es especialmente til porque ayuda a no equivocarse en los clculos, subndice son filas y superndice son columnas. Por ejemplo para el producto de dos vectores el producto interno (producto escalar) sera:xjyj=kdonde k es un escalar. En cambio si los dos vectores aparecen permutados obtenemos el producto externo de dos vectores en el que se amplan los ndices:yjxi=ajiPara entender perfectamente la siguiente operacin utilizaremos Matlab:>> y=ones(3,1);>> x=2*ones(1,3);>> x*yans = 6>> y*xans = 2 2 2 2 2 2 2 2 2 Todo ello es til cuando tenemos operaciones un poco ms complejas como la siguiente:ykzkakixibjCules sern las dimensiones de la matriz (o vector o escalar) resultante? Tenemos claro desde un principio queaes una matriz,xun vector columna eyyzson vectores fila. Sin hacer ningn clculo sabemos que la solucin tiene la forma:ykzkakixibj=cjdonde la repeticin de ndices en una misma posicin significa operacin escalar y no matricial. Vamos a exponer un ejemplo de la operacin para conseguir una idea ms grfica:octave:27> yy = 4 4 4octave:28> zz = 5 5 5octave:29> aa = 2 2 2 2 2 2 2 2 2octave:30> xx = 3 3 3octave:31> bb = 1 1 1octave:32> (((y.*z)*a)*x)*bans = 1080 1080 1080 Pero centrmonos ms en el problema importante:aijxj=biEs el problema ms importante del anlisis numrico. Casi todos algoritmos se reducen al planteamiento de un sistema de ecuaciones lineales. Los ms interesantes son los que vienen de una ecuacin en derivadas parciales. Las diferencias finitas, volmenes finitos, elementos finitos y mtodos espectrales terminan en la formulacin de un problema de este tipo. El anlisis preciso de estos problemas es una parte esencial de cualquier curso de anlisis numrico y tiene muchos claros y oscuros dependiendo siempre de la forma de la matriza. El siguiente rbol clasifica la mayora de los problemas con su tratamiento: aes cuadrada y regular ano tiene ninguna estructura reconocible La mayora de los elementos deason no nulos. Mtodos directos o iterativos dependiendo del tamao de la matriz. La mayora de los elementos deason nulos. Matrices sparse. Resolucin por mtodos iterativos. atiene una estructura determinada aes tridiagonal. Resolucin por un mtodo directo. No disponible en Matlab. aes hermitiana (a=a). Resolucin por un mtodo directo. No disponible en Matlab. aes subcondicionada o cuadrada singular. Descomposicin en valores singulares o SVD (pseudoinversa) aes sobrecondicionada, es decir, rectangular con ms filas que columnas. Es un problema de mnimos cuadrados o de SVD (pseudoinversa)No disponible en Matlab significa que no tiene ningn mtodo especial para resolver este tipo de problemas. Debemos utilizar el mtodo de resolucin general,el operador resolucin de sistemas de ecuaciones lineales.

3.3.1Matrices cuadradas regulares.Este es el nico caso en el que es rigurosamente cierto que:x=a-1bMatlab proporciona con el operador resolucin de ecuaciones lineales un mtodouniversalpara resolver estos problemas. Aunque Matlab sea capaz de darnos herramientas sencillas para resolver estos sistemas con matrices de gran tamao, en la referencia [3] se mencionan algunas de las dificultades planteadas.

Cuando hemos clasificado los mtodos de resolucin segn las caractersticas de la matriz del sistema hemos diferenciado tres maneras de resolver el problema: los mtodos directos, los iterativos y el uso de las matrices sparse. Los tres mtodos no son exclusivos ya que las matrices sparse no son ms que un mtodo de almacenamiento alternativo para ahorrar memoria; hablaremos de ellas ms adelante. Cundo utilizar un mtodo directo o un mtodo iterativo es una sutileza necesaria puesto que ahorra quebraderos de cabeza importantes.

3.3.2Mtodos directos.Utilizaremos un mtodo directo de resolucin, es decir, el equivalente aINVERTIRla matriza, cuando la matriz est densamente llena de nmeros (no estemos utilizando almacenamientos alternativos como matrices sparse) y el tamao sea moderado. En el caso de las matrices sparse los algoritmos de resolucin algebraicos (descomposicin LU o Cholesky) no pueden competir en la mayora de los casos con la rapidez de los mtodos iterativos. Cuando las matrices son demasiado grandes la resolucin puede ser muy costosa. Esto no significa que Matlab sea incapaz de resolver sistemas de ecuaciones de gran tamao. He aqu un ejemplo que demuestra lo contrario.>> a=rand(1000);>> b=rand(1000,1);>> tic;a\b;tocans = 1.0539Acaba de resolver un sistema de ecuaciones con mil grados de libertad en un segundo. Y la precisin?octave:4> max((a*(a\b))-b)ans = 2.3159e-13Pues del orden de la precisin de los propios clculos. Los algoritmos de resolucin de Matlab son muy sofisticados1y poseen mtodos para mantener un error aceptable en cualquier condicin. El uso de un mtodo directo o uno iterativo se debe principalmente a motivos de velocidad. Utilizaremos un mtodo iterativo slo en el caso de que exista una ganancia de velocidad apreciable o una ventaja en la aplicacin del algoritmo de resolucin como se ve en el ejercicio resuelto8.5.

3.3.3Mtodos iterativos.Utilizaremos un mtodo iterativo cuando no podemosINVERTIRla matriz del sistema en un tiempo razonable. Su objetivo es llegar a la solucin realizando menos operaciones a costa de perder precisin (los mtodos directos son exactos y su precisin depende nicamente de la precisin de la mquina)

Los mtodos iterativos llegan a la solucin mediante la evaluacin directa de la matriz del sistema. Si conseguimos llegar a una solucin evaluando la matriz cien veces el nmero de operaciones ser del orden de 100n, a tener en cuenta si el nmero de elementos es varios rdenes de magnitud mayor al nmero 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 cmo tratar las matrices ``casi vacas''.

Quizs el mtodo iterativo ms utilizado es el delGradiente Conjugado Precondicionadoo PGC que analizaremos tambin en la siguiente seccin. Una de las posibilidades de las funciones que implementan mtodos iterativos es que no necesitan una matriz del sistema, les basta con una funcin que calculeax. Esto es ideal para formular los problemas de un modo alternativo o para utilizar las matricessparsecomo veremos a continuacin.

3.3.4Matrices sparseSe dice que una matriz essparse2cuando el nmero de elementos no nulos es del orden densiendonel nmero de elementos de la matriz. Hay multitud de casos de sistemas cuyas matrices sonsparse, un ejemplo clsico 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 matricessparsees triple. Primero ahorrar la enorme cantidad de memoria desperdiciada almacenando una matriz llena de ceros; segundo, definir mtodos para realizar todas las operaciones disponibles en una matriz pero en el caso de almacenamiento sparse y por ltimo que los errores de resolucin no sean los equivalentes a resolver una matriz connelementos sino una conn. Resumiendo, aprovechar la forma de la matriz para obtener todas las propiedades beneficiosas que se pueda.

3.3.4.1Anlisis de matricesEs muy importante saber qu forma tiene nuestra matriz. Si es diagonal superior, tridiagonal o sparse influir en gran medida el algoritmo de resolucin del problema. Para analizar la naturaleza de las matrices contamos con funciones de gran utilidad:issparseSi la matriz argumento es sparse retornar verdadero.issquareSi la matriz es cuadrada retorna el nmero de elementos por lado, si no lo es retorna falso.spyEs un mtodo 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 seccin3.3.4.3.3.3.4.2Almacenamiento de matrices sparseExisten bastantes mtodos de almacenamiento de matricessparse.Los ms conocidos sonCompressed Row Storage,Compressed Column Storage,Block Compressed Row Storage,Compressed Diagonal Storage,Jagged Diagonal StorageySkyline Storage.

Probablemente el ms 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=10000-20

390003

078700

308750

0809913

04002-1

los vectores resultado del CRS son: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 ms 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 funcinsparse,sparse(M)Almacena la matrizMcomo 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) -> 3 Nos devolver la informacin de cmo ha almacenado la matriz3, nunca los vectores que est utilizando para almacenarla. Lo ms 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*bans = 19 12 7 3.3.4.3Creacin de matrices sparseEl uso de la funcinsparseno es nunca una buena opcin para crear una matriz sparse porque no tiene absolutamente ninguna lgica. El motivo es claro, si la existencia de las matrices sparse es ahorrar memoria y aumentar la velocidad no tiene ningn sentido que tengamos que crear antes la matriz original y luego tener que liberar memoria.

Entonces pensamos que sera muy cmodo contar con las mismas funciones de creacin de matrices convencionales para el caso sparse:speyeCrea la matriz unidad en formato sparsesponesCrea una matriz sparse con unos en los elementos que deseemossprandCrea una matriz sparse de nmeros aleatorios en posiciones aleatorias o siguiendo un determinado patrn.Por ejemplo, si queremos una matriz sparse de nmeros aleatorios con una densidad de 0.01 y de dimensiones 100100:>> randsparse=sprand(100,100,0.01);Para ver cmo se reparten aleatoriamente los nmeros en la matriz utilizaremos la funcinspyque ya conocemos:>> spy(randsparse)Y como resultado obtenemos la figura3.1:

Figure 3.1:Elementos no nulos de una matriz sparse creada consprand

spdiagsExtrae e introduce diagonales en matrices sparse.Al igual que su funcin homlogadiagtiene un gran sentido fsico. De hecho tiene mucha mayor utilidad en la funcin de creacinspdiagsquediagporque las matrices en banda entran en la definicin de matrices sparse. La nica consideracin a tener en cuenta es que las matrices en banda se resuelven rpidamente de manera exacta mediante una factorizacin LU diseada para matrices sparse. Utilizaremos los mtodos iterativos con matrices que no cumplan ningn patrn evidente o sencillo.

3.3.4.4Manipulacin y operaciones con matrices sparsesphcatConcatena horizontalmente dos matrices sparse de dimensiones compatibles.spvcatConcatena verticalmente dos matrices sparse de dimensiones compatibles.spabsCalcula el valor absoluto de los elementos complejos de una matriz sparse sin cambiar el carcter de la misma.sprealCalcula la parte real de los elementos complejos de una matriz sparse sin cambiar el carcter de la misma.spimagCalcula la parte imaginaria de los elementos de una matriz sparse sin cambiar el carcter de la misma.spinvCalcula la matiz inversa de una matriz sparse.Sin embargo no es recomendable que utilicemosspinvpara resolver sistemas de ecuaciones lineales. Como en el caso de sistemas de ecuaciones en matrices convencionales es recomendable que optemos por el operadorPARA INVERTIRla matriz del sistema.

3.3.5Matrices tridiagonales (Octave)Son uno de los tipos de matrices ms buscadas; tienen un gran sentido fsico (provienen de algunas ecuaciones en diferencias) y son muy adecuadas para el clculo numrico. Son un caso particular de matrices sparse pero al ser tan sencillas pueden resolverse directamente sin realizar demasiados clculos. La funcin de resolucin de sistemas de ecuaciones con una matriz tridiagonal estrisolve.

3.3.6Matrices no regulares.Una de las caractersticas ms importantes del operador es que tambin sirve en el caso de tener un sistema con exceso o defecto de ecuaciones. El operador ya no calcula la matriz inversa sino que busca una solucin cuyo error cumple la condicin de la norma mnima o busca el subespacio solucin. En este caso se dice que en vez de buscar la inversa se calcula una pseudoinversa que resuelve la ecuacin. El tratamiento de este tipo de problemas suele omitirse en la mayora de cursos de lgebra en ingeniera, sin embargo tiene una gran utilidad en clculo numrico y es una herramienta relativamente importante. Creo que no es en vano si intentamos entender con un poco ms de profundidad el tratamiento de estas aplicaciones lineales sin entrar en formalismos.

3.3.6.1Singular Value Decomposition (SVD)Se demuestra que cualquier matriz de aplicacin lineal acepta una descomposicin del tipo:A=UwVdondeUyVson matrices ortogonales ywuna matriz diagonal que contiene los valores singulares. Las matrices ortogonales, o ms concretamente las matrices de una transformacin ortogonal, tienen la propiedad de mantener los ngulos en la transformacin. Las matrices de giro son un ejemplo de transformacin ortogonal. La matrizw, al ser diagonal, no produce ningn tipo de rotacin, simplemente expande o contrae determinadas direcciones. La demostracin de la existencia de esta descomposicin va mucho ms all de los objetivos de este libro. Se utiliza la notacinVpara esta matriz debido a que para cumplir la definicin de matiz de transformacin ortogonal sus columnas deben ser vectores ortogonales; es slo un formalismo porque se demuestra tambin que la inversa de una transformacin ortogonal es tambin una transformacin ortogonal y es equivalente a su traspuesta. Ayuda escribir la anterior ecuacin especificando las dimensiones:A(NM)=U(NM)w(MM)V(MM)Vemos tambin que las dimensiones son correctas, aplicando queVV-1:A(NM)V(MM)=U(NM)w(MM)

Ahora, un repaso de lgebra de parte de [14]. Si nuestra matrizAes la matriz de una aplicacin linealf:STpodemos definir dos subespacios enS; el Ncleo, subespacio que forman los vectores que tienen al vector nulo como imagen y el subespacio Imagen que es el subespaciof(S)deT. Por ejemplo, seanSyTlos siguientes espacios reales: Sson las funciones polinmicasp(x)=a+bx+cx2+dx3 Tel espacio de funciones deRenRConsideremos la aplicacin lineal def:STdefinida porf(p(x))=p(x), la derivada segn la variable independiente; o sea:f(a+bx+cx2+dx3)=b+2cx+3dx2. La imagen defes el subespacio deTformado por los polinomios de grado menor o igual que 2. El ncleo defes el subespacio deTque forman los polinomios constantes.

Si volvemos a la forma de la descomposicin propuesta por la SVD,A(NM)=U(NM)w(MM)V(MM)

Tomemos la condicin deNMen el que tenemos ms ecuaciones que incgnitas.

Para obtener las tres matrices de la descomposicin Matlab proporciona la siguiente funcinsvdDescomposicin en valores singulares. Si se pide nicamente un valor de salida retorna un vector con valores singulares. Con tres valores de salida retorna la descomposicin matricial completa.Para ejemplificar el uso de la ecuacin vamos a intentar descomponer en valores singulares la ya clsica matriz...octave:3> a=[1,2,3;4,5,6;7,8,9];octave:4> [U,w,V]=svd(a)U = -0.21484 0.88723 0.40825 -0.52059 0.24964 -0.81650 -0.82634 -0.38794 0.40825w = 16.84810 0.00000 0.00000 0.00000 1.06837 0.00000 0.00000 0.00000 0.00000V = -0.479671 -0.776691 -0.408248 -0.572368 -0.075686 0.816497 -0.665064 0.625318 -0.408248Efectivamente, uno de los valores singulares es nulo debido a que la matriz es singular. Para analizar ms en profundidad las matrices, Matlab proporciona adems las siguientes funciones:rankCalcula el rango de la matriz, es decir, el nmero de valores singulares no nulos o la dimensin de la imagen de la aplicacin.nullCalcula la dimensin del ncleo de una aplicacin lineal expresada por una matriz, esto es, el nmero de valores singulares nulos.3.3.6.2Problemas con defecto o exceso de ecuaciones.Si lo nico que nos interesa es resolverax=b... Qu calcula exactamente el operador cuandoN> x = [-1 3 4]x =-1 3 4>> x = [-1;3;4]x =-134

Para introducir unamatriz, se define por filas separndolas por intro o punto y coma. Por ejemplo, el comando siguiente define la matriz A de dimensin 2X3>> A = [3 4 1; 5 -3 0]A =3 4 15 -3 01. Vectores y sus operacionesLasoperacionesde adicin, diferencia yproductopor escalares, se realizan mediante una sintaxis simple.Ejemplo 1.1Consideremos los vectores x = < -1,3 -2>, y = y el escalar 2.>> x = [-1 3 -2]x =-1 3 -2>> y = [0 -2 -1]y =0 -2 -1

>> z = x + yz =-1 1 -3>> w = x - yw =-1 5 -1

>> u = 2*xu =-2 6 -4

Geomtricamente, MATLAB ilustra mediante laleydel paralelogramo la suma de vectores. Para ello, se utiliza un programa escrito por T. A. Bryan [4], asequible en la pgina. Aqu, simplemente, se utilizarel lenguajedegrficaspara puntualizar la operacin.Ejemplo 1.2Sean los vectores x = , y = y x + y = >> u = [0 1 1 0 0];>> v = [0 2 5 3 0];>> w = [0 0 2 2 0];>> plot3(u,v,w)Figura 1.1

Para determinar la norma de un vector se utiliza el comando siguiente:Ejemplo 1.3>> norm(x)ans =2.2361Para el producto punto, el producto vectorial y el ngulo entre vectores, se procede de la forma siguiente:Ejemplo 1.4>> dot(x,y)ans =6>> cross([1,2,0],[0,3,2])ans =4 -2 3>> theta = acos(dot(x,y)/(norm(x)*norm(y)))theta =0.7314>> 360*theta/(2*pi)ans =41.9088

Para obtener el rea del paralelogramo cuyos lados adyacentes son los vectores x y y.>> area = norm(x)*norm(y)*(sin(theta))^2area =3.5970

Ejercicio propuesto 1.1 Obtener un vector unitario en ladireccindel vector x. Hallar la proyeccin del vector x sobre el vector y.2.SistemasdeecuacioneslinealesPara resolver unsistemade ecuaciones lineales AX = B de n ecuaciones con m (n puede ser igual a m) incgnitas, se introduce la matriz A del sistema y el vector columna B de los trminos independientes, no es preciso considerar el vector columna X de las incgnitas (x1, x2, x3).Ejemplo 2.1Consideremos el siguiente sistema:

>> A = [0 3 -4; 6 -3 -4; 6 -9 4; 1 1 1]A =0 3 -46 -3 -46 -9 41 1 1>> B = [0; 0; 0; 1]B =0001

>> X = A\BX =0.36360.36360.2727

Observacin. El operador matricial de MATLAB "\" divisin izquierda equivale a la solucin de sistemas lineales mediante X = inv(A)*B. este operador es ms poderoso de lo que parece, puesto que, suministra la solucin aunque la matriz A no tenga inversa.MATLAB, proporciona la solucin grfica de un sistema de ecuaciones lineales, mediante loscomandossiguientes:>> [x,y] = meshgrid(-4:0.5:5);>> z = 3*y/4;>> surf(x,y,z)Se obtiene el plano de la figura 2.1

Figura 2.1>> hold on>> z = (6*x - 3*y)/4;>> surf(x,y,z)Se obtiene los planos de la figura 2.2 interceptados en una lnea recta

Figura 2.2

>> z = (-6*x + 9*y)/4;>> surf(x,y,z)Se obtiene los planos de la figura 2.3 interceptados en una lnea rectaFigura 2.3>> z = 1- x - y;>> surf(x,y,z)Se obtienen los cuatro planos de la figura 2.4 interceptados en el punto (4/11. 4/11, 3/11) de IR3

Figura 2.4

Ejemplo 2.2Consideremos, ahora un sistema lineal incompatible.x + z = 1x y + 3z = -3x + y z = 1>> A = [1 0 1; 1 -1 3; 1 1 -1]A =1 0 11 -1 31 1 -1>> B = [1; -3; 1]B =1-31

>> X = A\BWarning: Matrix is singular to working precision.X =InfInfInf

Ejemplo 2.3Consideremos, ahora un sistema lineal compatible indeterminadox1 - x2 = 4x1 3x2 2x3 = -6x1 + 2x2 3x3 = 1>> A = [1 -1 0; 1 3 -2; 4 2 -3]A =1 -1 01 3 -24 2 -3>> B = [4; -6; 1]B =4-61

>> X = A\BWarning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 9.251859e-018.X =3.5000-0.50004.0000

El paquete desoftwareMATLAB permite la solucin de estos sistemas utilizando elMtodode Gauss-Jordan. Veamos la solucin del ejemplo 2.3.>> A = [1 -1 0 4; 1 3 -2 -6; 4 2 -3 1]% A es la matriz ampliadaA =1 -1 0 41 3 -2 -64 2 -3 1>> A(2,:) = A(2,:)-A(1,:)A =1 -1 0 40 4 -2 -104 2 -3 1

>> A(3,:) = A(3,:) - 4*A(1,:)A =1 -1 0 40 4 -2 -100 6 -3 -15>> A(2,:) = A(2,:)/-2A =1 -1 0 40 -2 1 50 6 -3 -15

>> A(3,:) = A(3,:) + 3*A(2,:)A =1 -1 0 40 -2 1 50 0 0 0Luegox1 = 4 + x2 y x3 = 5 + 2x2Infinitassoluciones!Por ejemplo, para x2 = - 0,5x1 = 3,5 y x3 = 4

Nota. Elcarcter% indica comienzo de un comentario, es ignorado por MATLAB.Ejercicio propuesto 2.1Resolver los sistemas de ecuaciones dados. Ver ejemplo 2.1.a)b)3. MatricesPara multiplicar una matriz por un vector y para el producto de matrices, se procede de la forma siguiente:Ejemplo 3.1SeanA =y x =Se introducen la matriz y el vector como se indic.Luego el comando>> B = A*xB =1313-2Ejemplo 3.2Consideremos las matricesA = [5 -2 0; 1 -1 4; 3 -3 2; 0 -5 -1] y B = [-2 3 4 6 7 -1; 0 9 -2 1 0 3; 1 -1 0 3 2 7]>> C = A*BC =-10 -3 24 28 35 -112 -10 6 17 15 24-4 -20 18 21 25 2-1 -44 10 -8 -2 -22>> D = B*A??? Error using ==> *Inner matrix dimensions must agree.

Ejemplo 3.3Para la matriz A = [1/2 -1/3 1/6; 1/4 -2 -1/7; 1/5 1/3 0]. Hallar:Transpuesta>> A'ans =0.5000 0.2500 0.2000-0.3333 -2.0000 0.33330.1667 -0.1429 0Inversa>> inv(A)ans =0.4181 0.4878 3.3449-0.2509 -0.2927 0.99304.2439 -2.0488 -8.0488Determinante>> det(A)ans =0.1139

Rango>> rank(A)ans =3Traza>> trace(A)ans =-1.5000Valores singulares>> svd(A)ans =2.07950.53670.1020

A3>> A^3ans =0.2651 -1.1341 -0.04600.8310 -7.4944 -0.6200-0.0945 1.4008 0.1353log(A)>> logm(A)ans =-0.8432 + 0.3607i 0.3205 + 0.2443i0.5082 - 0.9163i-0.2189 - 0.2087i 0.7465 + 3.1599i0.2834 - 0.0688i0.6527 - 1.1504i -0.6112 + 0.1011i-2.0758 + 2.7626i

sqrt(A)>> sqrtm(M)ans =0.6528 + 0.0198i -0.0574 + 0.1899i 0.2151 - 0.0822i0.0490 - 0.1420i -0.0043 + 1.4284i 0.0161 + 0.0841i0.2701 - 0.0979i -0.0237 - 0.1954i 0.0890 + 0.2721i

Para determinar eA., se utilizan las variantes mediante autovalores, aproximantes de Pad, desarrollos deTaylory condicin de la matriz A [6].>> expm(A)ans =1.6408 -0.1807 0.23100.1400 0.1146 -0.04670.2862 0.1195 1.0092>> expm1(A)ans =1.6408 -0.1807 0.23100.1400 0.1146 -0.04670.2862 0.1195 1.0092

>> expm2(A)ans =1.6408 -0.1807 0.23100.1400 0.1146 -0.04670.2862 0.1195 1.0092>> expm3(A)ans =1.6408 -0.1807 0.23100.1400 0.1146 -0.04670.2862 0.1195 1.0092

Como es posible observar, la matriz exponencial coincide en todos losmtodosempleados.Diagonalizar la matriz A calculando la matriz de paso V.>> [V,J] = jordan(A)

V =0.8852 -0.0169 + 0.0000i 0.1317 + 0.0000i0.0664 + 0.0000i -0.1131 - 0.0000i 0.0467 - 0.0000i0.3662 0.0212 + 0.0000i -0.3874 - 0.0000i

J =0.5439 - 0.0000i 0 00 -1.9358 00 0 -0.1082

Para matrices especiales MatLab ofrece ciertos comandos [6].Ejemplo 3.4(a) Obtener las formas de Smith y de Hermite de la inversa de la matriz de Hilbert de orden 2 en la variable t. Adems, determinar las matrices de paso.>> maple('with(linalg):H:=inverse(hilbert(2,t))');>> pretty(simple(sym(maple('H'))))[ 2 ][ -(-3 + t) (-2 + t) (-3 + t) (-2 + t) (-4 + t)][ ][ 2 ][(-3 + t) (-2 + t) (-4 + t) -(-3 + t) (-4 + t) ]

>> maple('B:=smith(H,t,U,V);U:=eval(U);V:=eval(V)');>> pretty(simple(sym(maple('B'))))[-3 + t 0 ][ ][ 2 ][ 0 (-2 + t) (t - 7 t + 12)]

>> pretty(simple(sym(maple('U'))))[ -1 -1 ][ ][ 2 2][10 - 13/2 t + t 9 - 13/2 t + t ]

>> pretty(simple(sym(maple('V'))))[-7/2 + t -4 + t][ ][-3/2 + t -2 + t]

>> maple('HM:=hermite(H,t,Q);Q:=evalm(Q)');>> pretty(simple(sym(maple('HM'))))[ 2 ][6 + t - 5 t 0 ][ ][ 2 ][ 0 t - 7 t + 12]

>> pretty(simple(sym(maple('Q'))))[3 - t 2 - t][ ][4 - t 3 - t]

(b) Determinar los autovalores de la matriz de Wilkinson de orden 6, de la matriz mgica de orden 6 y de la matriz de Rosser.>> [eig(wilkinson(8)),eig(rosser),eig(magic(8))]ans =1.0e+003 *-0.0010 -1.0200 0.26000.0002 -0.0000 0.05180.0011 0.0001 -0.05180.0017 1.0000 0.00000.0026 1.0000 0.0000 + 0.0000i0.0028 1.0199 0.0000 - 0.0000i0.0042 1.0200 -0.0000 + 0.0000i0.0043 1.0200 -0.0000 - 0.0000i

(c) Determinar la matriz y el determinante jacobianos de la transformacin:x = eusen(v), y = eucos(v).>> pretty(sym(maple('jacobian(vector([exp(u)*sin(v),exp(u)*cos(v)]),[u,v])')))[exp(u) sin(v) exp(u) cos(v) ][ ][exp(u) cos(v) -exp(u) sin(v)]

pretty(simple(sym(maple('det('')'))))2-exp(u)

4. Dependencia eindependencialineal, base y dimensinPara determinar si unafamiliao conjunto de vectores de un espacio vectorial es L.I. o L.D., se calcula el determinante de la matriz conformada por las componentes de los vectores.Ejemplo 4.1Dados los espacios vectoriales IR3, IR4 y IR5, determinar si los siguientesconjuntosde vectores son L.I. o L.D.(a) {(3,2,-1), (1,0,0), (0,1,2)}>> A = [3 2 -1; 1 0 0; 0 1 2]A =3 2 -11 0 00 1 2>> det(A)ans =-5Puesto que, el determinante es diferente de cero, la familia de vectores es L.I.(b) {(2,1,-3,4), (-1,3,2,1),(-5,1,8,-7), (3,2,1,-1)}>> B = [2 1 -3 4;-1 3 2 1;-5 1 8 -7;3 2 1 -1]B =2 1 -3 4-1 3 2 1-5 1 8 -73 2 1 -1>> det(B)ans =0Puesto que, el determinante es igual a cero, la familia de vectores es L.D.

(c) {(1,2,-1,5,3), (1,-1,4,-2,0), (1,1,-1,3,12), (0,4,3,1,-1)}>> C = [1,2,-1,5,3;1,-1,4,-2,0;1,1,-1,3,12;0,4,3,1,-1]C =1 2 -1 5 31 -1 4 -2 01 1 -1 3 120 4 3 1 -1>> rank(C)ans =4Puesto que, los 4 vectores pertenecen al espacio IR5, no es posible aplicar el determinante, sin embargo son L.I., ya que el rango de la matriz, que conforman, es 4. Si es menor que 4 es L.D.

Ejemplo 4.2Dadala familiade vectores en el espacio IR5{(1,2,-1,5,3), (1,-1,4,-2,0), (1,1,-1,3,1,2), (0,4,3,1,-1)}Determinar la dimensin de la variedad lineal engendrada por ella y una base.La dimensin de la variedad lineal engendrada por una familia de vectores es el rango de la matriz conformada por los vectores. Puesto que, el rango es 4, la dimensin es 4.Para la base se utiliza el siguiente comando>> pretty(sym(maple('basis({vector([1,2,-1,5,3]),vector([1,-1,4,-2,0]),vector([1,1,-1,3,12]),vector([0,4,3,1,-1])})')))[1 2 -1 5 3][ ][1 -1 4 -2 0][ ][0 4 3 1 -1][ ][1 1 -1 3 12]Ejemplo 4.3Dada la familia de vectores en el espacio IR3{(1,2,-1), (1,-1,4), (1,1,-1)}Determinar si conforman una base de IR3, en caso positivo, obtener las componentes del vector v = (2,4,1) en dicha base.>> det([1,2,-1;1,-1,4;1,1,-1])ans =5Constituyen una base, puesto que el determinante de la matriz es diferente de cero. La familia es L.I.Las componentes de v, se obtienen con el comando>> inv([1 1 1;2 -1 1;-1 4 -1])*[2 4 1]'ans =3.20000.6000-1.8000Ejemplo 4.4Consideremos las bases del espacio vectorial IR3A = {(1,0,0), (0,1,0), (0,0,1)} y B = {(2,1,0), (1,0,-1), (-1,1,1)}Hallar la matriz delcambiode base de A a B. Adems, calcular las componentes del vector (1,2,3) en base A, en la base B.>> A = [1,0,0;0,1,0;0,0,1];>> B = [2,1,0;1,0,-1;-1,1,1];

>> C = inv(B')*A'C =0.5000 0 0.5000-0.5000 1.0000 -1.5000-0.5000 1.0000 -0.5000>> sym(C)ans =[ 1/2, 0, 1/2][ -1/2, 1, -3/2][ -1/2, 1, -1/2]

Una vez obtenida la matriz del cambio de base, se determinan las componentes del vector (1,2,3) en la base B.>> sym(inv(B')*A'*[1,2,3]')ans =[ 2][ -3][ 0]

5. Transformaciones linealesEjemplo 5.1Consideremos la transformacin lineal, cuya matriz est conformada por la familia de vectores{(0,-3,-1,3), (-3,3,-1,0), (-2,2,1,3)}Hallar una base de su ncleo. Adems, determinar laimagendel vector (1,2,3,4) por intermedio de la transformacin lineal.>> A = [0,-3,-1,3;-3,3,-1,0;-2,2,1,3]A =0 -3 -1 3-3 3 -1 0-2 2 1 3>> null(A)ans =0.64480.4690-0.52760.2931

>> maple('T:=x->multiply(array([[0,-3,-1,3],[-3,3,-1,0],[-2,2,1,3]]),x)')ans =T := proc (x) options operator, arrow; multiply(array([[0, -3, -1, 3], [-3, 3, -1, 0], [-2, 2, 1, 3]]),x) end>> pretty(sym(maple('T([1,2,3,4])')))[3 0 17]

Ejemplo 5.2Sea la transformacin lineal f entre los espacios vectoriales V y W, tal quef(e1) = u1 u2, f(e2) = u2 u3 y f(e3) = u3 u4dondeA = {e1, e2, e3} es una base de VB = {u1, u2, u3} es una base de WDeterminar la matriz asociada a la transformacin lineal f. Adems, hallar la imagen en W del vector (1,2,3) de V a travs de la transformacin.>> M = [1 0 0;-1 1 0;0 -1 1;0 0 1]M =1 0 0-1 1 00 -1 10 0 1>>maple('T1:=x->multiply(array([[1,0,0],[-1,1,0],[0,-1,1],[0,0,1]]),x)');>> pretty(sym(maple('T1([1,2,3])')))[1 1 1 3]

Ejemplo 5.3Sea la transformacin lineal f entre los espacios vectoriales V y W, dondef(x,y,z) = (x + y, y + z, x + z), (x,y,z) es un punto cualquiera de V.Hallar la matriz asociada a las transformaciones lineales f, f3 y ef.>> maple('T:=(x,y,z)->[x+y,y+z,x+z]');

>> [maple('T(1,0,0)'),maple('T(0,1,0)')maple('T(0,0,1)')]ans =[1, 0, 1][1, 1, 0][0, 1, 1]Ntese, que para hallar la matriz de f, es necesario considerar los vectores transformados por f de la base cannica de IR3.Tambin, es posible hallarla con este comando>>[maple(T(1,0,0),maple(T(0,1,0)),maple(T(0,0,1))]A =1 1 00 1 11 0 1

>> A^3 % la matriz asociada a f3 es A3.ans =2 3 33 2 33 3 2

>> expm(A) % la matriz asociada a ef es eAans =3.1751 2.8321 1.38191.3819 3.1751 2.83212.8321 1.3819 3.1751

6. OrtonormalizacinLa tcnica bsica utilizada en MATLAB para determinar la descomposicin QR se sustenta en transformaciones de Householder. La sintaxis es sumamente simple [4].Ejemplo 6.1Consideremos la matriz con los vectorescomo columnas.Luego la matriz viene dada por>> A = [1 2 0; 0 -1 0; -1 0 -2; 0 1 1]A =1 2 00 -1 0-1 0 -20 1 1>> [Q R] = qr(A,0)

Q =-0.7071 0.5000 -0.45230 -0.5000 -0.15080.7071 0.5000 -0.45230 0.5000 0.7538R =-1.4142 -1.4142 -1.41420 2.0000 -0.50000 0 1.6583

7. DiagonalizacinEjemplo 7.1Diagonalizar la matriz A>> A = [-1 2 -1;3 0 1;-3 -2 -3]A =-1 2 -13 0 1-3 -2 -3Se determina los eigenvalores de A>> eig(A)ans =-42-2

Ahora se determinan los eigenvectores v1, v2, v3 asociados con los eigenvalores 4, 2 , -2>> v1=null(A+4*eye(3));>> v2=null(A-2*eye(3));>> v3=null(A+2*eye(3));>> S=[v1 v2 v3]S =0.5774 0.5774 0.5774-0.5774 0.5774 -0.57740.5774 -0.5774 -0.5774>> inv(S)*A*Sans =-4.0000 0 -0.0000-0.0000 2.0000 -0.0000-0.0000 0.0000 -2.0000A es diagonizable, puesto que existe una matriz S, tal que=

Leer ms:http://www.monografias.com/trabajos16/algebra-lineal/algebra-lineal.shtml#ixzz3bRaJ4JqW