第四讲：应用MATLAB解决高等代数问题 1.交换矩阵中的两个行向量的位置； 2.用一个非零数乘以矩阵的某一行向量 3.把矩阵的某一个行向量乘以实数并加到矩阵的另一行上

例：求解线性方程组 线性代数方法用增广矩阵初等变换即消元法过程

经过初等行变换将矩阵 变为矩阵 这时矩阵对应的方程组 此方程组的解为

方法之一：初等变换法 A=[3 -1 5 3;1 -1 2 1;1 -2 -1 2] %输入矩阵的数据 A([1 3],:)=A([3 1],:) %交换第一行和第三行数据 A(2,:)=A(2,:)-A(1,:) %将第一行乘以-1加到第二行 A(3,:)=A(3,:)-3*A(1,:)%将第一行乘以-3加到第三行 A(3,:)=A(3,:)-5*A(2,:) %将第二行乘以-5加到第三行

A=[3 -1 5 3;1 -1 2 1;1 -2 -1 2] %输入矩阵的数据 format rat %分数数据格式 rref(A) %化简矩阵 方法之二：Cramer法则 A=[3 -1 5;1 -1 2 ;1 -2 -1] %输入矩阵的数据 B=[3 1 2]‘; %输入线性方程组的常数项 S=[0 0 0]‘; %给解向量赋初值 for i=1:3 %for循环

C=A; %将矩阵A赋给临时矩阵C C(:,i)=B; %将常数项赋给矩阵C的第i列即求Ai S(i)=det(C)/det(A); %求xi end format rat %数据格式说明为分数形式 S %显示S

二、线性方程组的解结构 方法之三：利用矩阵的左除“\” A=[3 -1 5;1 -1 2 ;1 -2 -1] ； b=[3 1 2]‘； x=A\b x = 10/7 -1/7 -2/7 二、线性方程组的解结构 1。齐次方程组的解结构AX=0 求解方程过程如下

提取方程组右端各自由求知量的系数形成的向量组即为基础解系 根据最简行阶梯形矩阵写出简化方程组 确定自由求知量 整理方程组为向量形式量 提取方程组右端各自由求知量的系数形成的向量组即为基础解系 将系数矩阵化为最简行阶梯形矩阵

例：解线性方程组： 应用MATLAB计算过程如下： A=[1 –1 1 –1;1 –1 –1 1;1 –1 –2 2] %输入矩阵 rref(A) %将矩阵化为最简阶梯形矩阵

运行结果为： A = 1 -1 1 -1 1 -1 -1 1 1 -1 -2 2 ans = 1 -1 0 0 0 0 1 -1 1 -1 1 -1 1 -1 -1 1 1 -1 -2 2 ans = 1 -1 0 0 0 0 1 -1 0 0 0 0

由运行结果知化简的等价方程组为： 取x2，x4为自由求知量，得方程组的解的向量形式为

所以基础解系为： 所以齐次方程组的通解为

2.非齐次方程组的解的结构 求解非齐次线性方程组的通解的步骤如下： 1）、写出非次方程组的增广矩阵； 2）、将增广矩阵化为最简行阶梯形矩阵； 3）、观察增广矩阵与系数矩阵的秩是否相等，若相等则方程组有解，若不相等则方程组无解； 4）、写出对应的简化的线性方程组； 5）、确定自由求知量 6）、整理方程组为向量形式。

例：求解下列非齐线性方程组

A=[1 2 3 1;1 4 5 2;2 9 8 3;3 7 7 2]; b=[3;2;7;12] ； format rat 在MATLAB中输入的命令如下 A=[1 2 3 1;1 4 5 2;2 9 8 3;3 7 7 2]; b=[3;2;7;12] ； format rat c=[A b]； rref(c); 计算结果如下

ans = 1 0 0 -1/2 31/6 0 1 0 0 2/3 0 0 1 1/2 -7/6 0 0 0 0 0 所以简化方程组为：

取x4为自由求知量 所以原线性方程组的通解为：

三、向量组的线性相关性判定 1.向量组线性相关与线性无关的定义： 设有m个向量 如果存在m个不全为零的一组数k1,k2, …,km使 成立，则称向量组 是线性相关的。 如果仅当k1=k2= …=km=0时

