基于MATLAB的稀疏矩阵与符号求解线性方程组（附完整代码）

唠嗑！

已于 2024-01-08 10:15:50 修改

阅读量8.5k

点赞数 12

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

于 2022-04-19 17:23:47 首次发布

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

版权

MATLAB 专栏收录该内容

45 篇文章 103 订阅

订阅专栏

前言

前面文章已经更新了至少三种求解线性方程组的MATLAB指令，分别是：

逆矩阵法：x=inv(A)*b，或者写成 x=A^(-1)*b
伪逆法：x=pinv(A)*b
左除法：x=A\b

本节将更新另外两种方法，符号矩阵法与稀疏矩阵求解法。

一. 符号解法

在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，MATLAB格式如下：

linsolve(sym(A),sym(b)) %结果默认小数表示

还有另外一种函数也可以，形式如下：

X=sym(A)\sym(b) %默认分数形式结果显示

又或者是直接写成多个方程的形式，如下：

solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')

例题1

利用两种不同的方法求解如下方程组：

$\begin{pmatrix}10&-1&0\\ -1&10&-2\\ 0&-2&10\\ \end{pmatrix}\begin{pmatrix}x_1\\x_2\\x_3 \end{pmatrix}=\begin{pmatrix}9\\7\\6 \end{pmatrix}$

解：

MATLAB代码如下：

clc;clear;
A=[10 -1 0;-1 10 -2;0 -2 10];
b=[9;7;6];

ans1=linsolve(A,b) %方法一

ans2=sym(A)\sym(b)  %方法二

运行结果：
ans1 =

0.9958
0.9579
0.7916

ans2 =

473/475
91/95
376/475

二. 稀疏矩阵的生成

2.1 基本方法生成稀疏矩阵

稀疏矩阵，在数值分析中，是其元素大部分为零的矩阵。反之，如果大部分元素都非零，则这个矩阵是稠密的。在科学与工程领域中求解线性模型时经常出现大型的稀疏矩阵。

稀疏矩阵在MATLAB建立的格式，如下：

S=sparse(i,j,s,m,n)

此代码可以生成 $m\times n$ 阶的稀疏矩阵，以向量i和j为坐标的位置上对应的元素值为s。

例题2

利用MATLAB生成任意 $5\times 5$ 的稀疏矩阵。

解：

MATLAB代码如下：

clc;clear;
n=5;
A=sparse(1:n,1:n,4*ones(1,n),n,n)%生成的结果中，第一列代表位置，第二列代表元素

运行结果：

A =

(1,1) 4
(2,2) 4
(3,3) 4
(4,4) 4
(5,5) 4

例题3

利用MATLAB生成任意 $5\times 5$ 稀疏矩阵，要求里面非零元素仅有四个，且利用full()函数将结果表示成完整的矩阵。

解：

MATLAB代码如下：

clc;clear;
n=5;
A=sparse(2:n,1:n-1,ones(1,n-1),n,n);
full(A)

运行结果：

ans =

0 0 0 0 0
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0

例题4

n=5，利用稀疏矩阵建立主对角线上元素为4，两条次对角线为1的三对角矩阵。

解：

MATLAB代码如下：

clc;clear;
n=5;
A1=sparse(1:n,1:n,4*ones(1,n),n,n); %主对角线
A2=sparse(2:n,1:n-1,ones(1,n-1),n,n); %上对角线
A=full(A1+A2+A2') %分别代表三条对角线

运行结果：

A =

4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4

2.2 指定列生成稀疏矩阵

生成 $m\times n$ 阶稀疏矩阵，且把B的列放在由d指定的位置。MATLAB代码如下：

A=spdiags(B,d,m,n)

例题5

利用spdiags()函数生成例题4中的矩阵。

解：

MATLAB代码如下：

clc;clear;
n=5;
A=spdiags([ones(n,1),4*ones(n,1),ones(n,1)],[-1,0,1],n,n);
%此例子B有三列，0代表正对角线
full(A)

运行结果：
ans =

4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4

2.3 无规律的稀疏矩阵

对于无规律的稀疏矩阵，可以使用以下MATLAB指令将外部数据转化为稀疏矩阵：

A=spconvert(dd)

备注：上述格式中的dd为文件名，可根据实际情况进行修改

调用的形式：先用load函数加载以行表示对应位置和元素值的.dat文本文件，再利用上述命令将其转化为稀疏矩阵。

三. 稀疏矩阵的计算

同满矩阵相比，稀疏矩阵在算法上有很大的不同，其中就包括存储空间减少，计算时间减少。

例题6

利用稀疏矩阵和满矩阵两种方法来比较求解下面方程组时的区别。（n=1000)

解：

MATLAB代码如下：

clc;clear;
n=1000;
A1=sparse(1:n,1:n,4*ones(1,n),n,n);
A2=sparse(2:n,1:n-1,ones(1,n-1),n,n);
A=A1+A2+A2';
b=ones(1000,1);

%稀疏矩阵求解
tic;x=A\b;
t1=toc  %运行时间

%满矩阵求解
a=full(A);
tic;x=a\b;
t2=toc

运行结果：

t1 =

0.0012

t2 =

0.0180

备注：不同的MATLAB版本，每次运行的时间可能会略微有区别

