gepp.m

使用matlab进行矩阵的高斯分解,可以对高阶线性方程进行数值求解

function [amax,P,L,U,x]=gepp(A,n,b)
% GEPP   amax is the n-dim vector the maxmum element of every step
% aa is the max element of A;  x is the solution;
% P*L*U=A: 
B=A;        % Ax = b
P=1:n;      % P*A = L*U;P为变换阵的逆
U=zeros(n); % n*n 0矩阵
L=zeros(n); % n*n 0矩阵
amax(1)=max(max(abs(A))); %找出矩阵中绝对值最大的值
for k=1:n-1         
    [s1,s2]=max(abs(B(:,1))); % s1,B中第一列中绝对值最大的值.s2,取得最大值的行数
    B([1,s2],:)=B([s2,1],:);  % 把第一列中绝对值最大处对应的列和第一行交换位置
    m=k+s2-1;                 % 进行交换的列在原来矩阵中对应的列
    P([k m])=P([m k]);        % 把矩阵中(m,k)元素和(k,m)元素交换位置
    U([m k],:)=U([k m],:);    % 把第k行和第m行交换
    L([m k],:)=L([k m],:);    % 把第k行和第m行交换
    piv=B(2:n-k+1,1)/B(1,1);  % 把B中第一列第二到n-k+1元素取出来除以第一个元素赋给piv,共n-k个列元素
    u=B(1,2:n-k+1);           % 把B中第一行第二到n-k+1元素取出来给u，共n-k个行元素
    U(k,k:n)=B(1,:);          % 把B中第一列的元素赋给U中k的后n-k+1个列元素，共n-k+1个列元素
    L(k+1:n,k)=piv;           % 把piv赋给L的第k列的后n-k个列元素
    B(1,:)=[];                % 使B矩阵减少一列
    B(:,1)=[];                % 使B矩阵减少一行
    B= B-piv*u;               % 进行列消去对矩阵中其他元素的改变
    amax(k+1)=max(max(abs(B)));
end
%以上的循环是找到每次变换过程中的最大值一次赋给amax(2),......,amax(n)
U(n,n)=B;                     % 把B中的唯一一个元素赋给U(n,n)
L=L+eye(n);                   % L的对角元素加上1
% Permutation and forward elimination
y = zeros(n,1);               % n维列向量
for k = 1:n
    j = 1:k-1;
    y(k) = b(P(k)) - L(k,j)*y(j);
end

% Back substitution
x = zeros(n,1);
for k = n:-1:1
    j = k+1:n;
    x(k) = (y(k) - U(k,j)*x(j))/U(k,k);
end