才有上面的等式成立，则称向量组线性无关 2.应用MATLAB进行向量组的线性相关性的判定步骤： 1）将给定的m个向量组的写成列向量形式，组成一个n×m阶的矩阵 2)判定是否存在不全为零的一组数k1,k2, …,km使得

即判定线性方程组

是否有非零解，从而有 这说明向量组线性相关。如果方程组只有零解，则说明该向量组线性无关。 3）用命令rref将矩阵A化为最简行阶梯形矩阵； 4）观察最简行阶梯形矩阵中非零行向量的数目是否小于向量组全部向量数目m,若小于m则向量组线性相关；否则线性无关。

例 判断下列向量组的线性相关性 1)、a1=[4 3 –1 1 –1],a2=[2 1 –3 2 –5] a3=[1 –3 0 1 –2],a4=[1 5 2 –1 6] 2）a1=[1 0 0 1 4],a2=[0 1 0 2 5],a3=[0 0 1 3 6] a4=[1 2 3 14 32],a5=[4 5 6 32 77] (北京大学数学力学系《高等代数》p151 16-4）

解：1)先在MATLAB中将上面四个向量以行向量数据形式输入，再转置为列向量组成的矩阵，然后用rref命令将其化为最简行阶梯形矩阵，命令如下 A=A' rref(A)

ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 最简行阶梯形矩阵的变量名为ans,它的不全为零为行向量数为4，而向量组中向量数也是4,所以向量组是线性无关的。

2)可以应用矩阵拼接命令 a1=[1 0 0 1 4],a2=[0 1 0 2 5],a3=[0 0 1 3 6] a4=[1 2 3 14 32],a5=[4 5 6 32 77] A=[a1;a2;a3;a4;a5] rref(A)

ans = 1 0 0 1 4 0 1 0 2 5 0 0 1 3 6 0 0 0 0 0 得非零行数为3，所以该向量组线性相关。

四、向量组的最大无关组 1.极大无关组的定义：对于一个相关向量组T中最多有多少个向量是线性无关的，这就是极大无关组， 即一向量组的一个部分组本身是线性无关的，并且从这向量组中任意添加一个向量（如果还有的话），所得的部分向量组都线性相关的。

2.秩的定义：极大线性无关组所含向量个数r称为向量组的秩。 3.应用MATLAB求向量组的极大无关组的方法 借助向量组线性相关性分析的方法，可得求向量组 的极大无关组的步骤如下： 1）.将向量组中每个向量以列的形式排成矩阵A=[a1 a2…am]

2).把矩阵A化为最简行阶形矩阵 3).确定最简行阶梯形矩阵中非零行向量数目r（即向量组T的秩），在最简行阶梯形矩阵中寻找r个无关的列向量 4).根据 所在位置确定矩阵A中列向量位置即得T的极大无关组 在最简行阶梯形矩阵中寻找r个线性无关的列向量

时，只须在仅有一个非零元素的列向量中 寻找，非零元素不在同一位置的这类向量是线性无关的。 例：求下列向量组的秩和一个极大无关组并将其余向量用极大无关组线性表出 1）a1=[1 2 1 3];a2=[4 -1 -5 -6];a3=[1 -3 -4 -7];a4=[2 -1 1 0];A=[a1;a2;a3;a4]

a3=[3;0;7;14]; a4=[1;-1;2;0];a5=[2;1;5;6]; 北大《高等代数》P151 9-2 解：1）输入向量及命令如下 a1=[1 2 1 3];a2=[4 -1 -5 -6];a3=[1 -3 -4 -7]; a4=[2 -1 1 0]; A=[a1;a2;a3;a4] A=A' rref(A)

得简化的行阶梯形矩阵为 ans = 1 0 -11/9 0 0 1 5/9 0 0 0 0 1 0 0 0 0

最简矩阵中的有三个不全为零的行向量，所以向量组的秩为3，显然第一列、第二列、第四列线性无关，所以对应于原向量一个极大无关组为a1,a2 a4，最简矩阵中第三列向量有两个非零元素-11/9，5/9，它们是方程组 的解（x1=-11/9,x2=5/9)，也是方程组

