大名鼎鼎的LU分解——Matlab解线性方程组(3)

2024-02-06 19:40

本文主要是介绍大名鼎鼎的LU分解——Matlab解线性方程组(3),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

前言

一、伟大的LU分解法

1.前向消去的步骤

2.LU分解有啥用?

二、多出的一步

1.为什么要选主元?

2.为什么会这样捏?

总结


前言

        前文已经完成了铺垫,这里我们开始进行矩阵的LU分解。第一节我们看到了,因为求逆矩阵会导致运算时间增长、运算量增大,所以在一般求解线性方程组时,常用高斯消元的方法。将高斯消元法凝练、规范化,给出矩阵的LU分解法,将系数矩阵分解成L(乘子矩阵)、U(消元后的矩阵)和P(排列矩阵)——自己起的名字,反正就这么个玩意。


一、伟大的LU分解法

        其实之前已经接触到了LU分解的内容,最内核的就是高斯消元法,这也是在大学线性代数课程中教给大家的方法。

        总体上看来,高斯消元法包括了两个阶段:前向消去(forward elimination)和回代(back substitution)。

        前向消去要进行n-1步,其中第k步的工作是:将第k个方程以下的方程分别减去第k个方程的某个倍数,从而消去第k个变量

        回代有n步,意思是当所有前向消去都进行完成后,第n个方程此时应该只剩下一个变量,这时候就可以通过一步除法得出该变量的值,进而反代回第n-1个方程。此时n-1个方程只剩下两个变量,其中一个变量的值由上一步给出,同理循环即可求出所有变量

        这里需要注意一点,在前向消去时,如果第k行第k个变量的系数(主元)很小,则需要将k以下的某个方程中,第k个变量系数最大的那个换上来。这种方法看上去有点多余,但是在计算机运算中是必须的,下一条会讲到。

1.前向消去的步骤

        停!

        这个时候,需要回顾一下高斯消元整体的过程,要包括到每个细节:

        首先找到所有方程第一个变量系数最大的那个方程是哪个,然后把这个变量所在的方程交换到第一个来。这个过程用矩阵语言描述,就是:找到系数矩阵所有行第一列元素中最大的元素所在行,然后将这一行换到第一行来,形成一个新的系数矩阵。        进行这一步,需要进行一次行变换,用矩阵语言表示,就是需要堆系数矩阵A左乘一个排列矩阵P1:

A_{1}=P_{1}A

        然后,确认第一个方程第一个变量的系数为主元,依次消去剩下所有行的第一个变量的系数。具体做法为,先将主元乘以a1倍,然后与第二个方程第一个变量的系数加减抵消成0;再将主元乘以a2倍,与第三行第一个变量的系数加减抵消成0……,直到把第一个方程以下 所有方程 第一个变量 的系数全部 变成0.  用矩阵语言描述为:将第一行第一个元素确定为主元,然后把第一行的某个倍数加给第二行,用以消去第二行第一个元素,以此类推将第一行以下所有行的第一列元素变成0

        这里需要展开讲解:

        用第一节那个例子。假设矩阵为3*3的方阵,第一行主元是10,第二行第三行的第一个元素分别为-3和5 。 想用第一行的10消去第二行的-3,需要给第一行成0.3倍再加到第二行去。   想用第一行的10消去第三行的5,需要给第一行成0.5倍再减到第三行去。用矩阵表示为:

M_{1}=\begin{pmatrix} 1 & 0 &0 \\ 0.3 & 1& 0\\ -0.5 & 0& 1 \end{pmatrix}

(第一行的0.3倍加在第二行,第一行的0.5倍减在第二行)

        这是一个下三角矩阵。可以说,所有进行这种变换的M矩阵都是下三角矩阵。

        那么此时,经过一番变换后,只剩下第一行的第一列有非零元素在,其他行第一列都是0.

        也即:进行了两个矩阵变化后,第一行的第一列(或者说第一个变量)被保留了下来。

A_{1}=M_{1}P_{1}A

        这个式子表示,在消第一列(或者说第一个变量)时,进行了一次排列矩阵左乘变换(行变换)和一次下三角矩阵左乘变换(行变换)。这便是前向消去的第一步。总共有n-1步。

        这n-1步消元全部完成后,会得到一个崭新的系数矩阵U。

U=M_{n-1}P_{n-1}......M_{2}P_{2}M_{1}P_{1}A

值得注意的是,此时产生的U矩阵,由于第1行以下第1列元素全是0,……第k行以下第k列元素全是0……,第n-1行以下第n-1列元素全是0,所以U矩阵是一个妥妥的上三角矩阵

        此时,将所有的M变换矩阵全都挪到方程右边去,经过这个变换后的M矩阵全变成了L矩阵。

