gausselimcomppivot.m

这四个程序分别为高斯消去法、列主元消去法、全主元消去法解线性方程组和Gauss-Jordan消元法求矩阵的逆。 程序采用MATLAB语言开发

%Gaussian Elimination-Complete Pivoting
%Description: This procedure solves a system of linear equations using
%Gaussian Elimination Algorithm with full pivoting strategy. Full
%pivoting is the interchanging of rows and columns in order to place 
%a particularly good element(usually the largest) in the diagonal position.
%Input: the augmented matrix of the system [A b]
%Output: solution of the system and also the debug information if
%necessary.
%Note: This procedure can return the correct solution of the system
%mostly(as far as I can test), it needs to be tested on more input.
%This program is written by Xiaoke Yang @ School of Automation Science and
%Electrical Engineering, Beihang University. You can copy, modify or 
%redistribute it freely. Welcome reports of bugs, suggestions, etc. 
[e-mail]
%Last modified by Xiaoke Yang, Oct.10,2007

% A=[2 -1 1;-3 2 -1;-2 2 1];
% b=[8 -11 -3]';
% A=[1e-12 1;1 2]
% b=[1 3]'
% A=[1 1e12;1 2]
% b=[1e12 3]'
A=[0.012 0.01 0.167;1 0.8334 5.91;3200 1200 4.2]
b=[0.6781 12.1 981]'
N=length(b);
P=1:N;
%forward elimination
for i=1:N-1
    %pivoting
    [Amax,c]=max(max(abs(A(i:N,i:N))));
    [Amax,r]=max(abs(A(i:N,i+c-1)));
    disp(sprintf('The maximum of A= %f, and its position[%d,%d]\n',Amax,i+r-1,i+c-1));
    %swap the row
    T=A(i+r-1,:);
    A(i+r-1,:)=A(i,:);
    A(i,:)=T;
    T=b(i+r-1);
    b(i+r-1)=b(i);
    b(i)=T;
    %end swap the row
    %swap the column
    T=A(:,i+c-1);
    A(:,i+c-1)=A(:,i);
    A(:,i)=T;
    %record the position of the elements in the solution
    P(i)=i+c-1;
    P(i+c-1)=i;
    %end swap the column
    %end pivoting
    disp(sprintf('after %d pivoting,\n [A b]=',i));
    disp([A b]);
    m=-1/A(i,i);
    for j=i+1:N     %j indicating the row
        for k=i+1:N %k indicating the column
            A(j,k)=A(j,k)+A(j,i)*m*A(i,k);
        end
        b(j)=b(j)+A(j,i)*m*b(i);
        A(j,i)=0;
    end
    disp(sprintf('End of %d elimination,\n [A b]=',i));
    disp([A b])
end
%backword substitution
x(N)=b(N)/A(N,N);
for i=N-1:-1:1
    s=0;
    for j=i+1:N;
        s=s+x(j)*A(i,j);
    end
    x(i)=(b(i)-s)/A(i,i);
end
disp(sprintf('The solution to this system is \nx='));
disp(x);
disp(sprintf('The corresponding position of x is \nP='));
disp(P);