LU分解求线性方程组

解一维平板非稳态导热隐式格式时，需要求解线性方程组。LU分解适合线性方程组有唯一解的小规模求解。

function x=LUsolve(A,B,x0,eps,M)
%LU分解法求方程组的解（矩阵公式求解）
%A为方程组的系数矩阵；b为方程组的右端项
%x为线性方程组的解；x0为迭代初值
%eps为误差限；M为迭代的最大次数
if nargin==2
    x0=zeros(size(A,1),1);
    eps= 1.0e-3;%默认精度
    M = 10000;%参数不足时默认后两个条件
elseif nargin==3
    eps= 1.0e-3;%默认精度
    M = 10000;%参数不足时默认后两个条件
elseif nargin ==4
    M = 10000;%参数的默认值
elseif nargin<2
    errordlg('参数不足','错误');
    return
end
[n,m]=size(A);
nb=length(B);
%当方程组行与列的维数不相等时，停止计算，并输出出错信息
if n~=m
    errordlg('矩阵A行数和列数必须相等!','错误');
    return;
end
%当方程组与右端项的维数不匹配时，停止计算，并输出出错信息
if n~=nb
    errordlg('矩阵A的行数必须和b的长度相等!','错误');
    return;
end 
L =zeros(n,n);
U =zeros(n,n);
D =zeros(n,n);
L=-tril(A,-1);
U=-triu(A,1);
D=diag(diag(A));
DLU=(D-L)\U;        %DLU为迭代矩阵
DLB=(D-L)\B;        %DLB为右端项
pr=max(abs(eig(DLU))); %求迭代矩阵谱半径
if pr>=1
    errordlg('迭代矩阵谱半径大于1迭代法不收敛','错误');
    return;
end
k=0;
tol=1;
while tol>=eps 
    x = DLU*x0+DLB;
    k = k+1;         %迭代步数
    tol = norm(x-x0);%前后两步迭代结果的误差
    x0 = x;
    if(k>=M)
        disp('Warning: 迭代次数太多，可能不收敛！');
        return;
    end
end

也可以采用高斯赛德尔迭代求解。