解线性方程组——直接解法：LU分解、PLU分解(类似列主元消去法) | 北太天元 or matlab

文章对应视频讲解：

• PLU分解
• LU分解

  ---

  L: lower triangular 下三角

  U: upper triangular 上三角

  LU 分解，顾名思义，为 把一个 矩阵 分成 一个下三角矩阵 乘上一个上三角矩阵的形式。

  Example

  

  为什么可以这样

  几个基本的初等行变换，可以自己验算一下，等式的左边与右边是相等的

  

  用上面这几个等式，重新看一下 第一个例子，

  

  对A进行了三次行变换，得到上三角矩阵U，

  两边同时左乘初等矩阵的逆，表示成 A = 啥啥啥 乘 U

  再用 Fact4 和 Fact 3 得到 下三角矩阵 L

  LU分解

  

  有了这个形式后，利用矩阵相乘，元素对应相等，便可求出 L 和 U

  得到 L 和 U 后，

  

  这样便可得到 x

  所以关键是怎么得到 L 和 U

  计算顺序

  

  如果自己来算

  就会发现是先算出第一层，才能算出第二层，再算出第三层，等等

  因为要用计算机实现，所以需要知道，具体是怎么算的

  

  

  在算的过程中可以发现，只在一个矩阵 A 上便可以发生这些变化

  也就不需要开 A L U 三个矩阵的存储空间

  LU分解算法

  先单独求出 L 和 U

  

  这样对于 系数矩阵A 相同, 右端常数项 b 不相同的情况下,都可以使用同样的 L,U 进行计算.

  所以我把这里写出单独的一步,不然也体现不出 LU 分解 的优势所在.

  

  北太天元 or Matlab 实现

  LU分解

  function  [L,U] = LU_factorization(A)
  % LU分解
  % A ： 系数矩阵
  % A = LU
  %   Version:            1.0
  %   last modified:      09/25/2023
  n = length(A);
  A([2:n],1) = A([2:n],1) * (1/A(1,1)); 
  for r = 2:1:n
      for k = r:1:n
          A(r,k) = A(r,k) - A(r,[1:r-1])*A([1:r-1],k);
      end
      for m = r+1:1:n
          A(m,r) = (A(m,r) - A(m,[1:r-1])*A([1:r-1],r))*(1/A(r,r));
      end
  end
  L = tril(A,-1)+eye(n);
  U = triu(A,0);
  end
  

  保存为LU_factorization.m文件

  两次回代

  function [X] = back_substitution_two(L,U,b)
  % Ly=b, Ux=y
  % b : 列向量
  % X : 解向量
  %
  %   Version:            1.0
  %   last modified:      09/25/2023
      y = push_ltm(L,b);
      X = reg_utm(U,y);
  end
  

  保存为back_substitution_two.m文件

  简单使用一下

  clc,clear all;
  A = [1 2 -1;2 1 -2; -3 1 1];
  b1 = [3 3 -6];
  [L,U] = LU_factorization(A);
  X1 = back_substitution_two(L,U,b1)
  

  得

  L = 
     1.000000000000000   0.000000000000000   0.000000000000000
     2.000000000000000   1.000000000000000   0.000000000000000
    -3.000000000000000  -2.333333333333333   1.000000000000000
  U = 
     1   2  -1
     0  -3   0
     0   0  -2
  X1 = 
     3.000000000000000
     1.000000000000000
     2.000000000000000
  

  对于b不同的情况，可以这样

  clc,clear all;
  A = [1 2 -1;2 1 -2; -3 1 1];
  b1 = [3 3 -6];
  [L,U] = LU_factorization(A)
  b = [3 3 -6;1 2 5;4 9 8;10 2 5];
  m = length(b); X = cell(1,m);
  for i = 1:1:length(b)
      X{i} = back_substitution_two(L,U,b(i,:)')
  end
  

  有3个不同的 b 得 3个不同的 解向量，这里是元胞数组的表示

  X = 
      {3x1 double}    {3x1 double}    {3x1 double}    {3x1 double}
  

  正常情况下,使用 Gauss消去法的话, Ax=b下,

  相同的A 不同的 b,我们对于每一个b 都需要进行一套完整的消元过程,最后再进行一次回代.

  计算量相当于: k 次完整消元+ k次回代

  如果使用 LU分解, 则只需要进行一次完整的消元过程,加 2k 次回代

  计算量相当于: 一次完整消元 + 2k 次 回代

  显然 LU分解使用起来会更方便一些.

  当然,上面的LU分解还没有达到列主元消去法那样的精度,只是相当于基础版的Gauss消去法

  下面来简单介绍一下 PLU 分解,相当于 列主元消去法

  PLU分解

  主要是 通过 P 来达到一个 列主元消去法的效果,

  在计算每一层之前,先把列中最大的那个元素换到相对的第一行, 主要就这一个特点

  

  

  北太天元 or Matlab 实现

  PLU分解

  function  [L,U,P] = PLU_factorization(A)
  % PA = LU分解
  % Input: A
  % output: L,U,P
  %   Version:            1.0
  %   last modified:      09/27/2023
      n = length(A);
      % 第一次行交换
      [~,s]= max(A(1:n,1)); % s 表示第一列最大元素的位置
      P = eye(n);
      P([1,s],:) = P([s,1],:); 
          A = P*A; % 用初等矩阵左乘A 对 A 作行交换
      A([2:n],1) = A([2:n],1) * (1/A(1,1)); % 求第一层
      for r = 2:1:n
          % 先有 行交换
          p=eye(n);  % 用 p 记录每一次的初等矩阵
          [~,s]= max(A(r:n,r));
           s =  s + r-1;
          p([r,s],:) = p([s,r],:);  
          A = p*A; % A的改变
          P=p*P; % 记录P的变化
              % 求第 r 层
          for k = r:1:n
              A(r,k) = A(r,k) - A(r,[1:r-1])*A([1:r-1],k);
          end
          for m = r+1:1:n
              A(m,r) = (A(m,r) - A(m,[1:r-1])*A([1:r-1],r))*(1/A(r,r));
          end
      end
      L = tril(A,-1)+eye(n);
      U = triu(A,0);
  end
  

  例子

  % PA = LU test
  clc;clear all;
  A = [1 2 -1;2 1 -2; -3 1 1];
  b1 = [3 3 -6]';
  [L,U,P] = PLU_factorization(A)
  X1 = back_substitution_two(L,U,P*b1)
  

  运行后得

  L = 
     1.000000000000000   0.000000000000000   0.000000000000000
     0.500000000000000   1.000000000000000   0.000000000000000
    -1.500000000000000   1.666666666666667   1.000000000000000
  U = 
     2.000000000000000   1.000000000000000  -2.000000000000000
     0.000000000000000   1.500000000000000   0.000000000000000
     0.000000000000000   0.000000000000000  -2.000000000000000
  P = 
     0   1   0
     1   0   0
     0   0   1
  X1 = 
     3
     1
     2
  

  文中两次回代所用到的: 解上三角、下三角