%% Obtención de la matriz L en la factorización FA=LU
% En este "Live Script" de Matlab se justifica la forma en que se obtiene la matriz L en la factorización FA=LU de una matriz A en un caso general para n=6, suponiendo que solamente hay dos cambios de filas: en el pivoteo de la columna 2 (2 \leftrightarrow 4) y la columna 4 (4 \leftrightarrow 6).
% Ultano Kindelán Bustelo 9/2022
%% Matrices elementales
% Antes de empezar la explicación conviene recordar qué son las matrices elementales. Hay tres tipos de matrices elementales, cada una asociada a una de las tres operaciones elementales que se pueden realizar sobre las filas de una matriz:
%% Matriz de cambio de filas: I_{i\leftrightarrow j}\in M_{n}(\Re)
% Es el resultado de intercambiar las filas i, j en la matriz unidad.
%% Ejemplo en M_6(\Re)

n=6;
I6=eye(n)
I35=I6;
I35([3,5],:)=I35([5,3],:)

% Si se premultiplica una matriz A\in M_{6,m}(\Re) por I_{3\leftrightarrow5} se obtiene la misma matriz A con las filas 3 y 5 intercambiadas:

A = sym('a',[n,n])
A35=I35*A

% En general si se premultiplica una matriz A\in M_{n,m}(\Re) por I_{i \leftrightarrow j} se obtiene la misma matriz A con las filas i y j intercambiadas.
% Observación: si se intercambian las columnas i y j de la matriz unidad se obtiene la misma matriz I_{i \leftrightarrow j}:

Ic35=I6;
Ic35(:,[3,5])=Ic35(:,[5,3])

% y si se postmultiplica una matriz A\in M_{m,n}(\Re) por I_{i \leftrightarrow j} se obtiene la misma matriz A con las columnas i y j intercambiadas:

Ac35=A*I35

%% Matriz que multiplica una fila por un escalar: I_{\alpha i} \in M_n(\Re)
% Es el resultado de multiplicar la fila i por \alpha en la matriz unidad.
%% Ejemplo en M_6(\Re)

I2_4=I6;
I2_4(4,:)=2*I2_4(4,:)

% Si se premultiplica una matriz A\in M_{6,m}(\Re) por I_{24} se obtiene la misma matriz A con las fila 4 multiplicada por 2:

A2_4 = I2_4*A

% En general si se premultiplica una matriz A\in M_{n,m}(\Re) por I_{\alpha i} se obtiene la misma matriz A con las fila i multiplicada por \alpha.
% Observación: si se multiplica la columna i de la matriz unidad por \alpha se obtiene la misma matriz I_{\alpha  i}:

Ic2_4=I6;
Ic2_4(:,4)=2*Ic2_4(:,4)

% y si se postmultiplica una matriz A\in M_{m,n}(\Re) por I_{\alpha i} se obtiene la misma matriz A con la columna i multiplicada por \alpha:

Ac2_4=A*Ic2_4

%% Matriz que reemplaza una fila por ella misma más otra multiplicada por un escalar: I_{j+\alpha i}\in M_n(\Re)
%  Es la matriz obtenida al reemplazar la fila j de la matriz unidad por ella misma más la fila i multiplicada por \alpha.
%% Ejemplo en M_6(\Re)

I4_32=I6;
I4_32(4,:)=I4_32(4,:)+3*I4_32(2,:)

% Si se premultiplica una matriz A\in M_{6,m}(\Re) por I_{4+3\times2} se obtiene la misma matriz A con la fila 4 reemplazada por ella misma más la 2 multiplicada por 3:

A4_32=I4_32*A

% En general si se premultiplica una matriz A\in M_{n,m}(\Re) por I_{j+\alpha i} se obtiene la misma matriz A con la fila j sustituida por ella misma más la i multiplicada por \alpha.
% Observación: si en la matriz unidad se reemplaza la columna j de la matriz unidad por ella misma más la columna i multiplicada por \alpha que denominaremos I_{j+\alpha i}^c que verifica I_{j+\alpha i}^c=I_{i+\alpha j}= (I_{j+\alpha i})^t:

I4c_32=I6;
I4c_32(:,4)=I4c_32(:,4)+3*I4c_32(:,2)

% y se postmultiplica una matriz A\in M_{m,n}(\Re) por I_{j+\alpha i}^c se obtiene la misma matriz A con la columna j reemplazada por ella misma más la i multiplicada por \alpha:

Ac4_32=A*I4c_32

%%  En resumen
% Las operaciones elementales sobre las filas de A son equivalentes a premultiplicar A por matrices elementales obtenidas haciendo las mismas operaciones sobre las filas de I_n.
% Las operaciones elementales sobre las columnas de A son equivalentes a postmultiplicar A por matrices elementales obtenidas haciendo las mismas operaciones sobre las columnas de I_n.
%% Factorización FA=LU
%% Notación
% A= Matriz cuadrada que se quiere factorizar y que en este ejemplo supondremos de dimensión 6.
% U= Matriz triangular superior obtenida al realizar la triangularización gaussiana con pivote parcial en A.
% F= Matriz que almacena los cambios de filas (al buscar el pivote de mayor valor absoluto) realizados en el proceso de triangularización gaussiana de A.
% L= Matriz triangular inferior de la que se quiere hallar su expresión en función de los multiplicadores m_{ik} utilizados en la triangularización.
% P^k= Matriz de pivoteo en la etapa k.
% F^k= Matriz de cambio de filas en la etapa k.
%% Triangularización gausssiana con pivote parcial
% Se van a expresar las operaciones elementales realizadas sobre la matriz A mediante premultiplicaciones por matrices elementales. La matriz L se obtendrá a partir del producto de estas matrices elementales. La matriz A es una matriz genérica de dimensión 6 y se supondrá que para diagonalizarla con pivote parcial hay que intercambiar las filas 2 y 4 en la etapa k=2 y las filas 4 y 6 en la etapa k=4.
%% k=1
% En esta etapa no hay cambio de filas por lo tanto F^1=I_6.

