liegauss.m

来自「1、 写出此两点边值问题的积分形式； 2、 构造有限元方程； 3、 编程计算」· M 代码 · 共 36 行

function [L,U,P]=liegauss(A)
L=[];
U=[];
n=length(A(:,1));
P=eye(n);
for k=1:n-1
    q=0;
    m=0;
    v=[];
    q=A(k,k);
    for i=k:n
        if q<A(i,k)
            q=A(i,k);
            m=i;
        end
    end
    v(1:n)=A(k,1:n);
    A(k,1:n)=A(m,1:n);
    A(m,1:n)=v(1:n);
    v(1:n)=P(k,1:n);
    P(k,1:n)=P(m,1:n);
    P(m,1:n)=v(1:n);
    if A(k,k)~=0
        A(k+1:n,k)=A(k+1:n,k)/A(k,k);
        A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
    else
        '奇异矩阵'
    end
end
L=tril(A);
U=triu(A);
for i=1:n
    L(i,i)=1;
end