Labsesion 6 - Eigenvalor y Eigenvector PDF

Preview

Summary

Eigenvalor y Eigenvector...

Description

Pr´actica de Matlab. Autovalores y autovectores ´ Algebra Lineal. Curso 2010-11

En esta pr´actica vamos a ver algunas funciones que nos permiten calcular los autovalores y autovectores de una matriz. Adem´as, expondremos un m´etodo num´erico para calcularlos y lo programaremos en r Matlab .

1.

Algunos comandos ´ utiles La orden

>> poly(a) aplicada a una matriz a devuelve el polinomio caracter´ıstico. Y hemos visto que la orden roots permite r calcular sus ra´ıces (los autovalores). Pero Matlab dispone de ´ordenes m´as directas. Escribiendo >> eig(a) encontramos los autovalores de a, y con una orden del tipo >> [x,d]=eig(a) los autovectores y autovalores. En general, la orden anterior devuelve una matriz x y una matriz diagonal d, con la condici´on de que x*d=a*x. Pero esta matriz x no tiene por qu´e ser inversible, es decir, realiza esta computaci´on aunque a no sea diagonalizable. Esto se debe al car´acter num´erico de los c´alculos r que hace Matlab y a la alta sensibilidad de la condici´on de diagonalizaci´on de una matriz por peque˜ nos cambios en los coeficientes. Ejercicio 1 Escribe una matriz A aleatoria, del tama˜ no que quieras. Calcula una matriz P y una matriz diagonal D tales que A = P DP −1 , y compru´ebalo. r

Matlab ofrece la posibilidad de realizar c´alculos exactos de manera simb´olica. Para ellos, podemos definir la matriz empleando la orden >> b=sym([2 1 0 4 7; 0 2 4 1 -3; 0 0 1 8 -2; 0 0 0 1 5;0 0 0 0 3]) b = [2 [ [0 [ [0 [ [0 [ [0

1

0

4

2

4

1

0

1

8

0

0

1

0

0

0 1

7] ] -3] ] -2] ] 5] ] 3]

Observemos que la matriz est´ a expresada de manera diferente a la habitual. Ahora, la orden >> [x,d]=eig(b) x = [105/2 [ [71/2 [ [ 9 [ [ 5/2 [ [ 1

1

4] ] -4] ] 1] ] 0] ] 0]

0 0 0 0

d = [3 [ [0 [ [0 [ [0 [ [0

0

0

0

2

0

0

0

2

0

0

0

1

0

0

0

0] ] 0] ] 0] ] 0] ] 1]

proporciona una lista m´axima de autovectores independientes.

2.

M´ etodo de las potencias

Describimos en esta secci´on un m´etodo num´erico para calcular el autovalor de m´odulo m´aximo de una matriz. Supongamos que A es una matriz diagonalizable de tama˜ no n, con autovalores verificando |λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λn | , es decir, con un u ´nico autovalor de m´odulo m´aximo, y con multiplicidad algebraica 1. Esta situaci´on es gen´erica si se toma una matriz arbitraria de tipo num´erico. Sea x0 un vector “casi” arbitrario, y xk = Ak x0 . En las condiciones del enunciado, puede demostrarse que la direcci´on de xk se aproxima a la de un autovector correspondiente al autovalor λ1 , a condici´on de que, si {v1 , v2 , . . . , vn } es una base de autovectores, vi asociado a λi , y escribimos x0 = α1 v1 + α2 v2 + · · · + αn vn , se tenga que α1 6= 0. Asimismo, xk+1 ∼ λ1 xk . Esto es la base del m´etodo de las potencias que pasaremos a implementar. Este vector xk puede tender a cero o a infinito, con lo que hay que normalizarlo para que esto no ocurra. Para ello hay varias posibilidades: una de ellas es reducir al valor 1 la componente de x0 de mayor r valor absoluto, dividiendo por ella. Nos pueden ser u ´tiles los comandos abs y max de Matlab . Otra es dividir por la norma de x0 . La orden norm es entonces necesaria. Ejercicio 2 r

Busca en la ayuda de Matlab informaci´on sobre los comandos abs, max, y norm aplicados a r vectores o a matrices. Escribe un programa en Matlab que tome un vector y lo divida por una componente cuyo m´odulo sea m´aximo.

2

Vamos a programar el m´etodo de las potencias siguiendo el siguiente esquema, que corresponde a la primera de las t´ecnicas de normalizaci´on que hemos comentado: Dados A, x0 , numiter Desde k = 1 hasta numiter xk = Axk−1 Definimos αk como una componente de xk de m´ odulo m´aximo. xk = xk /αk % Normalizamos el vector xk . Fin del bucle. Devuelve xnumiter, αnumiter. Es un esquema muy simple, donde numiter es el n´ umero de iteraciones que efectuamos, que es fijo en r este programa. En c´odigo de Matlab tenemos lo siguiente: function [x,l]=potencia(a,v) % Metodo de las potencias aplicado a una matriz a, a partir de un vector v % Devuelve un autovector x y el autovalor l de modulo maximo numiter=100; x=v(:); %Esta orden hace de v un vector columna, y lo llama x for k=1:numiter x=a*x; [m,p]=max(abs(x)); l=x(p); x=1/l*x; end return Como hemos dicho, una forma alternativa de programar es dividir en cada paso el vector por su norma. Es decir, a partir del vector yk , de norma 1, calculamos xk+1 = Ayk . Para valores de k grandes se tiene que xk+1 ∼ λ1 yk , y como yk tiene norma 1, λ1 resulta ser aproximadamente el producto escalar r de xk+1 y yk . En c´odigo de Matlab : function [y,l]=potencia2(a,v) % Metodo de las potencias % En cada paso se divide el vector que se obtiene por su norma % El numero de iteraciones es fijo numiter=100; x=v(:); y=x/norm(x); for k=1:numiter x=a*y; y=x/norm(x); end l=x.’*y; Fijar el n´ umero de iteraciones no es habitualmente una buena estrategia. Es m´ as eficaz poner un l´ımite a dicho n´ umero, y fijar una precisi´on. Cuando se considere que esta precisi´on se ha alcanzado, detenemos el proceso. Una posibilidad es ir calculando los valores del autovalor en cada paso, y detenernos cuando r este valor apenas se modifique. Modificamos el c´odigo de Matlab : 3

function [y,l]=potencia3(a,v) % Metodo de las potencias % En cada paso se divide el vector que se obtiene por su norma % numiter es el numero m\’{a}ximo de iteraciones % eps es la precision fijada % % Si la diferencia entre dos valores consecutivos calculados de l % es menor que eps, el programa se detiene numiter=100; eps=1e-10; x=v(:); y=x/norm(x); l=0; for k=1:numiter x=a*y; y=x/norm(x); lc=x.’*y; if abs(l-lc)...