前言
在前面的博客中已经更新了基础的线性方程组的求解方法,可以参看这篇文章:
基于MATLAB的求解线性方程组(附完整代码和例题)_唠嗑!的博客-CSDN博客
但是这种基础的方法,在计算大型方程组时,在计算速度、磁盘空间、内存上很不友好。所以本博客提供另外一种思路:先分解+再求解。
一,二,三包括矩阵的三种分解方式,四,五解释如何利用分解的方法求解线性方程组。
一. LU 分解
将A矩阵可以分解为L和U矩阵相乘,如下形式:
观察发现L是一个下三角矩阵,U是一个上三角矩阵。根据L和U矩阵的形式,它们里面的元素与A矩阵元素的关系,可以得到如下方程:
L和U矩阵里面的元素可以形成递推计算公式如下:
其中计算式子中的递推初值如下:
在MATLAB中调用的格式,如下:
[L,U,P]=lu(A)
同样,在此表达式中L是一个下三角矩阵,U是一个上三角矩阵,P代表选主元的置换矩阵。此置换矩阵满足如下等式:
如果在调用时函数时,省去P,如下形式:
[l,u]=lu(A)
那么上面L代表,也就是。
例题1
对以下矩阵A进行LU分解。
解:
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)
对于正交矩阵的相关解释,可以参看这篇文章:
例题2
将以下矩阵A分解成QR矩阵形式。
解:
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分解算法可见如下公式:
在MATLAB中,调用的格式如下:
D=chol(A)
例题3
对矩阵A进行Cholesky分解。
解:
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分解求解方程组
原方程组形式:
现在可以转变为:
所以最终的解的形式:
X=U\(L\b)
之所以这样转变是因为可以大大提高运算速度。
例题4
利用LU分解,计算以下方程组是否有解。
解:
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分解成上三角矩阵和其转置的乘积。
那么,最初方程组如下:
由此可变成:
所以最终的解可以写成:
X=R\(R'\b)
备注:在MATLAB中,转置矩阵表示为R‘
5.2 QR分解
理论上来讲,对于任何长方形矩阵A,都可以进行QR分解,其中Q为正交矩阵,R为上三角矩阵的初等变换形式,也就是如下分解:
那么最初的方程:
可以变形为:
所以,最终的解,可以写成:
X=R\(Q\b)
总结
以上三种分解方式都是为了解决大型方程组。它们的优点包括:运算速度快、节省磁盘空间、节省内存。