分析：稀疏矩阵求解运行的时间会更短

四. 四种稀疏矩阵求解方程法

4.1 最小残差法

最小残差法可以求取线性方程组Ax=b的最小范数解x。MATLAB格式如下：

x=minres(A,b,tol,maxit)

此函数要求 $n\times n$ 的稀疏矩阵A必须为对称的，但无需是正定的。列向量b的长度要求为n。此函数中的“tol"代表精度，”maxit"代表预计最大迭代次数。

例题7

将例题6中的稀疏矩阵A的1全改为-2，n的值改为100，方程组右侧b的第一项为2，最后一项为2，其余均为零。剩下的条件与例题6保持不变，利用minres函数求解此方程组。

解：

MATLAB代码如下：

clc;clear;

n=100;
on=ones(n,1);
A=spdiags([-2*on 4*on -2*on],-1:1,n,n);
%主对角全为4，上下对角线全为-2
b=sum(A,2);
tol=1e-10;
maxit=100;
x=minres(A,b,tol,maxit);
%结果较长就不显示了，看看最后的实际迭代次数与相对误差

运行结果：

minres 在解的迭代 50 处收敛，并且相对残差为 4.7e-14。

4.2对称的LQ方法

对称的LQ方法尝试针对x对线性方程组Ax=b求解。前提要求和最小残差法要求一样， $n\times n$ 的稀疏矩阵A必须为对称的，但无需是正定的。列向量b的长度要求为n。此函数中的“tol"代表精度，”maxit"代表最大迭代次数。MATLAB调用格式，如下：

x=symmlq(A,b,tol,maxit)

例题8

此题条件和例题7完全一样，将调用求解函数改为symmlq()，求解此函数的实际迭代次数和相对误差。

解：
MATLAB代码如下：

clc;clear;

n=100;
on=ones(n,1);
A=spdiags([-2*on 4*on -2*on],-1:1,n,n);
%主对角全为4，上下对角线全为-2
b=sum(A,2);
tol=1e-10;
maxit=100;
x=symmlq(A,b,tol,maxit);
%结果较长就不显示了，看看最后的实际迭代次数与相对误差

运行结果：

symmlq 在解的迭代 49 处收敛，并且相对残差为 4.3e-15。

4.3 双共轭梯度法

第三个可调用的函数形式如下：

x=bicg(A,b,tol,maxit)

例题9

此题条件将例题7中上三角的-2改为-1，其他不变，利用函数bicg()，求解此函数的实际迭代次数和相对误差。

解：

MATLAB代码如下：

clc;clear;

n=100;
on=ones(n,1);
A=spdiags([-2*on 4*on -on],-1:1,n,n);
%主对角全为4，上下对角线全为-2
b=sum(A,2);
tol=1e-8;
maxit=100;
x=bicg(A,b,tol,maxit);
%结果较长就不显示了，看看最后的实际迭代次数与相对误差

运行结果：

bicg 在解的迭代 36 处收敛，并且相对残差为 5.9e-09。

4.4 LSQR方法

LSQR方法尝试为x对线性方程Ax=b求解，否则将尝试计算使得norm(b-A*x)最小的最小二乘解x。其中 $m\times n$ 系数矩阵A无需为方阵，要求为大型稀疏矩阵。列向量b的长度要求为m。MATLAB的调用形式如下：

x=lsqr(A,b,tol,maxit)

例题10

此题条件和例题9完全一样，将调用求解函数改为lsqr()，求解此函数的实际迭代次数和相对误差。

解：

MATLAB代码如下：

clc;clear;

n=100;
on=ones(n,1);
A=spdiags([-2*on 4*on -on],-1:1,n,n);
%主对角全为4，上下对角线全为-2
b=sum(A,2);
tol=1e-8;
maxit=100;
x=lsqr(A,b,tol,maxit);
%结果较长就不显示了，看看最后的实际迭代次数与相对误差

运行结果：

lsqr 在解的迭代 62 处收敛，并且相对残差为 9.5e-09。

总结

到目前为止，有关线性方程组的解法已经更新完了，现在做一点点总结。

通常情况下，利用逆矩阵法和左除法得到的均是数值解。若要得到解析解，需要将数值矩阵转化为符号矩阵。
左除法比逆矩阵法具有更好的数值稳定性，更快的运算速度。而且左除法还适合A不是方阵的情况
当矩阵A接近奇异矩阵时，逆矩阵法和左除法可能会得出错误的结果。这时需要参考其他方法

唠嗑！

关注

12
点赞
踩
81

收藏

觉得还不错? 一键收藏
打赏
1
评论
基于MATLAB的稀疏矩阵与符号求解线性方程组（附完整代码）

前言前面文章已经更新了至少三种求解线性方程组的MATLAB指令，分别是：逆矩阵法：x=inv(A)*b，或者写成 x=A^(-1)*b 伪逆法：x=pinv(A)*b 左除法：x=A\b本节将更新另外两种方法，符号矩阵法与稀疏矩阵求解法。一. 符号解法在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，MATLAB格式如下：linsolve(sym(A),sym(b)) %结果默认小数表示还有另外一种函数也可以，形式如下：X=s..
复制链接

扫一扫