Eliminacion Gaussiana CON pivoteo

Aqui les presento una nueva version del codigo de Matlab que ya antes habia publicado que implementa eliminacion Gaussina para resolver sistemas lineales del tipo Ax = b. Esta nueva version no es mas rapida que la anterior, pero el error relativo en los resultados es mucho menor que en la version anterior. El error que se obtiene con este metodo se debe en gran medida a que los resultados fueron obtenidos mediante calculos en representacion finita.

%**********************************************************
%Intercambia dos elementos de un arreglo en una dimension
%Recibe el arreglo y los indices de los elementos a ser intercambiados
%**********************************************************

function array = Swap(array,this,with)
temp = array(this);
array(this) = array(with);
array(with)= temp;
return

---------------------------------------------------------
%*********************************************************
%Implementa la sumatoria desde 'from' hasta 'to'
%Devuelve la suma
%*********************************************************
function [suma] = Sumatoria(from,to, A, x, i,index)

suma = 0;

for k=from:to
suma = suma + A(index(i),k)*x(index(k),1);
end

return

-------------------------------------------------
%**********************************************
%Implements backward substitution over an
%upper triangular matrix
%
%Input: Una matriz A triangular superior y un vector b
%Output: Solution vector
%*********************************************

function [flag,x] = solve_UTSP(A,b,index)
[x_rows,x_cols] = size(b);
[A_rows,A_cols] = size(A);
flag =1;
x = zeros(x_rows,x_cols);

for i=x_rows:-1:1
if abs(A(index(i),i)) < eps
flag = 0;
'Division by zero error'
break;
end
x(index(i),1) = (b(index(i),1)-SumatoriaP(i+1,A_cols,A,x,i,index))/A(index(i),i);
end
return

%************************************************
%Resuelve sistemas lineales utilizando eliminacion
%de Gauss con Pivoteo.
%Input: La matriz A y un vector b
%Output: El vector solucion
%************************************************

function [flag,x] = GE(A,b)

%Para trabajar con la matriz aumentada
A(:,end+1) = b(:,1);

flag =1; %Bandera indicando que no ha ocurrido ningun error
col=1; %Columna pivote
row=1; %inicializar fila pivote
nextRow = row+1;
[b_rowCant,b_colCant] = size(b); %tamano de b
[A_rowCant,A_colCant] = size(A); %tamano de A
index(1:A_rowCant,1) = 1:A_rowCant; %tamano del indice de filas de A

tic; %Comenzar a contar
for col=1:(A_colCant-1)
if abs(A(index(row),col)) < 0.01
[largestPivot,largestPivotIndex] = max(A(index(row):end,col));
largestPivotIndex=largestPivotIndex+(index(row)-1);
index = Swap(index,largestPivotIndex,row);
end

for rowIndex=1:(A_rowCant-row)
if abs(A(index(row),col)) < eps
flag = 0;
'Division by zero error'
return;
end
A(index(nextRow),:) = A(index(row),:)*(-A(index(nextRow),col)/A(index(row),col))+A(index(nextRow),:);
nextRow = nextRow + 1;
end
row = row+1;
col = col+1;
nextRow = row+1;
end

%Para regresar al sistema original
b = A(:,end);
A(:,end) = [];

%Resolver sistema triangular superior
[succesful,x] = solve_UTSP(A,b,index);

x(:,1) = x(index(:),1);
toc;

return

Claro este metodo no lo inventé yo este es solo mi versión del código. Este método fue inventado por Liu Hui y se hace referencia a él en el importantísimo libro chino "Jiuzhang suanshu" o "Los nueve capítulos en el arte de la matemática" que data del 263A.D.. Por otra parte el nombre que conocemos del método, se atribuye al famosísimo matemático alemán, C.F.Gauss, que hizo el mismo popular a finales del siglo 18.

Bueno, pues no se es que es mi opinión que estas cosas tambien son bonitas, tan bonitas como un poema o como una hermosa pintura, pero por alguna razon, no mucha gente se toma el tiempo de sentarse a apreciarlas. Pero, bueno por lo menos quisiera que la mayor cantidad de gente posible se tope al menos con la idea, que sepan al menos que estas cosas existen. Y ya que existen por que no compartirlas? Es como si todos pudieramos tener un Quijote, original, escrita con el mismo puño y letra de Cervantes o una la MonaLisa en la sala de nuestra casa. Quizas no lo parezca pero, si se detienen un segundo a ver lo que hay detras, se daran cuenta de lo inmensamente maravilloso que estas cosas pueden llegar a ser. Remember, "Life moves pretty fast, if you don't stop to look around once in a while, you could miss it".