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

也可以采用高斯賽德爾迭代求解。