的解，所以 被最大无关组 线性表出 2）输入向量及命令如下：

a1=[1;-1;2;4];a2=[0;3;1;2];a3=[3;0;7;14]; a4=[1;-1;2;0];a5=[2;1;5;6]; A=[a1 a2 a3 a4 a5] rref(A) 得最简行阶梯形矩阵 ans = 1 0 3 0 1 0 1 1 0 1 0 0 0 1 1 0 0 0 0 0

由此可知向量组的秩为3，第1列，第2列，第4列的向量是线性无关的，所以a1,a2,a4是极大无关组。最简矩阵中第三列向量有两个非零元素3,1，它们是方程组 的解（x1=3,x2=1)，也是方程组

的解，所以 被最大无关组 线性表出 最简矩阵中第五列向量有三个非零元素1,1,1，它们是方程组 的解(x1=1,x2=1,x4=1)，也是方程组

的解，所以 被最大无关组 线性表出 注意：两个例子中输入的向量和命令有所不同，请同学们思考为什么？

五.矩阵的特征值和特征向量 1 矩阵的特征值和特征向量 设A是n阶方阵，k是一个数，如果存在一非零的列向量X使得AX=kX成立，则称数k为A的征值，非零列向量X称为方阵A的属于特征值K的一个特征向量。 用MATLAB的命令 eig可以求出矩阵A的特征值和特征向量的方法有两种

法一）只求A的特征值命令为eig(A) 法二）同时求出特征值和特征向量用命令[p d]=eig(A) 例　求方阵 特征值和特征向量。

解：先输入矩阵的数据，然后用eig的两种使用方法求解，命令如下 A=[3 0 4;0 6 0;4 0 3]; eig(A) [p d]=eig(A) ans = -1 6 7 第一个命令eig(A)的结果为

命令［p d]=eig(A)计算结果为 p = 0.7071 0 0.7071 0 -1.0000 0 -0.7071 0 0.7071 d = -1 0 0 0 6 0 0 0 7

北京大学《高等代数》Ｐ２８４ 求矩阵 应用eig(A)得 ans = -1 5

应用[p d]=eig(A)结果为 p= 2131/3543 709/1284 780/1351 2131/3543 709/1284 780/1351 408/2299 -369/463 780/1351 -747/959 294/1201 780/1351 d = -1 0 0 0 -1 0 0 0 5

２ 矩阵的相似对角化 设三阶矩阵Ａ有三个线性无关的特征向量a1,a2,a3，对应的特征值为k1,k2,k3,现定义两个矩阵 Ｐ＝［a1 a2 a3] Aa1=k1a1,Aa2=k2a2,Aa3=k3a3 矩阵形式为AP=P∧或Ａ＝Ｐ∧Ｐ-1

说明矩阵A与对角矩阵相似。利用特征矩阵向量和特征值的方法可以求矩阵A的相似对角矩阵。矩阵的相似对角化方法可用计算一矩阵的方幂。 例设矩阵 求A10 解：先求出A的特征值和特征向量，得A的对角相似矩阵∧和可逆矩阵P，由等式Ａ＝Ｐ∧Ｐ-1

得A10=Ｐ∧10Ｐ-1 MATLAB命令如下 A=[4 6 0;-3 -5 0;-3 -6 1];[p d]=eig(A) p*d^10*inv(p) 结果为 ans = -1022 -2046 0 1023 2047 0 1023 2046 1

为了验证这一结果也可以直接输入命令A ^10也得这一结果。 例 判断二次型的类型（正定型、负定型、半正定型、半负定型），将结果化为标准形式。 解：首先写出二次型的矩阵，然后求特征值，由特征值的符号判断二次型的类型， 根据二次型的系数得其矩阵为

在MATLAB中输入矩阵A的数据并求特征值，所用命令如下： eig(A)

ans = -8 -5 -2 计算结果为： 说明A有三个负特征值，所以该二次型为负定二次型，它的标准形式为： 为了求得其变换矩阵C的数据，可由命令 [p d]=eig(A)得矩阵A的特征向量矩阵

