下载掌阅APP，畅读海量书库
立即打开
线性方程组的求解在日常生活中的应用较多，特别是解决企业规划、任务分配等问题。线性方程组的求解一般分为两类：一类是求唯一解或求特解，另一类是求通解。可以通过由MATLAB求解线性方程组系数矩阵的秩来判断：
若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；
若系数矩阵的秩r<n，则可能有无穷解；
线性方程组的通解（无穷解）=对应齐次方程组的通解 + 非齐次方程组的一个特解，其特解的求法属于解的第一类问题，通解部分属第二类问题。
在MATLAB中，函数null( )用来求解零空间，即满足A·X=0的解空间，实际上是求出解空间的一组基。
格式 z=null % z的列向量为方程组的正交规范基，满足Z′×Z=I
z=null(A,’r’) % z的列向量是方程A·X=0的有理基
【例1-7】求解方程组
的通解。
解：MATLAB求解程序代码如下。
>>A=[1 2 2 1;
2 1 -2 -2;
1 -1 -4 -3]; %原始系数矩阵
format rat %指定有理式格式
B=null(A,'r') %求解空间的有理基
B=
2 5/3
-2 -4/3
1 0
0 1
或通过最简行得到基：
>> B=rref(A)
B=
1 0 -2 -5/3
0 1 2 4/3
0 0 0 0
则相应地写出线性方程组的通解：
syms k1 k2 %定义符号变量
X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解
%运行结果显示
X=
2*k1 + (5*k2)/3
- 2*k1 - (4*k2)/3
k1
k2
pretty(X) %让通解表达式更精美
+- -+
| 5 k2 |
| 2 k1+ ---- |
| 3 |
| |
| 4 k2 |
| - 2 k1 - ---- |
| 3 |
| |
| k1 |
| |
| k2 |
+- -+
需要先判断非齐次线性方程组是否有解，若有解，然后求通解，步骤如下。
Step1：判断A·X=b是否有解，若有解，则进行第二步，否则终止求解；
Step2：求A·X=b的一个特解；
Step3：求A·X=0的通解；
Step4：A·X=b的通解等于 A·X=b的通解加上A·X=b的一个特解。
【例1-8】求解方程组
解：在MATLAB中建立脚本M文件：
A=[1 -2 3 -1;3 -1 5 -3;2 1 2 -2];
b=[1 2 3]';
B=[A b];
n=4;
RA=rank(A)
RB=rank(B)
format rat
if RA==RB&RA==n %判断是否有唯一解
X=A\b
elseif RA==RB&RA<n %判断是否有无穷解
X=A\b %求特解
C=null(A,'r') %求AX=0的基础解系
else X='equition no solve' %判断无解
end
运行后结果显示为：
RA=
2
RB=
3
X=
equition no solve
【例1-9】求解方程组：
解法一：在MATLAB编辑器中建立M文件：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n
X=A\b
elseif R_A==R_B&R_A<n
X=A\b
C=null(A,'r')
else X='Equation has no solves'
end
运行后结果显示为：
R_A=
2
R_B=
2
Warning: Rank deficient,rank=2,tol=3.826647e-15.
X=
0
0
-8/15
3/5
C=
3/2 -3/4
3/2 7/4
1 0
0 1
所以，原方程组的通解为
解法二：用rref( )求解：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增广矩阵的行最简
运行后结果显示为：
C=
Columns 1 through 5
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
对应齐次方程组的基础解系为：
；非齐次方程组的特解为：
所以，原方程组的通解为X=k 1 ×ξ 1 +k 2 ×ξ 2 +η * 。
函数symmlq的格式如下：
x=symmlq(A,b) %求线性方程组A·X=b的解X。A必须为n阶对称方阵，b为n元列向量。A可以是由afun定义并返回A × X 的函数。如果收敛，将显示结果信息；如果收敛失败，将给出警告信息并显示相对残差norm(b−A·X)/norm(b)和计算终止的迭代次数
symmlq(A,b,tol) %指定误差 tol，默认值是1e−6
symmlq(A,b,tol,maxit) %maxit指定最大迭代次数
symmlq(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子
symmlq(A,b,tol,maxit,M1,M2) %M=M 1 ×M 2
symmlq(A,b,tol,maxit,M1,M2,x0) % x 0 为初始估计值，默认值为0
[x,flag]=symmlq(A,b,…) % flag的取值为：0表示在指定迭代次数内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M 为坏条件的预处理因子；3表示两次连续迭代完全相同；4表示标量参数太小或太大；5表示预处理因子不是对称正定的
[x,flag,relres]=symmlq(A,b,…) % relres表示相对误差norm(b−A·x) /norm(b)
[x,flag,relres,iter]=symmlq(A,b,…) % iter表示计算 x的迭代次数
[[x,flag,relres,iter,resvec]=symmlq(A,b,…) % resvec表示每次迭代的残差：norm(b−A·x 0 )
[x,flag,relres,iter,resvec,resveccg]=symmlq(A,b,…) %resveccg 表示每次迭代共轭梯度残差的范数