L_{1}L_{2}......L_{n-1}U=P_{n-1}......P_{2}P_{1}A

        因为M矩阵和P矩阵都是对A进行行变换的矩阵,所以只要按顺序挪过去,不会影响最终结果。值得注意的是,M矩阵是倍加矩阵(通过标准矩阵进行倍加变换得到的矩阵,名字也是我编的),倍加矩阵移到等式另一边,可以直接给非对角线元素添上负号。所以移到等式左边的L矩阵依然是下三角矩阵

        不妨令:

L=L_{1}L_{2}...L_{n-1}

P=P_{1}P_{2}...P_{n-1}

        方程可以简化为:

LU=PA

        其中:乘子矩阵L包含了所有的乘子,是一个下三角矩阵;U矩阵是消元结束后剩下的系数矩阵,是一个上三角矩阵;P矩阵是初等变换经过排列得到的排列矩阵。

        这个过程,称为矩阵A的LU分解(LU factorization)或者三角分解(triangular decomposition)。

2.LU分解有啥用?

        表面上看,LU分解进行了大批量的矩阵初等变换,得到了一大堆没有什么用的东西。

但是,如果我们把这东西带回去,会发生什么事呢?

对于方程:

Ax=b

由LU分解可知:

LU=PA

A=P^{-1}LU

原方程变为:

P^{-1}LUx=b

LUx=Pb

不妨令Ux=y,则又有:

\left\{\begin{matrix} Ly=Pb\\ Ux=y \end{matrix}\right.

 其中,Pb看作对于已知结果构成的列向量进行行交换变换。通俗点就是把方程结果放一个列向量,起名叫b,对这个b的行与行按照P矩阵规定的那样去交换,得到一个交换过的列向量Pb。

        通过上述式子,已知系数矩阵L和结果向量Pb,可以求出未知向量y。再把未知向量y当成已知向量,结合第二个系数矩阵U求出未知数向量x。

        所以这有什么用呢?

        回想一下上一节的内容——虽然对于一般系数矩阵求解方程有难度,但是对于三角形矩阵当作系数矩阵的线性方程组,有非常简洁快速的算法去求其解。

        而求y时候的系数矩阵L和求x时候的系数矩阵U都是三角形矩阵。

        这样,就把求一般系数矩阵下的方程,化简为——①进行LU分解,②化成两个三角形矩阵的求解。

二、多出的一步

1.为什么要选主元?

        在LU分解中,有一个重要的步骤被我提到——那就是选主元。用人话说,就是在轮到第k行消元的时候,要选择第k行到第n行所有的第k个元素中最大的那个元素作为主元,并把他换到第k行来。

        but why?这一步显得有点多余。正常看来,不论第k行的第k个元素是大是小,经过矩阵基本初等变换后求解出的U矩阵都是唯一的,解出来的x也是唯一的。

        但是不要忘了,我们目前讨论的这些有的没的,都是在思考计算机的算法。包括未来对于大批量数据处理,得到的数据量绝对不是人类可以用纸笔能处理得了的。而用到计算机算法,就要考虑到编程软件自身的精确度。

        考虑上两节讲的那个例子,把这个例子中的系数矩阵简单变一下:

\begin{pmatrix} 10 & -7 &0 \\ -3 &2.099 &6 \\ 5 & -1& 5 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3} \end{pmatrix}=\begin{pmatrix} 7\\ 3.901\\ 6 \end{pmatrix}

易知,这个解仍然是

\begin{pmatrix} 0\\ -1\\ 1 \end{pmatrix}

但是如果我们计算用的电脑计算性能不佳,只能进行5位有效数字的浮点数运算,那么会怎样呢?

\begin{pmatrix} 10 & -7 &0 \\ 0 &-0.001 &6 \\ 0 &2.5& 5 \end{pmatrix}\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3} \end{pmatrix}=\begin{pmatrix} 7\\ 6.001\\ 2.5 \end{pmatrix}

如果不找主元,下一步需要将第二个方程乘以2.5*1000倍后倍加到第三个方程:

(5+(2.5\cdot 10^3)\cdot 6)x_{3}=(2.5+(2.5\cdot 10^3)\cdot 6.001)

因为6.001*2.5*1000=15002.5,而我们的电脑只能保留5位有效数字,所以只能保留下15002.

同理,2.5+15002.5=15005,但是在计算机中,因为有效位数不够,只能做:2+15002=15004这样的运算。

于是,在计算机中,上式变为:(5+15000)x3=(2+15002),即x3=15004/15005=0.99993

而前面说了,x3的理论解为1,看上去差距不是很大。但是当把x3带回去的时候,就会有一些不对劲的东西。

求x2:-0.001x2+6*0.99993=6.001————x2=0.0015/(-0.001)=-1.5

同理,x1=-0.35

最后,计算机返回给我们的解为x1=-0.35;x2=-1.5;x3=0.99993,与理论值有较大差距。