显然三个列向量相互正交的单位向量，可得变量的变换关系为 p= -2/3 1/3 2/3 2/3 2/3 1/3 1/3 -2/3 2/3 显然三个列向量相互正交的单位向量，可得变量的变换关系为

由此可以用此方法求《高等代数》二次型的变换矩阵、化简二次型及二次型的正定判断。

六、应用线性方程组求解数学模型 1.实际问题中方程组的类型： 适定方程组、不定方程组、超定方程组 适定方程组：方程数等于未知数个数时，这一类方程组称为适定方程组。如果其系数矩阵可逆，适定方程组有唯一的解，求解方法有克菜姆方法、消元法、矩阵分解法、迭代法。

不定方程组：实际方程的个数小于未知数的个数，这一类方程称为不定方程组。当系数矩阵的秩等于增广矩阵的秩时，不定方程组有无穷多组解，根据线性代数的理论和方法，可求得方程组的通解。 超定方程组：当方程组的数多于未知数的个数时，这一类方程组称为超定方程组。超定方程组没有准确解，但可以求广义解，例超定方程组的最小二乘解。

1)适定方程组：对于方程组AX=b,如果A为方阵，解适定方程组可以用方阵的系数矩阵的逆来求，即x=inv(A)\b 如果A是奇异方阵，则计算结果为INF， 则会给出警告信息。如果A为病态矩阵，也会出出警告信息。

此外，还可以用除法来解适定方程，X=A\B A=rand(100)+1.e10; %生成一个100阶的随机矩阵 x=ones(100,1); b=A*x; %求方程右边的向量

tic %开始计时 y=inv(A)*b; %用逆运算求解方程 toc %读计时时间 err=norm(y-x) %结果与精确解的2范数 res=norm(A*y-b) %方程的2范数误差 y=A\b; %除法运算求解方程

elapsed_time = err = 0.0457 res = 9.7113e+008 0.0360 0.0033

2）超定方程 对于方程Ax=b,A为n,m矩阵，如果A列满秩，且n>m,则方程没有精确解的，即其精确解的空间为零。然而在实际工程计算时，求得最小二乘解也是有意义的，方程的解除了用除法运算来求（x=A\b)外，还可以用广义逆来求：x=pinv(A),所求的解并不满足Ax=b，而x只是最小二乘意义上的解。 A=[3 4 5;6 1 2;4 -5 7;0 2 4]

B=[3 2 4 6]' x1=A\B A*x1-B x1 = 0.4149 0.0448 0.3737 ans = 0.2924 1.2815 0.0516 -1.0966

由此可见，x1不是方程Ax=B精确解，用x2=pinv(A) 由此可见，x1不是方程Ax=B精确解，用x2=pinv(A)*B所得的解与x1相同，用线性代数可以证明，列满秩的方程组Ax=B最小二乘解为X=inv(C’*C)*C*B,而广义逆pinv(A)=inv(C’*C)*C，如上例的结果可以这样计算 inv(A'*A)*A'*B ans = 0.4149 0.0448 0.3737

3)不定方程：理论上说有无穷多个解，如果用逆矩阵法和除法求解时，只能得到其中的一个解 A=[ 1 -2 3;0 1 –1;-1 0 –1;1 -3 4] B=[4 -3 -4 1]’ x=pinv(A)*B x= 2.2549 1.2157 1.0392

Warning: Rank deficient, rank = 2 tol = 4.6151e-015. y = 3.4706 x=A\B Warning: Rank deficient, rank = 2 tol = 4.6151e-015. y = 3.4706 -0.1765 x和y都是方程组的解，其中x=pinv(A)*B是方程解中最小的一个，norm(x)=2.7645.而norm(y)=3.4751,y=A\B是所有解中0最多的一个，也就是非零元素最多的一个。

2、交通流量问题 1）问题提出 如图所示给出了某城市部分单行街道的交通流量（每小时通过的车辆数）。图中有6个路口，已有9条路口记录了当天的平均车流量，另有7处的平均车流量未知，试用每个路口的进出车流量相等关系推算7处的平均车流量

