Matlab广义追赶法（Thomas法）

矩阵分解方法（广义追赶法）

对于一般形式的线性方程组

$Ax=b$

将系数矩阵A分解成下三角阵L和上三角阵U

$A=L{U}_1$

化归为两个三角方程组

$Ly=b$，$Ux=y$

对于3阶矩阵：

$$\begin{gathered} A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{11}& a_{12}& a_{13}\\ a_{11}& a_{12}& a_{13}\end{array}\right] \\ =\left[\begin{array}{ccc}{l}_{11}& 0& 0\\ l_{21}& l_{22}& 0\\ l_{31}& l_{32}& l_{33}\end{array}\right]\left[\begin{array}{ccc}1& u_{12}& u_{13}\\ 0& 1& u_{23}\\ 0& 0& 1\end{array}\right]=L{U}_1 \end{gathered}$$

展开

$\begin{array}{rl}{l}_{11}=a_{11},& u_{12}=a_{12}/l_{11},u_{13}=a_{13}/l_{11}\\ l_{21}=a_{21},& l_{22}=a_{22}-l_{21}{u}_{12},u_{23}=(a_{23}-l_{21}{u}_{13})/l_{22}\\ l_{31}=a_{31},& l_{32}=a_{32}-(l_{31}{u}_{12}),l_{33}=a_{33}-l_{31}{u}_{13}-l_{32}{u}_{23}\end{array}$

求解环节

1.预处理

实现矩阵分解$A=L{U}_1$:对$i=1,2,...,n$计算

$l_{i,j}=a_{i,j}-{\sum }_{k=1}^{j-1}{l}_{i,k}{u}_{k,j},j=1,2,...,i$

$u_{i,j}=(a_{i,j}-{\sum }_{k=1}^{j-2}{l}_{i,k}{u}^{k,j})/l_{i,i},j=i+1,i+2,...,n$

2.追的过程

解下三角方程组$Ly=b$即

${\sum }_{j=1}^i{l}_{i,j}{y}_j=b_i,i=1,2,...,n$

回代公式

$y_i=（b_i-{\sum }_{j=1}^{i-1}{l}_{i,j}{y}_j）/l_{i,i},i=1,2,...,n$

3.赶的过程

解单位上三角方程$U_{1}x=y$即

$x_i+{\sum }_{j=i+1}^n{u}_{i,j}{x}_j=y_j,i=1,2,...,n$

回代公式为

$x_i=y_i-{\sum }_{j=i+1}^n{u}_{i,j}{x}_j,i=n,n-1,...,1$

Matlab实现：

%% 追赶法（Thomas法）
A=[1  1  1;
   2 -1 -2;
   3  1 -4];
b=[17;1;3];
m=length(A);%矩阵大小
L=zeros(m);%下三角方程
U=zeros(m)+eye(m);%上三角方程
y=zeros(m,1);%下三角方程解初始化
x=zeros(m,1);%解向量初始化
%% 1.预处理
for i=1:m
    for j=1:i
        %计算矩阵L的各个元素，其中列坐标j=1,2,..,i
        K=0;
        if j<=1
            L(i,j)=A(i,j)-K;
        else
            for k=1:j-1
                K=L(i,k)*U(k,j)+K;
            end
            L(i,j)=A(i,j)-K;
        end
        %计算矩阵U的各元素，其中列坐标jj=i+1,i+2,...,n
        KK=0;
        if i<m;
            for jj=i+1:m
                if i==1;
                    U(i,jj)=(A(i,jj)-KK)./L(i,i);
                else
                    for k=1:i-2
                        KK=L(i,k)*U(k,i+1)+KK;
                    end
                    U(i,jj)=(A(i,i+1)-KK)./L(i,i);
                end
            end
        else
            
        end
        
    end
end
%% 追的过程 解Ly=b
for i=1:m
    u=0;
    j=1;
    if j>i-1
        y(i,1)=(b(i,1)-u)./L(i,i);
    else
        for j=1:i
            u=L(i,j)*y(j,1)+u;
        end
        y(i,1)=(b(i,1)-u)./L(i,i);
    end
end
%% 赶的过程 解Ux=y
for k=1:m
    i=m-k+1;
    jj=i+1;
    u=0;
    if jj>m
        x(i,1)=y(i,1)-u;
    else
        for j=1:m
            u=U(i,j)*x(j,1)+u;
        end
        x(i,1)=y(i,1)-u;
    end
end
%% 显示结果
disp('解向量为')
x
disp('下三角阵为')
L
disp('上三角阵为')
U

结果：
解向量为

x =

4.4706
7.9412
4.5882

下三角阵为

L =

1.0000         0         0
2.0000   -3.0000         0
3.0000   -2.0000   -5.6667

上三角阵为

U =

1.0000    1.0000    1.0000
     0    1.0000    0.6667
     0         0    1.0000

可以写成函数形式:function ($x$,$L$,$U_1$)=Thomas_Mathod($A$,$b$)

输入为$A$,$b$

输出$x$,$L$,$U_1$