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

本文主要是介绍解线性方程组——直接解法:LU分解、PLU分解(类似列主元消去法) | 北太天元,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

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 分解 的优势所在.
在这里插入图片描述

北太天元源代码

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:nfor k = r:1:nA(r,k) = A(r,k) - A(r,[1:r-1])*A([1:r-1],k);endfor m = r+1:1:nA(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/2023y = 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)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

正常情况下,使用 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/2023n = 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:nA(r,k) = A(r,k) - A(r,[1:r-1])*A([1:r-1],k);endfor m = r+1:1:nA(m,r) = (A(m,r) - A(m,[1:r-1])*A([1:r-1],r))*(1/A(r,r));endendL = tril(A,-1)+eye(n);U = triu(A,0);
end

例子

% PA = LU test
%
%   last modified:      09/27/2023
%% 1
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)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,P.*b(i,:))
end%% 2
clc;clear all;
A = [2 1 5 ; 4 4 -4 ;1 3 1];
[L,U,P] = PLU_factorization(A)%% 3 
clc;clear all;
A = [2 3 ; 3 2];
[L,U,P] = PLU_factorization(A)

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

这篇关于解线性方程组——直接解法:LU分解、PLU分解(类似列主元消去法) | 北太天元的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/925982

相关文章

poj 3974 and hdu 3068 最长回文串的O(n)解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大,用原来的O(n^2)会爆。 小白上的O(n^2)解法代码:TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO

[项目][CMP][直接向堆申请页为单位的大块内存]详细讲解

目录 1.系统调用 1.系统调用 Windows和Linux下如何直接向堆申请页为单位的大块内存: VirtualAllocbrk和mmap // 直接去堆上按页申请空间static inline void *SystemAlloc(size_t kpage){#ifdef _WIN32void *ptr = VirtualAlloc(0, kpage << 13,

PageOfficeCtrl支持直接打开服务器磁盘文件

一般来说,PageOfficeCtrl控件的WebOpen方法的第一个参数是待打开文档的URL,此URL可以是相对于当前页面的相对URL,也可以是相对于整个网站根的相对URL,还可以是http开头的完整URL,但是这个URL必须是当前网站的URL,不能跨域。 现在为了更加方便开发者编程,WebOpen支持打开服务器磁盘文件。也就是说,第一个参数可以写成服务器文件的绝对磁盘路径。例如: P

最直接显示 ubuntu 版本号的命令

有时候去看ubuntu版本号,去网上查,很多文章都列出一堆命令,复制命令运行一下,都是打印一些不相关的信息,我只是想看ubuntu版本号而已,能否直接列出版本号就可以了。 有,下面这条命令就是直接的打印出ubuntu版本号, 没有多余信息 lsb_release -a

特征值分解(EVD)和奇异值分解(SVD)—应用于图片压缩

特征值分解(EVD)和奇异值分解(SVD)—应用于图片压缩 目录 前言 一、特征值分解 二、应用特征值分解对图片进行压缩 三、矩阵的奇异值分解 四、应用奇异值分解对图片进行压缩 五、MATLAB仿真代码 前言         学习了特征值分解和奇异值分解相关知识,发现其可以用于图片压缩,但网上没有找到相应代码,本文在学习了之后编写出了图片压缩的代码,发现奇异值分

【UVA】10651-Pebble Solitaire(直接递归或者记忆化)

不知道这个题UVA的数据是怎么的,用2个方法交了,第一次直接递归,第二次记忆化剪枝,时间竟然一样!? 直接郁闷了,简单的二进制表示状态和二进制运算。 14145176 10651 Pebble Solitaire Accepted C++ 0.009 2014-09-04 09:18:21 #include<cstdio>#include<algorithm>#inclu

Anthropic 创始人 Dario Amodei 谈:关于护城河与风险,AI 大很难直接替代人

护城河的迷思   近期,Anthropic创始人Dario Amodei与投资人Erik Torenberg进行了一场引人关注的对话。他们探讨了AI的护城河与潜在风险。话说,护城河就像酒水的保质期,过了时间就得小心别翻车。Amodei提到,AI虽有强大的潜力,但短期内难以完全替代人类的智慧。这可让很多人松了一口气,毕竟机器发热总比人心复杂,听着都觉得不舒服。 聪明与控制的博弈   Dar

2015多校联合训练第一场Assignment(hdu5289)三种解法

题目大意:给出一个数列,问其中存在多少连续子序列,子序列的最大值-最小值< k 这题有三种解法: 1:单调队列,时间复杂度O(n) 2:RMQ+二分,时间复杂度O(nlogn) 3:RMQ+贪心,时间复杂度O(nlogn) 一:RMQ+二分 RMQ维护最大值,最小值,枚举左端点i,二分找出最远的符合的右端点j,答案就是ans += j - i+1;(手推一下就知道) 比如1 2 3

GitHub:代码是程序员沟通最直接的手段

如果不是 Andreessen horowitz 的投资,估计 GitHub 很难被福布斯、CNN、纽约时报等传统媒体注意到。普通大众之前不了解这个工具,是因为它距离记者的世界太远了——GitHub 是一个程序员所使用的托管项目的服务。 但在一些程序员眼里,它不仅是托管项目的地方,还是“开源”项目的大本营,而且是提高程序员“技术水平”和“技术品味”的地方,更是一个程序员社交的地方。

OpenGL/GLUT实践:流体模拟——数值解法求解Navier-Stokes方程模拟二维流体(电子科技大学信软图形与动画Ⅱ实验)

源码见GitHub:A-UESTCer-s-Code 文章目录 1 实现效果2 实现过程2.1 流体模拟实现2.1.1 网格结构2.1.2 数据结构2.1.3 程序结构1) 更新速度场2) 更新密度值 2.1.4 实现效果 2.2 颜色设置2.2.1 颜色绘制2.2.2 颜色交互2.2.3 实现效果 2.3 障碍设置2.3.1 障碍定义2.3.2 障碍边界条件判定2.3.3 障碍实现2.3.