2.为什么会这样捏?

        仔细观察可以发现,在求解x3时,因为有效位数不够,省去了一个0.5。但是实际上两个0.5相加可以进位为1。但是这两个0.5省去以后,x3就产生了一个小误差。

        而这个小误差在反代回x2x1的时候,误差就会被放大,进而对结果产生更大的影响。

        为什么呢?

        回溯过程,我们发现,由于第二个元素的“原主元”比第三行第二个元素小了太多,为了消去下方的元素,必须乘以一个巨大的乘子,而这个巨大的乘子,就会把最低位的有效数字抹去。

        反之如果主元比下方元素大了很多,这时候只需要乘以一个非常小的乘子,而乘一个很小的数并不会“吞掉”有效数字。

        所以应该选择最大的元素作为主元


总结

        本节讲了高斯消元法最关键的一步——LU分解。掌握LU分解之后,可以将求一个一般线性方程组的问题转化成求一次系数矩阵为上三角矩阵的线性方程组和一次系数矩阵为下三角矩阵的线性方程组。而这个代码,上一节已经有了。

这篇关于大名鼎鼎的LU分解——Matlab解线性方程组(3)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



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

相关文章

matlab读取NC文件(含group)

matlab读取NC文件(含group): NC文件数据结构: 代码: % 打开 NetCDF 文件filename = 'your_file.nc'; % 替换为你的文件名% 使用 netcdf.open 函数打开文件ncid = netcdf.open(filename, 'NC_NOWRITE');% 查看文件中的组% 假设我们想读取名为 "group1" 的组groupName

利用matlab bar函数绘制较为复杂的柱状图,并在图中进行适当标注

示例代码和结果如下:小疑问:如何自动选择合适的坐标位置对柱状图的数值大小进行标注?😂 clear; close all;x = 1:3;aa=[28.6321521955954 26.2453660695847 21.69102348512086.93747104431360 6.25442246899816 3.342835958564245.51365061796319 4.87

C# double[] 和Matlab数组MWArray[]转换

C# double[] 转换成MWArray[], 直接赋值就行             MWNumericArray[] ma = new MWNumericArray[4];             double[] dT = new double[] { 0 };             double[] dT1 = new double[] { 0,2 };

libsvm在matlab中的使用方法

原文地址:libsvm在matlab中的使用方法 作者: lwenqu_8lbsk 前段时间,gyp326曾在论坛里问libsvm如何在matlab中使用,我还奇怪,认为libsvm是C的程序,应该不能。没想到今天又有人问道,难道matlab真的能运行libsvm。我到官方网站看了下,原来,真的提供了matlab的使用接口。 接口下载在: http://www.csie.ntu.edu.

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数

Matlab/Simulink中PMSM模型的反电动势系数和转矩系数_matlab pmsm-CSDN博客

MATLAB层次聚类分析法

转自:http://blog.163.com/lxg_1123@126/blog/static/74841406201022774051963/ 层次聚类是基于距离的聚类方法,MATLAB中通过pdist、linkage、dendrogram、cluster等函数来完成。层次聚类的过程可以分这么几步: (1) 确定对象(实际上就是数据集中的每个数据点)之间的相似性,实际上就是定义一个表征

MATLAB的fix(),floor()和ceil()函数的区别与联系

fix(x),floor(x)和ceil(x)函数都是对x取整,只不过取整方向不同而已。 这里的方向是以x轴作为横坐标来看的,向右就是朝着正轴方向,向左就是朝着负轴方向。 fix(x):向0取整(也可以理解为向中间取整) floor(x):向左取整 ceil(x):向右取整 举例: 4个数:a=3.3、b=3.7、c=-3.3、d=-3.7 fix(a)=3 fl

MATLAB中的eig函数

在MATLAB中,计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有5种: E=eig(A):求矩阵A的全部特征值,构成向量E。 [V,D]=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的列向量。 [V,D]=eig(A,'nobalance'):与第2种格式类似,但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量,而格式3直接求矩阵A的特

MATLAB中的diag函数

diag函数功能:矩阵对角元素的提取和创建对角阵 设以下X为方阵,v为向量 1、X = diag(v,k)当v是一个含有n个元素的向量时,返回一个n+abs(k)阶方阵X,向量v在矩阵X中的第k个对角线上,k=0表示主对角线,k>0表示在主对角线上方,k<0表示在主对角线下方。例1: v=[1 2 3]; diag(v, 3) ans =      0     0     0

Matlab simulink建模与仿真 第十章(模型扩展功能库)

参考视频:simulink1.1simulink简介_哔哩哔哩_bilibili 一、模型扩展功能库中的模块概览         注:下面不会对Block Support Table模块进行介绍。 二、基于触发的和基于时间的线性化模块 1、Trigger-Based Linearization基于触发的线性化模块 (1)每次当模块受到触发时,都会调用linmod或者dlinmod函数