在图中的任何一个路口（十字路口或丁字路口）处，都有车辆流进和流出。当一天结束后，流进的流出的车辆数应该相等以达到平衡， x1 x5 200 300 400 500 x2 x3 x4 x6 x7 2）问题分析与数学模型 在图中的任何一个路口（十字路口或丁字路口）处，都有车辆流进和流出。当一天结束后，流进的流出的车辆数应该相等以达到平衡，

在图中有的长远街道车流量有数据记录，而有的没有数据记录，我们可以理解为有数据的街道有专人（或者设备）记录了当天的车流量情况，而没有记录的街道由于人力不足（或设备的经费还没到位）造成的。为了填补空白，设在没有数据的街道处假设车流量是未知数，在每一个路口可根据进出的车流量相等关系，建立一个线性方程，如图有六个路口，可以建立六个方程的线性组，问题的答案应该是在所列的线性方程组

通解中支寻找，将方程组写成矩阵向量形式为AX=b其中

显然是一个不定方程组，因为方程组的个数少于未知数的个数。所以当方程组的系数矩阵A的秩增广矩阵[A b]的秩相等时，该问题有无穷多解，由于 图

中的街道是单行道，每一街道上的车流量只能是正数或零，故应在方程组的解集合中寻找非负解，如果方程组没有解或没有非负解，说明问题所给的数据有误，求解问题分三步，第一步判断方程组是否有解，第二步，如果方程组有解则求出方程组的通解，第三步，在通解中找非负特解。 3）程序和计算结果 A=[1 0 1 0 0 0 0;1 -1 0 1 0 0 0;0 1 0 0 -1 0 0;0 0 1 0 0 1 0; 0 0 0 1 0 1 -1; 0 0 0 0 1 0 -1]

b=[700 200 200 500 0 -200]' >>rank(A) >> rank([A b]) ans =

>> rref([A b]) ans = 1 0 0 0 0 -1 0 200 0 1 0 0 0 0 -1 0 0 0 1 0 0 1 0 500 0 0 0 1 0 1 -1 0 0 0 0 0 1 0 -1 -200

即有简化方程组 x6、x7为自由未知量，直接可得原方程组的通解形式

由上面所得的方程组通解表达式，取适当的k1和k2使特解为非负数即可。求非负解的程序如下：

A=[1 0 1 0 0 0 0;1 -1 0 1 0 0 0;0 1 0 0 -1 0 0;0 0 1 0 0 1 0; 0 0 0 1 0 1 -1; 0 0 0 0 1 0 -1] ; b=[700 200 200 500 0 -200]'; x=zeros(7,1); s=rref([A b]); s=[s;zeros(1,8)]; k1=0;

while 1==1 k2=200; x=k1*((-1)*s(:,6)+[0 0 0 0 0 1 0]')+k2*((-1)*s(:,7)+[0 0 0 0 0 0 1]')+s(:,8); k2=k2+1; for i=1:7 if x(i)<0

break; else continue; end

end x=k1*((-1)*s(:,6)+[0 0 0 0 0 1 0]')+k2*((-1)*s(:,7)+[0 0 0 0 0 0 1]')+s(:,8) 事实上，可以通过前面的通解公式可得当k1=0,k2=200时，得 A=[1 0 1 0 0 0 0;1 -1 0 1 0 0 0;0 1 0 0 -1 0 0;0 0 1 0 0 1 0; 0 0 0 1 0 1 -1; 0 0 0 0 1 0 -1]

b=[700 200 200 500 0 -200]' x=zeros(7,1); s=rref([A b]); s=[s;zeros(1,8)] x1=(-1)*s(:,6)+[0 0 0 0 0 1 0]' x2=(-1)*s(:,7)+[0 0 0 0 0 0 1]' nonegativs=0*x1+200*x2+s(:,8)

nonegativs =[ 200 200 500 200 0 0 200 200 300 400 500

闭合经济问题 一个木工，一个电工，一个油漆工，三人相互同意彼此装修他们自己的房子。在装修之前，他们达成如下协议：（1）每人总共工作十天（包括给自己家干活在内）；（2）每人的日工资根据一般的市价在60~80元之间；（3）每人的日工资总数应使得每人总收入与总支出相等。下表是他们协议制定出的工作天数的分配方案