perturbacion.m

function [sk, sX, sY, sm, U, iU] = perturbacion(F, G, k, X, Y, problema, grado)
%perturbacion
%   [sk, sX, sY] = perturbacion(F, G, k, X, Y)
%
%   Series perturbativas para el problema de autovalores, para VN o para TE. 
%
%   [sk, sX, sY, sm, U, iU] = perturbacion(F, G, k, X, Y, problema, grado)
%      F         recetas en tiempo discreto con los procesos de producci�n
%      G         recetas en tiempo discreto con la perturbaci�n
%      k         autovalor o factor
%      X         autovector por la izquierda o intensidades
%      Y         autovector por la derecha o valores
%      problema  para autovalores 0, para VN 1 (por defecto), y para TE 2
%      grado     grado de la serie (por defecto 100)
%      sk        serie para el factor
%      sX        serie para las intensidades
%      sY        serie para los valores
%      sm        serie para el multiplicador de la normalizaci�n
%      U         matriz de coeficientes del sistema lineal
%      iU        inversa de la matriz de coeficientes
%
%   F y G deben tener la misma dimensi�n. Las series usan el formato de los
%   polinomios en Matlab-Octave. Podemos aproximar la soluci�n de VN (o
%   para los otros problemas) para las recetas F + h*G para un escalar h
%   dado con el comando
%      formapolyval(sk, h)     para obtener el factor
%      formapolyval(sX, h)'    para las intensidades
%      formapolyval(sY, h)     para los valores
%   siempre que estemos dentro del radio de convergencia y aplicaci�n de 
%   las series. 

if nargin < 6, problema = []; end
if nargin < 7, grado = []; end

if isempty(problema), problema = 1; end              %por defecto calcula las series para VN
if isempty(grado), grado = 100; end                  %por defecto el grado de las series es 100

[filas, columnas, ancho] = size(F);                  %establece el tama�o del problema
if problema == 0, fbas = (X == X); cbas = (Y == Y); else fbas = (X ~= 0); cbas = (Y ~= 0); end %determina los procesos y materias b�sicos
if problema == 1, q = 1; else q = 0; end             %si q es 1, k(h) - X(h)(A + h M)Y(h) = 0; si q es 0, 1 - X(h)(A + h M)Y(h) = 0
if sum(fbas) ~= sum(cbas), disp(['Atenci�n, la soluci�n no es cuadrada']); end %avisa si la soluci�n no es cuadrada
X0 = X(fbas);                                        %extrae las intensidades de los procesos b�sicos
for c2 = 1:ancho-2                                   %para el resto de pasos temporales
    X0 = [X0, X(fbas)./k^c2];                            %completa las intensidades que simbolizan la evoluci�n de los procesos
end                                                  %
if ancho > 2                                         %si hay m�s de dos pasos temporales
    Yr = F(fbas,cbas,ancho)*Y(cbas)./k; Y0 = Yr;         %calcula el Yr de la materia que simboliza el proceso en el �ltimo paso temporal y lo agrega a los valores
    for c2 = ancho-1:-1:3                                %para el resto de pasos temporales
        Yr = (F(fbas,cbas,c2)*Y(cbas) + Yr)./k; Y0 = [Yr;Y0]; %calcula Yr de la materia que simboliza el proceso para ese paso temporal y lo agrega a los valores
    end                                                  %
    Y0 = [Y(cbas); Y0];                                  %agrega el Y de las materias b�sicas a los valores
else                                                 %si s�lo hay dos pasos temporales
    Y0 = Y(cbas);                                        %los valores son los de las materias b�sicas
end                                                  %
[A, B] = formaab(F(fbas,cbas,:));                    %establece A y B para los procesos y materias b�sicos
[M, N] = formaab(G(fbas,cbas,:));                    %establece M y N para los procesos y materias b�sicos

[sk, sXn, sYn, sm, U, iU] = serieperturbacion(A, B, M, N, k, X0, Y0, q, grado, sum(fbas)); %calcula las series

sX = zeros(filas, grado+1);                          %define la variable para almacenar la serie de las X
sX(fbas, :) = sXn(1:sum(fbas), :);                   %extrae las series de las intensidades de los procesos b�sicos y las a�ade
sY = zeros(columnas, grado+1);                       %define la variable para almacenar la serie de los Y
sY(cbas, :) = sYn(1:sum(cbas), :);                   %extrae las series de los valores de las materias b�sicas y las a�ade




function [k, X, Y, m, U, iU] = serieperturbacion(A, B, M, N, k0, X0, Y0, q, grado, filasf)
%Resuelve sucesivamente el sistema lineal U dZn = Rn para cada grado n
%      |   0        -A*Y0     B-k0*A     1| |Xn|   |  |
%      |(-A*Y0)'      q       -X0*A      0| |kn| = |Rn|
%      |(B-k0*A)'  (-X0*A)'      0       0| |Yn|   |  |
%      |   1          0          0       0| |mn|   |  |
[filas, columnas] = size(A);                                       %establece el tama�o del problema
normaliza = [ones(filasf,1); zeros(filas-filasf,1)];               %establece que la suma de las X de los procesos que se inician sea cte
U = [zeros(filas), -A*Y0   , B-k0*A           , normaliza        ; %establece la matriz del sistema lineal
     (-A*Y0)'    ,   q     , -X0*A            , 0                ; %
     (B-k0*A)'   , (-X0*A)', zeros(columnas)  , zeros(columnas,1); %
     normaliza'  , 0       , zeros(1,columnas), 0               ]; %
iU = inv(U);                                                       %calcula U^-1
r = grado+1;                                                       %establece el n�mero de elementos de las series
X=zeros(filas,r); k=zeros(1,r); Y=zeros(columnas,r); m=zeros(1,r); %define dX, dk, dY, dm para almacenar las soluciones
X(:,r-(0))=X0'; k(r-(0))=k0; Y(:,r-(0))=Y0; m(:,r-(0))=0;          %establece los t�rminos de grado 0
for cn = 1:grado                                                   %para cada grado n
    ci = 0;                                                            %construye Rn con los t�rminos para i = 0
    Rn = [M*Y(:,r-(cn-ci-1))*k(r-(ci)) - N*Y(:,r-(cn-1))     ;         %
          X(:,r-(ci))'*M*Y(:,r-(cn-ci-1))                    ;         %
          (k(r-(ci))*X(:,r-(cn-ci-1))'*M - X(:,r-(cn-1))'*N)';         %
          0                                                 ];         %
    for ci = 1:cn-1                                                    %a�ade a Rn los t�rminos para i = 1:n-1
        Rn = Rn + [k(r-(ci))*A*Y(:,r-(cn-ci)) + M*Y(:,r-(cn-ci-1))*k(r-(ci))      ; 
                   X(:,r-(ci))'*A*Y(:,r-(cn-ci)) + X(:,r-(ci))'*M*Y(:,r-(cn-ci-1)); 
                   (k(r-(ci))*X(:,r-(cn-ci))'*A + k(r-(ci))*X(:,r-(cn-ci-1))'*M)' ;
                   0                                                             ];
    end                                                                %
    Zn = iU * Rn;                                                      %resuelve el sistema lineal con Zn = U^-1 Rn
    X(:,r-(cn)) = Zn(1:filas);                                         %extrae de la soluci�n Xn
    k(r-(cn)) = Zn(filas+1);                                           %extrae de la soluci�n kn
    Y(:,r-(cn)) = Zn(filas+2:filas+columnas+1);                        %extrae de la soluci�n Yn
    m(r-(cn)) = Zn(filas+columnas+2);                                  %extrae de la soluci�n mn
end                                                                %