基于MATLAB求解大型方程组的方法：LU、QR和乔里斯基算法（附完整代码和例题）

最新推荐文章于 2024-03-25 14:36:51 发布

唠嗑！

最新推荐文章于 2024-03-25 14:36:51 发布

阅读量4.4k

点赞数 6

分类专栏： MATLAB 文章标签： matlab 线性代数矩阵

本文链接：https://blog.csdn.net/forest_LL/article/details/124225194

版权

MATLAB 专栏收录该内容

43 篇文章 90 订阅

订阅专栏

前言

在前面的博客中已经更新了基础的线性方程组的求解方法，可以参看这篇文章：

基于MATLAB的求解线性方程组（附完整代码和例题）_唠嗑！的博客-CSDN博客

但是这种基础的方法，在计算大型方程组时，在计算速度、磁盘空间、内存上很不友好。所以本博客提供另外一种思路：先分解+再求解。

一，二，三包括矩阵的三种分解方式，四，五解释如何利用分解的方法求解线性方程组。

一. LU 分解

将A矩阵可以分解为L和U矩阵相乘，如下形式：

观察发现L是一个下三角矩阵，U是一个上三角矩阵。根据L和U矩阵的形式，它们里面的元素与A矩阵元素的关系，可以得到如下方程：

L和U矩阵里面的元素可以形成递推计算公式如下：

$l_{ij}=\frac{a_{ij}-\sum_{k=1}^{j-1}l_{ik}u_{kj}}{u_{jj}},\quad (j<i)$

$u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ij}u_{kj},\quad (j\geq i)$

其中计算式子中的递推初值如下：

$u_{1i}=a_{1i},\quad i=1,2,\cdots,n$

在MATLAB中调用的格式，如下：

[L,U,P]=lu(A)

同样，在此表达式中L是一个下三角矩阵，U是一个上三角矩阵，P代表选主元的置换矩阵。此置换矩阵满足如下等式：

$PA=LU$

如果在调用时函数时，省去P，如下形式：

[l,u]=lu(A)

那么上面L代表 $P^{-1}L$ ，也就是 $(P^{-1}L)U=A$ 。

例题1

对以下矩阵A进行LU分解。

$A=\begin{pmatrix}1&2&3\\ 2&4&1\\ 4&6&7 \end{}$

解：

MATLAB代码如下：

clc;clear;
A=[1 2 3;2 4 1;4 6 7];
[l,u,p]=lu(A)

运行结果：

l =

1.000000000000000 0 0
0.500000000000000 1.000000000000000 0
0.250000000000000 0.500000000000000 1.000000000000000

u =

4.000000000000000 6.000000000000000 7.000000000000000
0 1.000000000000000 -2.500000000000000
0 0 2.500000000000000

p =

0 0 1
0 1 0
1 0 0

二. QR分解

将矩阵A分解成一个正交矩阵和一个上三角矩阵的乘积。在MATLAB中，调用的格式如下：

[Q,R]=qr(A)

对于正交矩阵的相关解释，可以参看这篇文章：

MATLAB:矩阵基础_唠嗑！的博客-CSDN博客

例题2

将以下矩阵A分解成QR矩阵形式。

$A=\begin{bmatrix}1&2&3\\ 4&5&6\\ 7&8&9\\ 10&11&12\\\end{}$

解：

MATLAB代码如下：

clc;clear;
A=[1 2 3;4 5 6;7 8 9;10 11 12];
[Q,R]=qr(A)

运行结果：
Q =

-0.077615052570633 -0.833052161400749 0.533584619560035 0.123642443234394
-0.310460210282533 -0.451236587425403 -0.803603795405472 0.232853902715677
-0.543305367994433 -0.069421013450063 0.006453732130830 -0.836635135134536
-0.776150525706333 0.312394560525280 0.263565443714604 0.480138789184465

R =

-12.884098726725124 -14.591629883279062 -16.299161039832992
0 -1.041315201750939 -2.082630403501872
0 0 -0.000000000000003
0 0 0

