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

news2024/12/23 18:54:31

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: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)

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/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
%
%   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)

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1611926.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

spring高级篇(二)

1、Aware和InitializingBean Aware和InitializingBean都与Bean的生命周期管理相关。 Aware接口: 概念: Aware接口是Spring框架中的一个标记接口,它表示一个类能够感知到(aware of)Spring容器的存在及其特定的环境。Spring框架提供了多个Awar…

一周学会Django5 Python Web开发-Django5模型数据修改

锋哥原创的Python Web开发 Django5视频教程: 2024版 Django5 Python web开发 视频教程(无废话版) 玩命更新中~_哔哩哔哩_bilibili2024版 Django5 Python web开发 视频教程(无废话版) 玩命更新中~共计47条视频,包括:2024版 Django5 Python we…

小图标还不会设计!

ICON图标设计 hello,我是小索奇 image-20230805225451447 你有好奇过这样的图标如何设计的吗? 其实非常简单,仅需要一行代码即可完成,本篇文章就带伙伴们使用,每天看一篇,简单易懂,日久技长~…

5 款免费数据恢复工具,用于恢复误删数据

数据恢复工具是通过组装幸存的片段、从剩余的片段重建或使用备份来重新获得对因存储损坏、人为错误或意外中断而丢失的文件的访问权限的行为。它是将丢失、损坏、意外擦除或以其他方式无法访问的数据恢复到服务器、计算机、手机或存储设备的过程。 在大多数情况下,数…

kaggle 泰坦尼克使用xgboost 得分0.73684

流程 导入所要使用的包引入kaggle的数据集csv文件查看数据集有无空值填充这些空值提取特征分离训练集和测试集调用模型 导入需要的包 import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import warnings warnings.filterwarni…

人人都能玩赚数字人操作员 数字人直播搭建/多路开播/选品技巧/0-1开播流程

课程目录 01 数字人工业化直播车间打造 02 数字人直播规则及防封技巧 03 数字人直播间搭建步骤流程 04 数字人直播行业应用盘点 05 数字人直播多平台多路开播 06 数字人高成交循环话术运营 07 数字人直播选品及组品技I5 08 数字人直播0-1流程 09 工业化直播0-1流程 网…

Axure中的样式

样式 首先说一下Axure里面的原点位置 如下图: 还有一个办法是我们选中我们的按钮,如上图,然后打开右边的样式,可以看按钮的x,y属性,类似于游戏中unity软件的x,y属性,类似于html中…

【JVM系列】关于静态块、静态属性、构造块、构造方法的执行顺序

💝💝💝欢迎来到我的博客,很高兴能够在这里和您见面!希望您在这里可以感受到一份轻松愉快的氛围,不仅可以获得有趣的内容和知识,也可以畅所欲言、分享您的想法和见解。 推荐:kwan 的首页,持续学…

数字时代的智慧演奏

数字化时代,工业不再是孤独的机器运转,而是演绎着一场智能与数据的华丽交响。无数智能节点的联动,数据的涌动,成为工业的新活力,同时也是创新的源泉。 工业互联网将每个机器、设备连接在一起,打破了原本独立…

【数据结构练习题】堆——top-k问题

♥♥♥♥♥个人主页♥♥♥♥♥ ♥♥♥♥♥数据结构练习题总结专栏♥♥♥♥♥ ♥♥♥♥♥上一章:【数据结构练习题】二叉树(1)——1.相同的树2.另一颗树的子树3.翻转二叉树4.平衡二叉树5.对称二叉树♥♥♥♥♥ 文章目录 1.top-k问题1.1问题描述1.2思路分析1.3绘图分析…

simulink中怎么使用solve函数?

🏆本文收录于「Bug调优」专栏,主要记录项目实战过程中的Bug之前因后果及提供真实有效的解决方案,希望能够助你一臂之力,帮你早日登顶实现财富自由🚀;同时,欢迎大家关注&&收藏&&…

Sulley入门教学——简介、安装(Win7、VMware)

1、简介 Sulley 是由 Pedram Amini 和 Aaron Portnoy 开发的开源工具。它以 Python 编写,可以轻松地在不同平台上部署和使用。Sulley 提供了一个灵活且功能强大的框架,允许用户定义协议消息的结构、字段类型、边界条件和模糊测试策略。用户可以使用 Sul…

第三方应用类---Phpmyadmin 后台 Getshell 操作

免责声明:本节仅做技术交流学习. 目录 什么是Phpmyadmin? getshell前提条件: 详细步骤: 1-搜集到开放phpmyadmin的web,然后访问进去 2-执行SQL命令查看是否开启了读写权限 3-开启了读写权限-->继续 没有开读写权限--->鸡鸡 4-有读写权限之后,执行SQL语句导出文件…

解析 IP(IPv4)地址

IPv 4 地址 一、组成二、IPv4 的分类三、子网掩码四、特殊的地址五、私有 IP 地址六、全局 IP 地址七、私有 IP 地址和全局 IP 地址的关系八、广播地址九、网络地址十、IP 地址个数计算十一、查看电脑的 IP 地址(window)十二、手动设置电脑的 IP 地址 为…

第二篇:阿里为什么要禁用Java内置线程池?自定义线程池参数如何设置?

线程池基本概念不清楚的请移步:第一篇:Java中的线程池介绍(阿里为什么要禁用Java内置线程池?自定义线程池参数如何设置?)-CSDN博客 Java中为我们提供了内置线程池,为什么阿里却要严令禁止使用呢…

Python 数据结构和算法实用指南(一)

原文:zh.annas-archive.org/md5/66ae3d5970b9b38c5ad770b42fec806d 译者:飞龙 协议:CC BY-NC-SA 4.0 前言 数据结构和算法是信息技术和计算机科学工程学习中最重要的核心学科之一。本书旨在提供数据结构和算法的深入知识,以及编程…

kkFileView 任意文件上传致远程代码执行漏洞复现(QVD-2024-14703)

0x01 产品简介 kkFileView 是使用 spring boot 搭建的文件文档在线预览解决方案,支持主流办公文档的在线预览。 0x02 漏洞概述 2024年4月,互联网上披露kkFileView远程代码执行漏洞情报,攻击者可利用该漏洞上传恶意文件,获取操作系统权限。该漏洞利用简单,建议受影响的客…

【Linux开发 第六篇】Linux常用命令知识

常用命令知识 关机和重启用户管理用户组运行级别帮助指令文件目录类时间日期类搜索查找类压缩和解压类 关机和重启 shutdown -h now //立刻进行关机 shutdown -h 1 “1分钟后进行” //一分钟后进行关机 并向登录Linux的各个用户发送字符串 shutdown -r now //现在重新启动…

网络安全产品---态势感知EDR

态势感知 what SA,Situational Awareness 是对一定时间和空间内的环境元素进行感知,并对这些元素的含义进行理解,最终预测这些元素在未来的发展状态。 why 安全防护思想已经从过去的被动防御向主动防护和智能防护转变。如果不做到主动防御…

git工具简单使用

文章目录 git上传克隆README.gitignore常用指令冲突 git 进行版本控制的版本控制器。安装git yum install -y git 配置git git config --global user.email "youexample.com" 告诉git你的邮箱是什么?最好输入你的gitee的注册邮箱git config --global …