F1=I6;

% La matriz de pivoteo P^1 será de la forma

[P1] = gen_mat_P_l(6,1,1,1,1,1)

% en donde m_{i1}=-a_{i1}/a_{11}
% Por lo tanto la expresión matricial de la matriz A transformada después de la primera etapa será:
% P^1F^1A.
%% k=2
% En esta etapa se intercambian las filas 2 y 4  F^2=I_{2 \leftrightarrow4}.

F2=I6;
F2([2,4],:)=F2([4,2],:)

% La matriz de pivoteo P^2 será de la forma

[P2] = gen_mat_P_l(6,2,1,1,1,1)

% en donde m_{i2}=-a_{i2}/a_{22}.
% La expresión matricial después de la segunda etapa será:
% P^2F^2P^1F^1A.
%% k=3
% En esta etapa no hay cambio de filas:  F^3=I_6

F3=I6;

% La matriz de pivoteo P^3 será de la forma

[P3] = gen_mat_P_l(6,3,1,1,1,1)

% en donde m_{i3}=-a_{i3}/a_{33}.
% La expresión matricial después de la tercera etapa será:
% P^3F^3P^2F^2P^1F^1A.
%% k=4
% En esta etapa se intercambian las filas 4 y 6  F^4=I_{4 \leftrightarrow 6}.

F4=I6;
F4([4,6],:)=F4([6,4],:)

% La matriz de pivoteo P^4 será de la forma

[P4] = gen_mat_P_l(6,4,1,1,1,1)

% en donde m_{i4}=-a_{i4}/a_{44}.
% La expresión matricial después de la cuarta etapa será:
% P^4F^4P^3F^3P^2F^2P^1F^1A.
%% k=5
% En esta etapa no hay cambio de filas:  F^5=I_6

F5=I6;

% La matriz de pivoteo P^5 será de la forma

[P5] = gen_mat_P_l(6,5,1,1,1,1)

% en donde m_{i5}=-a_{i5}/a_{55}.
% La quinta etapa es la última y, por tanto, al acabarla se obtiene la matriz escalonada (triangular superior) U:
% P^5F^5P^4F^4P^3F^3P^2F^2P^1F^1A=U.
% Para obtener L, en el término de la izquierda, se van a mover todas las F^khacia la derecha utilizando la siguiente propiedad
% Propiedad 1: Suponiendo que F^k=I_{i\leftrightarrow j}, siempre se verifica
% F^kP^l=P^l_{i\leftrightarrow j}F^k \quad \text{para } k>l,
% en donde P^l_{i\leftrightarrow j} es igual a P^l con los elementos i, j de la columna l intercambiados.
% No se demuestra la propiedad pero se comprobará que es cierta, para los casos en los que la aplicamos:
% Se desplaza F^2 hacia la derecha:
% F^2P^1=P^1_{2 \leftrightarrow 4}F^2

[~,P1_24] = gen_mat_P_l(6,1,2,4,1,1)
F2P1=F2*P1
P1_24F2=P1_24*F2

% Se desplaza F^3hacia la derecha:
% F^3P^2=P^2F^3
% F^3P^1=P^1F^3
% Se desplaza F^4hacia la derecha:
% F^4P^3=P_{4 \leftrightarrow 6}^3F^4
% F^4P^2=P_{4 \leftrightarrow 6}^2F^4
% F^4P_{2 \leftrightarrow 4}^1=P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1F^4

[~,P3_46] = gen_mat_P_l(6,3,4,6,1,1)
[~,P2_46] = gen_mat_P_l(6,2,4,6,1,1)
[~,~,P1_24_46] = gen_mat_P_l(6,1,2,4,4,6)
F4P3=F4*P3
P3_46F4=P3_46*F4
F4P2=F4*P2
P2_46F4=P2_46*F4
F4P1_24=F4*P1_24
P1_24_46F4=P1_24_46*F4

% Se desplaza F^5hacia la derecha:
% F^5P^4=P^4F^5
% F^5P^3=P^1F^3
% F^5P^2=P^2F^5
% F^5P^1=P^1F^5
% Entonces podemos expresar el proceso de triangularización como
% P^5P^4P_{4 \leftrightarrow 6}^3P^2_{4 \leftrightarrow 6}P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1F^5F^4F^3F^2F^1A=U
% .
% LLamando F=F^5F^4F^3F^2F^1se tiene
% P^5P^4P_{4 \leftrightarrow 6}^3P^2_{4 \leftrightarrow 6}P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1FA=U
% con lo que
% FA=\left(P^5P^4P_{4 \leftrightarrow 6}^3P^2_{4 \leftrightarrow 6}P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1\right)^{-1}U
% y, por lo tanto,
% L=\left(P^5P^4P_{4 \leftrightarrow 6}^3P^2_{4 \leftrightarrow 6}P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1\right)^{-1}=\left(P_{\begin{array}{c} 2 \leftrightarrow 4\\ 4 \leftrightarrow 6 \end{array}}^1\right)^{-1}(P^2_{4 \leftrightarrow 6})^{-1}(P_{4 \leftrightarrow 6}^3)^{-1}(P^4)^{-1}(P^5)^{-1}
% .

P1_24_46_inv=inv(P1_24_46)
P2_46_inv=inv(P2_46)
P3_46_inv=inv(P3_46)
P4_inv=inv(P4)
P5_inv=inv(P5)
L=P1_24_46_inv*P2_46_inv*P3_46_inv*P4_inv*P5_inv

%

F=F5*F4*F3*F2*F1