三. 乔里斯基（Cholesky)分解

乔里斯基（Cholesky)分解将矩阵A分解成对角元素为正的三角阵D。当然此种分解的前提要求A为n阶对称正定阵，如下分解形式：

从数学的角度来讲，对称矩阵的Cholesky分解算法可见如下公式：

$d_{ii}=\sqrt{a_{ii}-\sum_{k=1}^{i-1}d_{ik}^2}$

$d_{ij}=\frac{a_{ij-\sum_{k=1}^{j-1}d_{ik}d_{jk}}}{l_{jj}},\quad (j<i)$

在MATLAB中，调用的格式如下：

D=chol(A)

例题3

对矩阵A进行Cholesky分解。

$A=\begin{pmatrix}16&4&8\\ 4&5&-4\\ 8&-4&22 \end{}$

解：
MATLAB代码如下：

clc;clear;
A=[16 4 8;4 5 -4;8 -4 22];
D=chol(A)

运行结果：
D =

4 1 2
0 2 -3
0 0 3

四. 利用LU分解求解方程组

原方程组形式：

$AX=b$

现在可以转变为：

$L*U*X=b$

所以最终的解的形式：

X=U\(L\b）

之所以这样转变是因为可以大大提高运算速度。

例题4

利用LU分解，计算以下方程组是否有解。

$\begin{cases}4x_1+2x_2-x_3=2\\3x_1-x_2+2x_3=10\\11x_1+3x_2=8 \end{}$

解：

MATLAB代码如下：

clc;clear;
A=[4 2 -1;3 -1 2;11 3 0]; %提取系数
B=[2;10;8];
[L,U]=lu(A)
X=U\(L\B)

%验证
b=A*X

运行结果：

L =

0.363636363636364 -0.500000000000000 1.000000000000000
0.272727272727273 1.000000000000000 0
1.000000000000000 0 0

U =

11.000000000000000 3.000000000000000 0
0 -1.818181818181818 2.000000000000000
0 0 0

警告: 矩阵为奇异工作精度。

X =

NaN
Inf
Inf

b =

NaN
NaN
NaN

分析：最终MATLAB显示无解，可能由于A矩阵的系数行列式接近0导致的。

五. Cholesky和QR分解来解方程组

5.1 Cholesky分解

回顾前文，如果A为对称正定矩阵，那么Cholesky分解可以将矩阵A分解成上三角矩阵和其转置的乘积。

那么，最初方程组如下：

$Ax=b$

由此可变成：

$R^TRX=b$

所以最终的解可以写成：

X=R\(R'\b)

备注：在MATLAB中，转置矩阵表示为R‘

5.2 QR分解

理论上来讲，对于任何长方形矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，也就是如下分解：

$A=QR$

那么最初的方程：

$Ax=b$

可以变形为：

$QRX=b$

所以，最终的解，可以写成：

X=R\(Q\b)

总结

以上三种分解方式都是为了解决大型方程组。它们的优点包括：运算速度快、节省磁盘空间、节省内存。

唠嗑！

关注

6
点赞
踩
32

收藏

觉得还不错? 一键收藏
打赏
0
评论
基于MATLAB求解大型方程组的方法：LU、QR和乔里斯基算法（附完整代码和例题）

前言在前面的博客中已经更新了基础的线性方程组的求解方法，可以参看这篇文章：基于MATLAB的求解线性方程组（附完整代码和例题）_唠嗑！的博客-CSDN博客但是这种基础的方法，在计算大型方程组时，在计算速度、磁盘空间、内存上很不友好。所以本博客提供另外一种思路：先分解+再求解。一，二，三包括矩阵的三种分解方式，四，五解释如何利用分解的方法求解线性方程组。一. LU 分解将A矩阵可以分解为L和U矩阵相乘，如下形式：观察发现L是一个下三角矩阵，U是一个上三角矩阵。根据L和U..
复制链接

扫一扫