《MATLAB智能算法超级学习手册》一一1.4 线性方程组的求解

本节书摘来自异步社区出版社《MATLAB智能算法超级学习手册》一书中的第1章，第1.4节，作者：MATLAB技术联盟 , 高飞 , 许玢更多章节内容可以访问云栖社区“异步社区”公众号查看。

1.4 线性方程组的求解

MATLAB智能算法超级学习手册线性方程组的求解在日常生活中的应用较多，特别是解决企业规划、任务分配等问题。线性方程组的求解一般分为两类：一类是求唯一解或求特解，另一类是求通解。可以通过由MATLAB求解线性方程组系数矩阵的秩来判断：

若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；

若系数矩阵的秩r

线性方程组的通解（无穷解） = 对应齐次方程组的通解 + 非齐次方程组的一个特解，其特解的求法属于解的第一类问题，通解部分属第二类问题。

1.4.1 齐次线性方程组的通解在MATLAB中，函数null( )用来求解零空间，即满足A·X=0的解空间，实际上是求出解空间的一组基。

格式　z = null　　　　 % z的列向量为方程组的正交规范基，满足Z’×Z=I 　　　z =null(A,’r’)　 % z的列向量是方程A·X=0的有理基

解：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　　　 | 　+-　　　　　　　 -+

1.4.2 非齐次线性方程组的通解需要先判断非齐次线性方程组是否有解，若有解，然后求通解，步骤如下。

Step1：判断A·X=b是否有解，若有解，则进行第二步，否则终止求解；

Step2：求A·X=b的一个特解；

Step3：求A·X=0的通解；

Step4：A·X=b的通解等于 A·X=b的通解加上A·X=b的一个特解。

解：在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

解法一：在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

1.4.3 线性方程组的LQ解法函数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=M1×M2 　　symmlq(A,b,tol,maxit,M1,M2,x0) %x0为初始估计值，默认值为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·x0) 　 [x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,…) %resveccg表示每次迭代共轭梯度残差的范数 相关资源：非线性方程组求解matlab程序