| 解线性方程组 | 您所在的位置：网站首页 › 95完全分解 › 解线性方程组 |
解线性方程组
|
文章对应视频讲解： PLU分解LU分解L: lower triangular 下三角 U: upper triangular 上三角 LU 分解，顾名思义，为 把一个 矩阵 分成 一个下三角矩阵 乘上一个上三角矩阵的形式。 Example
几个基本的初等行变换，可以自己验算一下，等式的左边与右边是相等的
得到 L 和 U 后，
所以关键是怎么得到 L 和 U 计算顺序
先单独求出 L 和 U 保存为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 来达到一个 列主元消去法的效果, 在计算每一层之前,先把列中最大的那个元素换到相对的第一行, 主要就这一个特点 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文中两次回代所用到的: 解上三角、下三角 |
用上面这几个等式，重新看一下 第一个例子， 
对A进行了三次行变换，得到上三角矩阵U， 两边同时左乘初等矩阵的逆，表示成 A = 啥啥啥 乘 U 再用 Fact4 和 Fact 3 得到 下三角矩阵 L
有了这个形式后，利用矩阵相乘，元素对应相等，便可求出 L 和 U
这样便可得到 x
如果自己来算 就会发现是先算出第一层，才能算出第二层，再算出第三层，等等 因为要用计算机实现，所以需要知道，具体是怎么算的 
在算的过程中可以发现，只在一个矩阵 A 上便可以发生这些变化 也就不需要开 A L U 三个矩阵的存储空间
这样对于 系数矩阵A 相同, 右端常数项 b 不相同的情况下,都可以使用同样的 L,U 进行计算. 所以我把这里写出单独的一步,不然也体现不出 LU 分解 的优势所在. 
【本文地址】