Matlab 7.0 Simulink 6.0应用实例仿真与高效算法开发(第四章).doc

第四章 线性代数方程组的数值实验矩阵计算是求解线性代数方程组最简单有效的方法。经典的线性代数教材中，对于矩阵运算都是基于手工推导的方法，为实现高阶矩阵的分析与计算，人们引入了计算机数学语言，更方便于求解高阶问题。随着对线性代数问题的研究逐步深入，为了研究线性代数问题的数值解法与解析解问题，科学家开发了许多计算机数学语言，如，MATLAB、EISPACK、Mathmetica和Maple等，这些都大大拓宽了人们解决大型工程问题的思路。随着计算机技术的发展，人们已经不再满足于解决矩阵分析与运算问题的数值线性代数方法，逐渐的也可以求解解析解问题。本章向读者介绍了求解线性代数方程组的一些基本思路，同时也对于一些典型的代数方程组给出了求解方法。本章的内容分以下几个部分： 特殊矩阵的输入 矩阵基本分析 线性代数方程的求解 稀疏矩阵的线性方程4.1 特殊矩阵的输入 4.1.1 数值矩阵的输入.在输入一个数值矩阵时，我们可以用最底层、最基本的赋值语句逐行输入，不过对于非常复杂和具有特殊结构的矩阵来说就非常繁琐了。例如输入一个1515的单位矩阵，用一般的赋值语句也未尝不可，不过该方法对于此类特殊形式的矩阵输入是既耗时又费力。Matlab中提供了现成的函数eye( )，很轻松的输入该矩阵。所以熟悉一些特殊矩阵的输入方法是很有必要的。通常我们所谓的特殊矩阵包括：零矩阵、幺矩阵、单位矩阵、对角元素矩阵、Hankel矩阵、Hilbert矩阵及其逆矩阵、Vandermonde矩阵、伴随矩阵和随机元素矩阵等。下面我们就介绍这些特殊矩阵的输入方法。4.1.1.1 零矩阵、幺矩阵、单位矩阵及对角元素矩阵矩阵理论中，所有元素为0的矩阵定义为零矩阵，而所有元素为1的矩阵为幺矩阵。单位矩阵的定义为主对角元素均为1，其他元素皆为0。可以说对角元素矩阵更为一般的矩阵类型为对角元素矩阵，它的定义是，主对角元素可为0或非0，而非对角元素的值均为0。Matlab中提供了这些特殊矩阵的输入指令，分别举例说明。【例4-1】 分别运用matlab提供的指令生成55的方阵，使该方阵分别为零矩阵、幺矩阵和单位矩阵。 A = zeros(5) %55零矩阵输入A = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B=ones(5)B = %55幺矩阵输入 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 C=eye(5) %55单位矩阵输入C = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1指的注意的是，以上命令所完成的功能事实上可用更标准的格式来实现，分别为zeros(5,5)、ones(5,5)和eye(5,5)。 指令zeros( )和ones( ) 还可用于多维数组的生成，如要生成一个234元素皆为1的多维数组可用命令ones(2,3,4)来完成。Matlab提供的对角矩阵生成命令为diag( )，该指令的调用函数为M = diag(X) 由已知向量生成对角矩阵X = diag(M) 由已知矩阵得到对角元素列向量M = diag(X, k) 生成的矩阵M的第k条对角线上的元素列向量为X，其余元素都为0。具体指令的操作，举例如下，【例 4-2】diag( )指令的操作。 A = 3 5 9; M=diag(A)M = 3 0 0 0 5 0 0 0 9输出的M矩阵的主对角线元素为向量A。 G = M+1G = 4 1 1 1 6 1 1 1 10A = diag(G)A = 4 610矩阵G为矩阵M每个元素加1得到的新矩阵，A列向量为矩阵G的对角线元素矩阵。B = diag(G, 1)B= 1 1所得向量B为矩阵G的第1条对角线元素。 C= 2 4 3; H = diag(C,2)H = 0 0 2 0 0 0 0 0 4 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0所得到的H矩阵为以C向量为第二条对角线元素其余元素为零的55矩阵。综合应用diag( )指令可以很方便的得到比较复杂的三对角矩阵。 I = diag(1 2 3 4 5) + diag( 2 3 4 5, 1) + diag(6 5 4 3, -1)I = 1 2 0 0 0 6 2 3 0 0 0 5 3 4 0 0 0 4 4 5 0 0 0 3 54.1.1.1 Hankel矩阵Hankel矩阵的通常的表达形式为：如要构造无穷型Hankel矩阵，可令。很明显，Hankel矩阵是对称矩阵，反对角线上所有的元素都相同。在matlab中， Hankel矩阵生成指令为：指令中A和R为向量，如H矩阵为阶矩阵，那么A为n阶向量，是H矩阵的第一列各个元素所组成的向量。R为m阶向量，为最后一行元素所构成的向量。【例 4-3】用matlab语句输入下面给出的Hankel矩阵H。由Hankel矩阵的生成指令，可以知道，如通过将A向量，R向量分别设置为该H矩阵的第一列和最后一行元素组成的向量，然后用下面语句很方便的得到H矩阵。 A = 1 2 3 4;R = 4 5 6 7 8 9 10;H = hankel(A,R)H = 1 2 3 4 5 6 7 2 3 4 5 6 7 8 3 4 5 6 7 8 9 4 5 6 7 8 9 10这样就很方便的得到了H矩阵。4.1.1.2随机元素矩阵很显然，随机元素矩阵的各个元素是随机产生的。假设矩阵的元素满足0，1区间上的均匀分布，那么可以由matlab指令rand( )生成，其调用格式为： 生成nn阶标准均匀分布为随机数矩阵 生成nm阶标准均匀分布随机数矩阵当然，函数rand( )同样可用于生成多维数组。组成随机元素矩阵的伪随机数，是通过某种数学公式生成的，可满足一定随机标准的数据。这些随机数与通过电子方法获得的不可重复的随机数是不同的，它们是可以重复的。此外，如果想得到一些更具一般性的随机元素矩阵，可首先用rand( )函数得到随机矩阵A ，然后以该矩阵为变量作线性函数，比如，可通过该形式的函数来得到满足要求的矩阵。举例说明【例 4-4】 试构造一个满足，其中，的5阶正态分布随机矩阵。 A = rand(5,5)A = 0.9501 0.7621 0.6154 0.4057 0.0579 0.2311 0.4565 0.7919 0.9355 0.3529 0.6068 0.0185 0.9218 0.9169 0.8132 0.4860 0.8214 0.7382 0.4103 0.0099 0.8913 0.4447 0.1763 0.8936 0.1389 V = 0.3+2*AV = 2.2003 1.8242 1.5309 1.1114 0.4158 0.7623 1.2129 1.8839 2.1709 1.0057 1.5137 0.3370 2.1436 2.1338 1.9263 1.2720 1.9428 1.7764 1.1205 0.3197 2.0826 1.1894 0.6525 2.0873 0.5778A矩阵就是55阶标准均匀分布矩阵，经过线性组合所得到的V矩阵就是满足要求的正态随机分布矩阵。4.1.2 符号矩阵的输入 对于已经存在的数值矩阵M，可由sym( )指令转换为符号矩阵，当然，这是最基本的生成方法，这样一来对得到的新矩阵，可利用符号运算工具箱获得更高精度的解。该函数的调用格式如下:对于具有特殊形式的矩阵，如Hankel矩阵及其伴随矩阵等，符号工具箱不直接支持它们。我们可以自行编写一些函数来满足我们对此类特殊函数的要求。我们可参考matlab语言中对数字矩阵生成的相应函数，写出适合符号运算的新函数。接下来就用伴随矩阵的生成函数accompan( )来说明。Function N = accompany(Vector)Vector = Vector(:). %如Vector为列向量，那么转换为行向量n = size(Vector); %得到Vector向量的长度N = diag(ones(1, n-2), -1); %生成n-1阶第-1条对角线元素都为1的矩阵N = sym(A); %符号化 N(1, :) = -Vector(2: n)./Vector(1);【例 4-5】应用函数accompany( )生成以下多项式的伴随矩阵程序如下： syms c1 c2 c3 c4 c5 c6 c7 c8 c9 c10; N = accompan(c1 c2 c3 c4 c5 c6 c7 c8 c9 c10)N = -c2/c1, -c3/c1, -c4/c1, -c5/c1, -c6/c1, -c7/c1, -c8/c1, -c9/c1, -c10/c1 1, 0, 0, 0, 0, 0, 0, 0, 0 0, 1, 0, 0, 0, 0, 0, 0, 0 0, 0, 1, 0, 0, 0, 0, 0, 0 0, 0, 0, 1, 0, 0, 0, 0, 0 0, 0, 0, 0, 1, 0, 0, 0, 0 0, 0, 0, 0, 0, 1, 0, 0, 0 0, 0, 0, 0, 0, 0, 1, 0, 0 0, 0, 0, 0, 0, 0, 0, 1, 04.2 矩阵基本分析 4.2.1 矩阵的特征值与特征向量 矩阵A与向量x相乘，即表示矩阵对向量的变换，这里所指的变换通常为旋转、 反射、放大和缩小。不过对于任何矩阵，总存在一些特殊的向量，在其作用下，矩阵的长短发生变换而方向不变。此类向量就是所谓的特征向量，它满足的方程为：式中称为特征值，该方程为特征方程。通常求解上式的数值计算方法是：首先对矩阵A施行一系列的Household变换，产生一个准上三角矩阵，然后运用QR方法迭代，使得前面所得到的准上三角矩阵对角化。Matlab中提供了计算特征值和特征向量的指令，如下d = eig(A) %求矩阵A的特征值d，以向量形式存放d。V, D = eig(A) %计算A的特征值对角阵D和 特征向量V，使AV=VD成立。V, D = eig(A, nobalance) %当矩阵A中有与截断误差数量级相差不远的值时，该指令可能更精确。nobalance起误差调节作用。d = eig(A,B) %A、B为方阵，求广义特征值d，以向量形式存放d。V, D = eig(A, B) %计算广义特征值向量阵V和广义特征值阵D，满足AV=BVD【例 4-6】求矩阵 的特征值和特征向量。 A = 23 12 34 5; 14 35 56 15;16 18 32 43;12 37 83 120; V,D= eig(A)V = -0.1359 -0.5012 -0.4703 -0.4189 -0.2653 -0.7198 0.8587 -0.6140 -0.3414 -0.0656 -0.1671 0.6395 -0.8914 0.4757 -0.1161 -0.1964D = 164.6305 0 0 0 0 39.9348 0 0 0 0 14.4058 0 0 0 0 -8.97114.2.2 矩阵的LU分解 矩阵的三角分解又称LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。LU分解是用Gaussian主元消去法进行的。其中L和U矩阵可以分别写成， 在matlab中，实现LU分解的指令调用格式是：L,U = lu(X) %U为上三角阵，L为下三角阵或其变换形式满足LU=X。L,U,P = lu(X) %U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，满足LU=PX。举例说明函数的使用，【例 4-7】试对矩阵进行LU分解。 X = 2 5 57 23;13 25 34 12; 78 23 46 26;17 19 34 56; L1,U1=lu(X)L1 =U1 =很明显，这样得到的L1矩阵并非下三角矩阵，这是由于分解过程中采用了主元素交换的方法。现在采用另一种调用方式： L,U, P=lu(X)L =U =P = 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1此处得到的P矩阵并不是一个单位矩阵，而是单位矩阵的置换矩阵。，表明将L1矩阵的第三行换到第一行，这样就可以得到一个真正的下三角矩阵。代入检验： inv(P)*L*Uans = 2.0000 5.0000 57.0000 23.0000 13.0000 25.0000 34.0000 12.0000 78.0000 23.0000 46.0000 26.0000 17.0000 19.0000 34.0000 56.00004.2.3 矩阵的QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。matlab提供的函数为qr( ) Q,R = qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。Q,R,E = qr(A) %求得正交矩阵Q和上三角E为单位矩阵的变换形式，R的对角线元素按大小降序排列，满足AE=QR。【例 4-8】试对矩阵A = 1 2 3; 4 5 6; 7 8 9; 10 11 12;Q,R = qr(A)Q = -0.0776 -0.8331 0.5444 0.0605 -0.3105 -0.4512 -0.7709 0.3251 -0.5433 -0.0694 -0.0913 -0.8317 -0.7762 0.3124 0.3178 0.4461R = -12.8841 -14.5916 -16.2992 0 -1.0413 -2.0826 0 0 0.0000 0 0 04.2.4 矩阵的奇异值分解 矩阵的奇异值可以看作是矩阵的一种测度。对任意的阶矩阵A，总会有并且，同时，和有相同的非负特征值，在数学上把这些非负特征值的平方根称作矩阵A的奇异值，记作，。矩阵奇异值分解的定义为，对于任意矩阵，存在酋阵，使，其中，。以上，的组合就称矩阵A的奇异值分解三对组，它们分别为矩阵A的第i个奇异值、左奇异值和右奇异值向量。从理论上讲，矩阵A奇异值恰等于ATA(或AAT)特征值的平方根，不过从数值上看，借助ATA(或AAT)特征值求取矩阵A奇异值的方法是不可取的，将丧失计算精度。看下面例子，考虑一个演示例子，假设矩阵A为，其中，求取矩阵A的秩。 A = 1 1; 6*eps 0; 0 6*eps; rank(A)ans = 2同时又注意到， rank(A*A)ans =1即ATA的秩为1，为什么会这样呢，我们可以得到，双精度运算中，由于为量级，如果在加上1，事实上就已经被忽略了，所以ATA矩阵将退化成幺矩阵，其秩为1。这与实际矛盾，故对这样的问题应该引入一个新的量做为矩阵的测度，这样引入奇异值的概念就显得很必要了。Matlab提供了直接求取矩阵奇异值分解的函数svd( )，其调用方式为D=svd(A) %只是计算矩阵的奇异值L,B,M=svd(A) %矩阵奇异值与变换矩阵其中，A为原始矩阵，返回的B为对角矩阵，而L和M均为正交变换矩阵，并满足。矩阵的奇异值大小通常决定矩阵的形态，如果矩阵的奇异值变化特别大，则矩阵中某个元素有一个很小的变化会严重影响到原矩阵的参数，又可将其称为病态矩阵。另一个比较重要的概念是矩阵的条件数，记作cond(A)，即条件数的大小，反映着对元素变化的敏感程度。Matlab中提供的求矩阵A的条件数的函数为cond(A)。下面举例对给定矩阵进行奇异值分解，【例 4-9】 对矩阵进行奇异值分解。 A = 17 2 4 7; 2 57 35 7; 28 93 21 20; 23 21 13 15; L, B, M=svd(A)L =-7.004381577366817e-0025.231149103298072e-0013.351701456897706e-001-7.804521945005988e-001-5.163440830512738e-001-6.124153900427823e-0015.881225200750654e-001-1.115709609048906e-001-8.144830700530317e-0011.649159889448149e-001-5.536165869539846e-001-5.411764810819996e-002-2.551469779691531e-0015.692986011046732e-0014.850582949175016e-0016.127948865201114e-001B =1.236690421626503e+00200002.974689953424784e+00100001.843184260490569e+00100002.087380382198967e+000M =-2.498389164233517e-0018.531859995077364e-0011.372198059348684e-001-4.368350843780604e-001-8.949427478559386e-001-2.208293946449815e-001-3.855715681894596e-001-4.057612998690056e-002-3.135244898525784e-001-2.850031956512153e-0019.008745706661361e-001-9.434288875363396e-002-1.958580816896453e-0013.769366089587355e-0011.446730748626815e-0018.936599499100979e-001对于长方形矩阵，也可用svd( ) 命令进行奇异值分解。【例 4-10】对以下长方形矩阵进行奇异值分解。 B=1 2 3;4 5 6;7 8 9;10 11 12;13 14 15 L,B1,M=svd(B)L =-1.013459963216538e-0017.679381414082588e-001-1.826540898671791e-002-4.182600740172534e-001-4.740515639860605e-001-2.485687533216272e-0014.880712805237905e-0015.366924394418480e-0012.763055233190366e-0015.793241607505693e-001-3.957915103216006e-0012.082044196393209e-001-8.133012257039489e-0011.110365817377026e-0013.552632739032574e-001-5.430142673215739e-001-7.166244124514866e-0028.958676902922533e-0026.220505626364989e-001-5.522927741139834e-001-6.902370243215471e-001-3.515293021296176e-0012.052874262195935e-001-5.911325936759846e-0019.175690344621651e-002B1 =3.518264833189422e+0010001.476907699980091e+0000009.860230903782525e-016000000M =-5.192726084435612e-001-7.507924423258088e-001-4.082482904638638e-001-5.755207206646645e-001-4.592639131217075e-002 8.164965809277259e-001-6.317688328857681e-001 6.589396597014696e-001-4.082482904638623e-001接下来验证一下运算的结果： B2=L*B1*MB2 =9.999999999999996e-0011.999999999999997e+0002.999999999999998e+0004.000000000000001e+0005.000000000000000e+0006.000000000000004e+0007.000000000000002e+0007.999999999999999e+0009.000000000000004e+0001.000000000000000e+0011.100000000000000e+0011.200000000000001e+0011.300000000000000e+0011.400000000000000e+0011.500000000000000e+001由于运算过程中的舍入误差，导致所得到的验证矩阵与原矩阵有很微小的误差。4.3 线性代数方程的求解 许多科学与工程问题都归结为一个线性代数方程组，例如，曲线拟合所得到的正规方程组，三次样条中常需要求解一个三对角线性方程组以及有限元离散问题中大型线性方程组等等。正因为线性方程组有着广泛的应用，所以对它的数值求解方法的研究一直都没有停止。本节针对常见的几类线性代数方程给出它的MATLAB解法。 4.3.1 线性代数方程组 考虑线性方程组的一般形式 （4.3.1）其中， （4.3.2）由此得到增广矩阵 （4.3.3）一般地，线性代数方程组的求解根据其解的存在定理可以分为以下三类：4.3.1.1当，且时，方程组（4.3.1）有唯一解 （4.3.4）这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 直接法；另一类是解大型稀疏矩阵 迭代法(将在4.4节里详细介绍)。利用矩阵除法求线性方程组的唯一解，语法如下：方程：AX=b解法：X=Ab【例 4-11】求方程组的解。方法1：直接利用矩阵除法，运行后结果如下 A=5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5; B=1 0 0 0 1; R_A=rank(A) %求秩R_A = 5 X1=AB %求解X1 = 2.2662 -1.7218 1.0571 -0.5940 0.3188这就是方程组的解。将上面得出的解代入原方程，误差为： Error1 =norm(A*X1 - B)Error1 = 1.9984e-015可以看出误差很小，达到10e-15数量级。方法2，对于该问题也可以利用函数rref求解： C=A,B %由系数矩阵和常数列构成增广矩阵CC = 5 6 0 0 0 1 1 5 6 0 0 0 0 1 5 6 0 0 0 0 1 5 6 0 0 0 0 1 5 1 R=rref(C) %将C化成行最简行R = 1.0000 0 0 0 0 2.2662 0 1.0000 0 0 0 -1.7218 0 0 1.0000 0 0 1.0571 0 0 0 1.0000 0 -0.5940 0 0 0 0 1.0000 0.3188这样，R的最后一列元素就是所求之解。再考虑误差： X2 = R(26:30) % 取出矩阵R的最后一行；X2 = 2.2662 -1.7218 1.0571 -0.5940 0.3188这就是方程组的解。将上面得出的解代入原方程，误差为： Error2 = norm(A*X2 - B)Error2 = 4.3785e-005 显然，该方法的误差要比上述方法1的误差大很多。事实上，如果也可以采用符号运算工具箱中的求逆函数得出该方程的精确解，方法如下。方法3：利用符号运算工具箱中的求逆函数求精确解 X3 = inv(sym(A)*BX3 = 1507/665 -229/133 37/35 -79/133 212/665 Error3 =norm(double(A*X3 - B)Error3 = 0该方法得到了精确解。利用矩阵的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)【例4-12】求方程组的一个近似解。解：用LU分解方法如下： A = 4 2 -1; 3 -1 2; 11 3 0; B = 2 10 8; D = det(A)D = 0 L,U=lu(A)L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0U = 11.0000 3.0000 0 0 -1.8182 2.0000 0 0 0.0000 X=U(LB)Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.018587e-017.说明 结果中的警告是由于系数行列式为零产生的。X = 1.0e+016 * -0.4053 1.4862 1.3511即为该方程的一个近似解。 norm(A*X - B)ans =2.8284即为所得的解的误差。说明 结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。（2）Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即： 其中R为上三角阵。方程 A*X=b 变成 所以 （3）QR分解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR方程 A*X=b 变形成 QRX=b所以 X=R(Qb)【例 4-13】例4.3.3， 同上例4.3.2， Q,R = qr(A)Q = -0.3310 0.4730 -0.8165 -0.2483 -0.8785 -0.4082 -0.9104 0.0676 0.4082R = -12.0830 -3.1449 -0.1655 0 2.0272 -2.2299 0 0 0.0000 X = R(QB)Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.594136e-017.X = 1.0e+015 * 1.3238 -4.8539 -4.4126即为该方程的一个近似解。 norm(A*X - B)ans =5.5902即为所得的解的误差。说明 在求解大型方程组中经常使用这三种分解方法，其优点是运算速度快、可以节省磁盘空间、节省内存。4.3.1.2当时，方程组（4.3.1）有无穷多解可以构造出线性方程组的个化零向量，原方程对应的齐次方程组的解可以由的线性组合来表示，即 （4.3.5）其中，为任意常数。在Matlab中，函数null用来求解零空间，即满足AX=0的解空间，实际上是求出解空间的一组基（基础解系）。Z = null(A), 求解矩阵A的化零矩阵(一组基)。其中Z的列数为n-r，而各列构成的向量又称为矩阵A的基础解系。Z = null(A,r)求解矩阵A的有理基.【例 4-14】例4.3.4 求解方程组的通解：方法1：首先判定方程的可解性： A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8A = 1 1 -3 -1 3 -1 -3 4 1 5 -9 -8 B=1 4 0b = 1 4 0 C=A BC = 1 1 -3 -1 1 3 -1 -3 4 4 1 5 -9 -8 0 R_A=rank(A)R_A = 2 R_C=rank(C)R_C = 2 显然，矩阵A和C的秩相同，都等于2，小于矩阵的阶数4，由此可知原线性方程组有无穷多组解。其次，先求出化零空间的Z（对应的齐次方程的基础解系）和满足方程的一个特解X0. Z = null(A,r)Z = 3/2 -3/4 3/2 7/4 1 0 0 1 X0 = pinv(A)*BX0 = 130/371 -34/371 -144/371 157/371 最后将方程全部的解表示出来 syms k1 k2 % k1和k2为任意常数。 X = k1*Z(:,1) + k2*Z(:,2) + X0 %写出方程组的通解。X = 3/2*k1-3/4*k2+130/371 3/2*k1+7/4*k2-34/371 k1-144/371 k2+157/371还可以验证所得的解是否精确： k1 =rand(1),k2=rand(1) %取任意随机变量。k1 = 1105/1163 k2 = 337/1458 X=k1*Z(:,1) + k2*Z(:,2) + x0 X = 1285/802 763/439 825/1468 4310/6587 norm(A*X - B) % 验证所得的解的精确性。ans = 1.275549143317629e-015显然，我们得到的解还是相当精确的。方法2：用函数rref求解 A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; b=1 4 0; C=A B; C_simple = rref(C) %求增广矩阵的行最简形，可得最简同解方程组。C_simple = 1 0 -3/2 3/4 5/4 0 1 -3/2 -7/4 -1/4 0 0 0 0 0 对应齐次方程组的基础解系为：， 非齐次方程组的特解为：所以，原方程组的通解为：X=k1+k2+。4.3.1.3 当时，方程组（4.3.1）无解此方程称为矛盾方程。矛盾方程可以利用Moore-Penrose广义逆求解出方程的最小二乘解： （4.3.6）该解不满足原方程，只能得到在范数的意义下使误差最小的近似解。【例 4-15】求解方程组：的解解： A = 1 2 3 4;2 2 1 1;2 4 6 8;4 4 2 2A = 1 2 3 4 2 2 1 1 2 4 6 8 4 4 2 2 B = 1;2;3;4B = 1 2 3 4 C = A BC = 1 2 3 4 1 2 2 1 1 2 2 4 6 8 3 4 4 2 2 4 rank(A),rank(C)ans = 2 3 由于，故原方程为矛盾方程，无解；只能利用函数pinv()求取Moore-Penrose广义逆，从而得到原方程的最小二乘解。 x = pinv(A)*B % x = 0.54656488549618 0.45496183206107 0.0