Method of Principal Elements

by

 

This script implements the method presented by B.P. Demidivich, I.A. Maron 1981

Demidopivot.m
function Demidopivot(a,b)
%%     THE METHOF of PRINCIPAL ELEMENTS
%
%    See, Computational Mathematics
%         B.P. Demidivich, I.A. Maron
%         Third Printing, 1981, MIR.
%         page 278.
%
% Implemented by Prof. Diomar Cesar Lobo PhD
% UFF - Volta Redonda, RJ, Brazil
% April 16th, 2013
%OBS: Matrix A must be Augmented matriz [a b]
%%
A = [a b]; atemp=A;
disp(A)
[n,m]=size(A);
%
%Initialize vector to keep track of column swaps.
for i = 1:n
    col(i) = i;
end
%
%Do Demidovich Method of Principal Elements, pg. 278. 
for k = 1:n-1
    [n,m]=size(A);
    max = 0; 
    isav = k;
    jsav = k;
    for i = k:n  % Find pivot a(p,q)
        for j=k:n
            mtest = A(i,j);
            if (abs(mtest)) > max 
              max = abs(A(i,j));
              isav = i;
              jsav = j;
            end
        end
    end
    out=sprintf('valor amax=%f, imax=%d, jmax=%d\n',max,isav,jsav);disp(out);
    %
    % calculate the multipliers mi=-a(i,q)/a(p,q) 
	for i =k:n
        if (i ~= isav)
		    factorm(i)=-A(i,jsav)/A(isav,jsav);
        else
            if (i == isav)
                factorm(i)=-A(isav,jsav)/A(isav,jsav);
            end
        end
    end
% calculate the M matrix 
for i = k:n
    if(i ~= isav)
		for j = k:m
			atemp(i,j)=A(i,j)+A(isav,j)*factorm(i);
        end
    end
end
%
atemp
A = atemp;
% rebuild matrix in upper triangular form
    if (isav ~= k)
        %switch rows.
        for j = 1:m
            temp = A(k,j);
            temp1 = A(isav,j);
            A(k,j) = temp1;
            A(isav,j) = temp;
        end 
            disp('swap rows');
            disp(A);
    end
    if (jsav ~= k)
        %switch columns.
        tempcol = col(jsav);
        col(jsav) = col(k);
        col(k) = tempcol;
        for i=1:n
            temp = A(i,k);
            temp1 = A(i,jsav);
            A(i,jsav) = temp;
            A(i,k) = temp1;
        end
            disp('swap columns');
            disp(A);
            disp('variables now in order: ');
            disp(col');       
    end
    atemp = A;
    %
end   %End of matrix M sequence
%
if (A(n,n) == 0) %Check for singular matrix.
    disp('matrix of coefficients is singular')
    disp(A);
    return;
end
%
%Do back substitution.
for k = 1:m-n
    x(k,n) = A(n,n+k)/A(n,n);
end
for k = 1:m-n
    for i = n-1:-1:1
        sum(k) = 0;
        for j = i+1:n
            sum(k) = sum(k) + A(i,j)*x(k,j);
        end
        x(k,i) = A(i,n+k) - sum(k);
        if (A(i,i) == 0) %Check for singular matrix.
            disp('matrix of coefficients is singular');
            disp(A);
            break;
        end
        x(k,i) = x(k,i)/A(i,i);
    end
end
%
%Display upper triangular form.
disp('upper triangular form:');
disp(A);
%
%Swap solution vector if columns have been exchanged.
%
    for i = 1:n
        if(col(i) ~= i)
          disp('current order of variables: ');
          disp(col');
          pause;
          disp('switching back to original order ');
          disp(' ');
          break;
        end
    end
%
%Display solution.
disp('solution:');
%
for j = 1:m-n 
    for i = 1:n
        xo(j,col(i)) = x(j,i);
    end
end
disp(xo');
return;
%

Contact us