autovalor.m

function [K, X, Y] = autovalor(F, positivos)
%autovalor
%   [K, X, Y] = autovalor(F)
%
%   Calcula los autovalores y los autovectores para el problema
%   polin�mico inverso
%      X (f0 + f1/k + f2/k^2 + f3/k^3 + ...) = 0
%      (f0 + f1/k + f2/k^2 + f3/k^3 + ...) Y = 0
%      1 + X (-f1/k^2 - 2 f2/k^3 - 3 f3/k^4 + ...) Y = 0
%
%   [K, X, Y] = autovalor(F, positivos)
%      F           recetas en tiempo discreto
%      positivos   si no es [] devuelve s�lo los autovalores positivos
%      K           autovalores soluci�n
%      X           autovectores por la izquierda
%      Y           autovectores por la derecha
%
%   Los autovectores por la izquierda son vectores-fila y los autovectores
%   por la derecha son vectores-columna. 

%   necesita forma1.m, formaab.m

if nargin < 2, positivos = []; end
cero = 1e-6;                                        %define la precisi�n del cero
[filas, columnas, ancho] = size(F);                 %determina el tama�o del problema
[A, B] = formaab(F, 1);                             %convierte F(k) a la forma -A + B/k
[X1, K1] = eig(B', A');                             %resuelve el problema de autovalores y autovectores por la izquierda
[Y2, K2] = eig(B, A);                               %resuelve el problema de autovalores y autovectores por la derecha
K = []; X = []; Y = [];                             %crea las variables para almacenar las soluciones
for c1 = 1:size(K1, 2)                              %para cada autovalor del problema por la izquierda
    Knn = K1(c1, c1);                                   %extrae el autovalor
    if or(isempty(positivos), and(abs(imag(Knn))<cero, real(Knn)>cero)) %si positivos es [] o si Knn es real y Knn > 0
        if and(Knn == Knn, and(abs(Knn)<Inf, abs(Knn)>cero))%si Knn >< NaN y Knn >< +-Inf y Knn >< 0
            DF = forma1(F, Knn);                                %calcula F'(k)
            Xnn = X1(1:filas, c1)';                             %extrae el X de los procesos que se inician
            s = sum(Xnn); if abs(s) > cero, Xnn = Xnn ./ s; end %normaliza X para que sum(X) = 1
            for c2 = 1:size(K2, 2)                              %para cada autovalor por la derecha
                if abs(Knn - K2(c2, c2)) < cero                     %si el autovalor por la izquierda y la derecha coinciden
                    Ynn = Y2(1:columnas, c2);                            %extrae el Y de las materias
                    q = -Xnn*DF*Ynn; if abs(q)>cero, Ynn = Ynn ./ q; end %normaliza Y para que 1 + X F'(k) Y = 0
                    K = [K, Knn]; X = [X; Xnn]; Y = [Y, Ynn];            %almacena la soluci�n
                end
            end
        end
    end 
end