1、第6章 线性方程组求解,线性方程的求解分为两类：,（1）方程组求唯一解或求特解， （2）方程组求无穷解即通解。秩=n，有唯一解；秩n，则可能有无穷解；非齐次线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解； 其特解的求法属于解的第一类问题， 通解部分属第二类问题。,一、求线性方程组的唯一解或特解（第一类问题）,这类问题的求法分为两类： （1）解低阶稠密矩阵 直接法； （2）解大型稀疏矩阵 迭代法。,1利用矩阵除法求线性方程组的特解（或一个解）,方程：AX=b 解法：（1）X=Ab（2）rref(A,b) % 求行最简形,例1 求方程组,的解,解：,A=5 6 0 0 01
2、5 6 0 00 1 5 6 0 0 0 1 5 60 0 0 1 5; B=1 0 0 0 1; R_A=rank(A) %求秩 X=AB %求解,运行后结果如下 R_A =5 X =2.2662-1.72181.0571-0.59400.3188,用函数rref求解：, C=A,B %增广矩阵C R=rref(C) %将C化成行最简形 R =1.0000 0 0 0 0 2.26620 1.0000 0 0 0 -1.72180 0 1.0000 0 0 1.05710 0 0 1.0000 0 -0.59400 0 0 0 1.0000 0.3188则R的最后一列元素就是所求之解。,例2
3、 求方程组,的一个特解。,A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; B=1 4 0; X=AB %由于系数矩阵不满秩，该解法可能存在误差。 X = 0 0 -0.5333 0.6000,若用rref求解，则比较精确： A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; B=1 4 0; C=A,B; %构成增广矩阵 R=rref(C) R =1.0000 0 -1.5000 0.7500 1.25000 1.0000 -1.5000 -1.7500 -0.25000 0 0 0 0 由此得解向量X=1.2500 0.2500 0 0（一个特解）。,2利用矩
4、阵的LU、QR和cholesky分解求方程组的解,（1）LU分解：LU分解又称Gauss消去分解，可把任意 方阵分解为下三角矩阵的基本变换形式（行 交换）和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则：A*X=b 变成L*U*X=b 所以X=U(Lb) 这样可以大大提高运算 速度。命令 L，U=lu (A),例3 求方程组的一个特解 A=4 2 -1;3 -1 2;11 3 0; B=2 10 8; D=det(A) L,U=lu(A) %A=LU X=U(LB),显示结果如下：,D =0 L =0.3636 -0.5000 1.00000.2727 1.0000 01.000
5、0 0 0 U =11.0000 3.0000 00 -1.8182 2.00000 0 0.0000,Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 2.018587e-017. X =1.0e+016 *-0.40531.48621.3511 A*X,ans =088,B =2108,（2）Cholesky分解,若A为对称正定矩阵，则有A=RR 其中R为上三角阵。 方程 A*X=b 变成 所以,（3）QR分解,对于任何长方矩阵A，都可以进行QR分解， 其中Q为正交
6、矩阵，R为上三角矩阵的初等变 换形式，即：A=QR 方程 A*X=b 变形成 QRX=b 所以 X=R(Qb) Q, R=qr(A) X=R(QB) 说明 这三种分解，在求解大型方程组时很有用。 其优点是运算速度快、可以节省磁盘空间、节省内 存。,二、求齐次线性方程组的通解,在Matlab中，函数null用来求解零空间， 即满足AX=0的解空间，实际上是求出解空 间的一组基（基础解系）。 格式 z = null(A) % z的列向量为方程AX=0的正交规范基 z = null(A,r) % z的列向量是方程AX=0的有理基,例4,求解方程组的通解： 解： A=1 2 2 1;2 1 -2 -2
7、;1 -1 -4 -3; format rat %指定有理式格式输出 B=null(A,r) %求解空间的有理基 B =2 5/3-2 -4/31 00 1, C=null(A) C =831/1157 -211/72443 -535/866 543/2167 280/2593 -584/941 575/1908 1882/2533,或通过行最简形得到基：, D=rref(A) D =1.0000 0 -2.0000 -1.66670 1.0000 2.0000 1.33330 0 0 0 即可写出其基础解系（与上面结果一致）。 写出通解： syms k1 k2 X=k1*B(:,1)+k2*
8、B(:,2),X =2*k1+5/3*k2-2*k1-4/3*k2k1k2 pretty(X) %让通解表达式更加精美2 k1 + 5/3 k2 -2 k1 - 4/3 k2 k1 k2 ,三、求非齐次线性方程组的通解,非齐次线性方程组需要先判断方程组是否有 解，若有解，再去求通解。 因此，步骤为： 第一步：判断AX=b是否有解，若有解则进行第二步 第二步：求AX=b的一个特解 第三步：求AX=0的通解 第四步：AX=b的通解= AX=0的通解 +AX=b的一个特解。,Solveequation.m,function solveequation(A,b) m,n=size(A); R_A=ra
9、nk(A) B=A b; R_B=rank(B) format rat if R_A=R_B & R_A=nX=Ab,elseif R_A=R_B & R_AnX=AbC=null(A,r) else X=Equation has no solves end, A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; b=1 4 0; solveequation(A,b) R_A =2 R_B =2,Warning: Rank deficient, rank = 2, tol = 8.8373e-015. In solveequation at 10 X =0 0 -8/15 3/5
10、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 =1 0 -3/2 3/4 5/4 0 1 -3/2 -7/4 -1/40 0 0 0 0,四、解方程组的其他方法（函数）,1.线性方程组的LQ解法 函数 symmlq 格式 x = symmlq(A,b) %求线性方程组AX=b的解X。A必须为n阶对称方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息；如果收敛失败，将给出
11、警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。,2.双共轭梯度法解方程组,函数 bicg 格式 x = bicg(A,b) %求线性方程组AX=b的解X。A必须为n阶方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息；如果收敛失败，将给出警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。,3.稳定双共轭梯度方法解方程组 bicgstab 4.复共轭梯度平方法解方程组 cgs 5.共轭梯度的LSQR方法 lsqr 6.广义最小残差法 qmres 7.最小残差法解方程组 minres 8.预处
12、理共轭梯度方法 pcg 9.准最小残差法解方程组 qmr,利用solve求解,函数格式： var1,var2,varN = solve(eqn1,eqn2,eqnN, var1,var2,varN) 其中eqn1，eqn2，eqnN表示代数方程，可以是线性方程，也可以是非线性方程。方程的系数可以是常数，也可以是符号常量；var1,var2,varN表示方程中的未知变量,例：求解,x =-12 y =2-1,x,y=solve(x2+y2=5,x+y=1,x,y),求解：,x1,x2,x3=solve(x1+2*x2+3*x3=5,2*x1+4*x2+x3=6,4*x1+6*x2+7*x3=15,x1,x2,x3),x1 = 8/5 x2 =1/2 x3 =4/5,
