第4章 线性代数 4.1 矩阵的生成 4.1.1 通过元素列表榆入 4.1.2 通过外部数据加载 4.1.3 在M文件中创建矩阵 4.1.4 通过函数产生矩阵 线性代数中有若干特殊意义的矩阵，在Matlab中可以很容易的通过函数的方式创建它们，见表4.1。其中：matlab函数名必须小写。

产生均值为0，方差为1的标准正态分布随机矩阵 魔方矩阵 magic(n) 生成一个n阶魔方阵（其每行、每列及两条对角线上的元素和都相等） 特殊矩阵名称 函数命令格式 功能简介 全零矩阵 zeros(m,n) zeros(n) zeros(size(A)) 创建m行n列的全零矩阵； 创建n行n列的全零矩阵； 产生一个与矩阵A同样大小的零矩阵 全1矩阵 ones(m,n) 创建m行n列的全1矩阵 单位矩阵 eye(m,n) 创建m行n列的对角线为1的矩阵 随机矩阵 rand(m,n) 创建m行n列的(0,1)均匀分布的随机矩阵 正态分布随机矩阵 randn(m,n) 产生均值为0，方差为1的标准正态分布随机矩阵 魔方矩阵 magic(n) 生成一个n阶魔方阵（其每行、每列及两条对角线上的元素和都相等） Vandermonde 矩阵 vander(V) 生成以向量V为基础向量的范得蒙矩阵（矩阵最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积） Hilbert矩阵 hilb(n) 生成希尔伯特矩阵（第(i,j)个元素为1/(i+j-1)） Hilbert逆矩阵 invhilb(n) 求n阶的希尔伯特矩阵的逆矩阵 Toeplitz矩阵 toeplitz(x,y) 生成一个以向量x为第一列，向量y为第一行的托普利兹矩阵（除第一行第一列外，其他每个元素都与左上角的元素相同） Hankel 矩阵 hankel(c) hankel(c,r) 生成一个以向量c为第一列，其它元素是0； 生成一个以向量c为第一列，最后一行是r向量，其中第一个元素以c的为准，r向量是从第二个元素开始的（汉克尔矩阵是反对角线上元素相等矩阵） 伴随矩阵 compan(p) p是一个多项式的系数向量，高次幂系数排在前，低次幂排在后 帕斯卡矩阵 pascal(n) 生成一个n阶帕斯卡矩阵（由杨辉三角形表组成的矩阵） Hadamard矩阵 hadamard(n) n，n/12或，n/20必须是2的乘方。 （哈达玛矩阵是由+1和-1元素构成的正交方阵） 空矩阵 [] 创建一个空矩阵

示例4.1 求特殊矩阵 (1) 在区间[20,50]内均匀分布的5阶随机矩阵。 x =20+(50-20)*rand(5) (2) 均值为0.6、方差为0.1的5阶正态分布随机矩阵。 y=0.6+sqrt(0.1)*randn(5) (3) 求4阶希尔伯特矩阵及其逆矩阵。 format rat %以有理形式输出 H=hilb(4) H=invhilb(4)

(4) 已知c = 1:3; r = 7:10;求其Hankel 矩阵。 hankel(c,r) (5) 求(x+y)5的展开式。 pascal(6) 矩阵次对角线上的元素1,5,10,10,5,1即为展开式的 系数.

注：数值矩阵A转换成符号矩阵B，可以通过 whos命令来详细查看这两类矩阵的区别。 4.1.5 符号矩阵的生成 生成方法 命令格式 功能简介 直接创建 sym('[ , , , ]') 矩阵元素可以是任何的符号变量、符号表达式及方程，且元素的长度可以不同。 间接创建 syms 在创建符号矩阵之前，先把符号矩阵的 所有变量定义为符号变量，然后按创建普通矩阵的格式输入矩阵 由数值矩阵转化为符号矩阵 B=sym(A) 将一个数值矩阵A转化为符号矩阵B 注：数值矩阵A转换成符号矩阵B，可以通过 whos命令来详细查看这两类矩阵的区别。

4.2 矩阵操作 4.2.1 元素变换操作 有已知矩阵A，由A的元素构成的各种矩阵分 别如表4.3所示。

元素操作名称 操作命令格式 功能简介 对角阵 diag(A); diag(A,k); diag(V); diag(V,k) 提取矩阵A主对角线元素产生一个具有min(m,n)个元素的列向量； 提取第k条对角线的元素，主对角线为0，向上为第1，2…条，向下为-1，-2…条； 产生一个m×m对角矩阵，其主对角线元素即为向量V的元素； 产生一个n×n(n=m+|k|)对角阵，其第k条对角线的元素即为向量V的元素 （只有对角线上有非0元素的矩阵称为对角矩阵，对角线上的元素相等的对角矩阵称为数量矩阵，对角线上的元素都为1的对角矩阵称为单位矩阵。） 上三角阵 triu(A); triu(A,k) 求矩阵A的上三角阵； 求矩阵A的第k条对角线以上的元素 （对角线以下的元素全为0的一种矩阵） 下三角阵 tril(A); tril(A,k) 求矩阵A的下三角阵； 求矩阵A的第k条对角线以下的元素 （对角线以上的元素全为0的一种矩阵） 转置矩阵 A′ transpose(S) 求矩阵A的转置（行列互换） 符号矩阵S的转置矩阵 旋转矩阵 rot90(A,k) 矩阵A旋转90º的k倍，当k为1时可省略 左右翻转矩阵 fliplr(A) 对矩阵A实施左右翻转（将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换，…，依次类推。） 上下翻转矩阵 flipud(A) 对矩阵A实施上下翻转 改变矩阵大小 reshape(A,m,n) 在矩阵总元素保持不变的前提下，将矩阵A重新排成m×n的二维矩阵 提取子阵 A(:); A(k,:); A(:,k); A([k l m],[r s]) 取出矩阵A成一列向量 取出矩阵A第k行成一行向量 取出矩阵A第k列成一列向量 取出矩阵A第k、l、m行及第r、s列成一新矩阵

示例4.2 矩阵扩展 (1) 矩阵的数值扩展 A=[1 4;7 8] A(3,3)=5 (2) 合成扩展 C=[1 5;2 4;2 7] d=[A C] (3) 提取子阵 d1=d([1 3],:) d2=d([1 2],[3 5]) (4) 提取矩阵中的某元素 d3=d(2,4)

示例4.3 矩阵重构 已知矩阵 做如下元素操作： ①将矩阵的每一个元素取绝对值构成矩阵Y1； ⑤取矩阵的2、1、3、l列元素构成子矩阵Y5； ⑥矩阵中绝对值大于等于3的元素赋值为0．

A＝[-1 3 -4 5;-3 4 -2 0;4 -2 1 6] Y1=abs(A) Y2=abs(A) Y3=A(abs(A)>2) Y4=A(find([1 0 1 1 0 1 0 1])) Y5=A(:,[2 1 3 1]) Y6=A(abs(A)>=3)=0 数据矩阵的变换操作函数对符号矩阵大部分都适用，可直接应用于符号矩阵。

示例4.4 符号矩阵元素变换操作 A=sym('[sinx(x),cos(x);acos(x),asin(x)]'); B=transpose(A) C=diag(A,1) D=triu(A) E=flipud(A)

4.2.2 数据操作 数据操作名称 操作命令格式 功能简介 最大元素 max(X); [y,k]=max(X); max(A); [Y,U]=max(A); max(A,[],dim); max(A,B); max(A,n) 求向量X的最大元素； y记录X最大元素，k记录序号； 求矩阵A每列的最大元素； Y记录A的每列的最大元素，U记录每列的最大元素的行号； dim=1同max(A)，按列求 dim=2，返回列向量，第i个元素是A上第i行上的最大元素； 两向量或矩阵比较，取对应元素较大者； 同标题n比较，取对应元素较大者 最小元素 min(A) 求矩阵A每列的最小元素 和 sum(X)) sum(A,dim) 求向量X或矩阵A各列元素的和，A可为向量或矩阵，dim用法同上 积 prod(X)) prod(A,dim) 求向量X或矩阵A各列元素的积；dim=1按列求，dim=1按行求 累计和 cumsum(X)) cumsum(A,dim) 求向量X或矩阵A各列元素的累计和；dim=1按列求，dim=1按行求 累计积 cumprod(X)) cumprod(A,dim) 求向量X或矩阵A各列元素的累计积；dim=1按列求，dim=1按行求 平均值 mean(X); mean(A,dim) 求向量X或矩阵A各列元素的平均值；dim=1按列求，dim=1按行求 中值 median(A) median(A,dim) 求向量X或矩阵A各列元素的中值；dim=1按列求，dim=1按行求 标准差 std(X) std(A,flag,dim) 求向量X或矩阵A各列元素的标准差；flag=0，无偏差方差，1为有偏差方差；dim=1按列求，dim=1按行求 相关系数 corrcoef(A) corrcoef(X,Y) 由矩阵A形成的一个相关系数矩阵，把每列作为一个向量； X,Y是向量 元素排序 sort(X) [Y,I]=sort(A,dim) 按升序对向量X中元素进行排序； 按升序对向量X中元素进行排序；Y是排序后的矩阵，I记录Y中的元素在A中位置 行升序排列 sortrows(A) 按升序排列矩阵各行 数值积分 cumtrapz(A,dim) 梯形法求累积数值积分

示例4.6 矩阵数据操作 已知矩阵 做如下计算： (1)求矩阵A每列的最大元素赋给max (2)求矩阵A每列的最小元素赋给min； (3)求矩阵A列元素的平均值赋给me； (4)求矩阵A列元素的中值赋给ded (5)求矩阵A元素的标准差赋给ndl (6)求矩阵A各列元素的和赋给sum

4.3 多项式运算 4.3.1 多项式运算 Matlab有专门的函数进行多项式计算。用赋值 语句将多项式计算结果赋给变量是处理多项式 计算的通常方式。多项式运算中常用命令如表 4.5所示，其中字母a,b表示多项式对应的系数向 量，按未知量次数由高向低排列，空的补零。

多项式运算名称 命令格式 功能简介 计算多项式的值 polyval(a,x) x可以为向量或矩阵 多项式的加法(减法) a+b(a-b) 如果两个多项式的阶数相同，则可直接进行加减，如果两个多项式阶 数不同，则需用旨零填补，使之具有和高阶的多项式一样的阶数。 多项式的乘法 conv(a,b) 矩阵A中每个元素加数k 多项式除法 [q,r]＝deconv(c,a) c(x)多项式除以a(x)多项式 q是商多项式的系数向量,r是余数多项式的系数向量 多项式求根 roots(c) 多项式的系数向量构成的方程的根 多项式的微分 polyder(a) 多项式的导函数 多项式的乘积进行求导 polyder(p1,p2) 求多项式p1与多项式p2乘积的导数 多项式数据拟合 P=polyfit(x,y,n) 给定向量x，y作为对应数据点拟合成n次多项式，P为所求拟合多项式的系数向量，向量x，y具有相同的维数，n为正整数 多项式的估值 po1yval(a,x) 多项式a(x)的，估x分别取值时多项式的值．

在多项式和有理分式的计算过程中使用符号计算 比较简便，常用运算指令见表4.6所示。 4.3.2 多项式符号运算 在多项式和有理分式的计算过程中使用符号计算 比较简便，常用运算指令见表4.6所示。 符号工具箱 命令格式 功能简介 因式分解 factor (s) 符号表达式因式分解函数，s为符号矩阵或符号表达式 符号矩阵的展开 expand (s) 符号表达式展开函数，s为符号矩阵或符号表达式，常用于多项式的因式分解中，也常用于三角函数，指数函数和对数函数的展开中。 同类项合并 collect (s) collect (s, v) s为矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合并 将s中的变量v的同幂项系数合并 符号简化 simple (s) pretty (s) s为矩阵或表达式, 寻找符号矩阵或符号表达式的最简型。 使表达式更加精美 代数式求值 sn=subs(a,'old','new) r=vpa(sn) a是变量替换前的符号表达式的变量名，sn是替换后，old为被替换变量，new为替换变量； r最终求得结果

4.4 矩阵运算 4.4.1 基本数值运算 矩阵运算除包括加、减、乘、除等基本算术运 算外，还包括逆运算、矩阵的行列式、范数、 条件数特征值运算、矩阵秩、迹的运算、求正 交矩阵运算、向量组的最大无关组等复杂的运 算，其功能及其Matlab命令形式见表4.6。

矩阵的函数运算名称 命令格式 功能简介 矩阵加法和减法 A+B（A-B） 将两个同型矩阵相加(减) 数加矩阵 k+A 矩阵A中每个元素加数k 数乘矩阵 k*A 矩阵A中每个元素乘数k 矩阵的乘幂 A^n 计算An 矩阵的乘法 A*B 将两个矩阵进行矩阵相乘 矩阵的左除 A\B 等价于A-1B 矩阵的右除 A/B 等价于AB-1 矩阵的点运算 A.^n(A.*B,A.\B) 矩阵A与矩阵B对应元素作相应运算 求行数、列数 size(A) 测矩阵A是几行几列的矩阵 行列式运算 det(A) 求矩阵A的行列式值 逆运算 inv(A)或A^(-1) 计算矩阵A的逆矩阵 矩阵秩 rank(A) 计算矩阵秩（矩阵线性无关的行数与列数称为矩阵的秩） 矩阵的迹 trace(A) 计算矩阵的迹（矩阵的迹等于矩阵的对角线元素之和，也等于矩阵的特征值之和。） 特征多项式 ploy(A) 求矩阵A的特征多项式 向量的点乘（内积） dot(X,Y) =sum(X.*Y) 维数相同的两个向量的点乘 向量叉乘 cross(X,Y) 两个向量的叉乘是一个过两相交向量的交点且垂直于两向量所在平面的向量 混合积 dot(X,cross(Y,c) 混合积由以上两函数实现 矩阵A的特征值和特征向量 E=eig(A); [V,D]=eig(A); [V,D]=eig(A,‘nobalance’) 求矩阵A的全部特征值，构成向量E； 求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量（先对A作相似变换后求矩阵A的特征值和特征向量）； 直接求矩阵A的特征值和特征向量 矩阵或向量范数 norm(V,1); norm(V)或norm(V,2); norm(V,inf) 计算向量V的1—范数； 计算向量V的2—范数； 计算向量V的∞—范数 矩阵条件数 cond(A,1); cond(A)或cond(A,2); cond(A,inf) 计算A的1—范数下的条件数； 计算A的2—范数数下的条件数； 计算A的 ∞—范数下的条件数 正交化 orth (A) 将矩阵A正交规范化

示例4.7 矩阵运算的应用 某农场饲养的动物所能达到的最大年龄为15岁，将其分为三个年龄组：第一组0~5岁；第二组6~10岁；第三组11~15岁。动物从第二年龄组起开始繁殖后代，经过长期统计，第二年龄组的动物在其年龄段平均繁殖4个后代，第三组在其年龄段平均繁殖3个后代，第一年龄组和第二年龄组的动物能顺利进入下一个年龄组的存活率分别是1/2和1/4。假设农场现有三个年龄段的动物各1000头，问15年后农场饲养的动物总数及农场三个年龄段的动物各将达到多少头?指出15年间动物总增长多少头及总增长率。

即： 其中： 解 年龄组为5岁一段，故将时间周期也取5年。15年经 过3个周期。用k=1，2，3分别表示第一、二、三个周期， 解 年龄组为5岁一段，故将时间周期也取5年。15年经 过3个周期。用k=1，2，3分别表示第一、二、三个周期， xi(k)表示第i个年龄组在第k个周期的数量。由题意，知每 阶段有如下递推关系，方程组表示为： 可化为如下矩阵递推关系： 即 即： 其中： 其中：

Matlab程序设计： x0=[1000 1000 1000]'; L=[0 4 3;1/2 0 0 ;0 1/4 0]; x3=L^3*x0 %15年后三个年龄组的数量 x3 = 14375 1375 875 total=sum(x3) %15年后农场动物总数量 total = 16625 per=x3/total %15年后三个年龄组占总数量的百分数 per = 0.8647 0.0827 0.0526 pie(x3) %绘出三个年龄组的饼形图

结果分析：15年后，农场饲养的动物总数将达 到16625头，其中0~5岁的有14375头，占总数的 86.5%；6~10岁的有1375头，占总数的8.3%； 1~15岁的有875头，占总数的5.2%。

示例4.8 可逆矩阵的应用 密码问题：为了保护您的数据在传递过程中不 被别人窃听或修改，您必须对数据进行加密（加 密后的数据称为密文），这样，即使别人窃取了 您的数据（密文），由于没有密钥而无法将之还 原成明文（未经加密数据），从而保证了数据的 安全性，接收方因有正确的密钥，因此可以将密 文还原成正确的明文。

一种经典的加密方法是在26个英文字母与数字 间建立起一一对应，例如可以是 A B C D E 。。。 X Y Z 1 2 3 4 5 24 25 26 若要发出信息“SEND MONEY”，使用上述代码，则此信息的编码是19，5，14，4，13，15，14，5，25，其中5表示字母E。不幸的是，这种编码很容易被别人破译。在一个较长的信息编码中，人们会根据那个出现频率最高的数值而猜出它代表的是哪个字母，比如上述编码中出现最多次的数值时5，人们自然会想到它代表的是字母E，因为统计规律告诉我们，字母E是英文单词中出现频率最高的。

方法：利用矩阵乘法来对信息的数字编码“明文”SEND MONEY进行加密，让其变成“密文”后再进行传送， 以增加非法用户破译的难度，而让合法用户轻松解密。 取矩阵A作为“密钥”（矩阵A的元素均为整数，而且其 行列式|A|=±1，那么由 可知 A-1的元素也均为整数）。可以利用这样的矩阵A乘以 “明文”所构成矩阵B方法加密生成“密文”矩阵C，使加 密之后的密文很难破译。解密则将“密文”矩阵C左乘 矩阵“密钥”矩阵逆矩阵A-1即可解密得到明文B。

如取 明文“SEND MONEY”对应的9个数值3列被排成以下 矩阵 作矩阵乘 对应着将发出去的密文编码：43，105，81，45， 118，77，49，128，93

合法用户用A-1去左乘上述密文编码构成矩阵C 即可解密得到明文B编码。 为了构造“密钥”矩阵A，我们可以从单位阵I开 始，有限次地使用第三类初等行变换，而且只 用某行的整数倍加到另一行，当然，第一类初 等行变换也能使用。这样得到的矩阵A，其元素 均为整数，而且由于|A|=±1可知，A-1的元素 必然均为整数。

Matlab程序设计： A=[1 2 1;2 5 3;2 3 2 ]; ming=[19 5 14 4 13 15 14 5 25]; B=reshape(ming,3,3); C=A*B; B1=inv(A)*B; mi=reshape(B1,1,9) 矩阵密码法是信息编码与解码的技巧，其中的 一种是基于利用可逆矩阵的方法， 称为Hill（希尔密码）

4.4.3 符号矩阵运算 1、符号矩阵的四则运算同数值矩阵的四则运算 完全相同。 2、符号矩阵的其它一些基本运算也与数值矩阵 的运算格式相同。这些运算包括矩阵的行列式（det）、逆(inv)、秩（rank）、幂（^）和 指数（exp和expm等）运算。 3、符号工具箱中还提供了符号矩阵因式分解、 展开、合并、简化及通分等符号操作函数，它们 同多项式中符号运算完全一样，只是将多项式中 符号运算中变量化为矩阵运算。

4.5 线性方程组的求解 线性方程组是线性代数研究的主要问题，而且 很多实际问题的解决也归结为线性方程组的 求解。在matlab中求解线性方程组主要方法如 表4.8所示。

求解方法名称 求解方法 功能简介 最简行阶梯法 U=rref([A,b]) x=U(:,end) 将其增广矩阵[A,b] 化为最简行阶梯型,等式右端的四个数也就是方程的解 直接法（解低阶稠密矩阵） 方程：AX=b 求逆法 inv(A)*b 如果A为奇异方矩阵，则会给出警告信息 左除 A\b 如果A为方矩阵，则除了算法不同，A\b与inv (A)*b大致相当； 如果A是一个n×n的矩阵，b是一个具有n个元素的列向量或具有多个此类列向量的矩阵，则X=A\b为用Gauss消元法得到的方程AX=b的解; 如果A为m×n的矩阵（m≠n），且b为具有m个元素的列向量或具有多个此类列向量的矩阵，则X=A\b为等式系统AX=b的最小二乘意义上的解。 右除 b/A 当方程形如XA=b时，XAA-1=bA-1 分解法(解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存) LU分解 [L,U]= lu(A) X=U\(L\ b) 高斯消元法求系数阵.A*X = b,A = LU(其中A分解成为下三角阵L和上三角阵U),L*U*X = b,X=U\(L\ b) Cholesky分解 R = chol (A) X = R\ (R’\b) A为对称正定矩阵，A*X = b,A = R’*R(其中A分解成上三角矩阵R和其转置的乘积),R’*R*X = b, X = R\ (R’\b) QR分解 [Q,R]=qr (A) X = R\ (Q\b) A为长方矩阵，A*X = b, A = QR(其中A分解为正交矩阵Q和上三角矩阵R),QRX = b ,X = R\ (Q\b) 线性齐次方程组的通解 z = null (A) z = null(A, ‘r’) z的列向量为方程组的正交规范基，满足于Z’×Z = E z的列向量为方程组A X = 0的有理基 非齐次线性方程组的通解 第一步：判断AX = b是否有解，若有解则进行第二步； 第二步：求AX = b的一个特解； 第三步：求AX = 0的通解； 第四步：AX = b的通解 ( AX = 0的通解+AX = b的一个特解) 符号方程组求解 线性方程组AX=B X =linsolve(A,B) 只给出特解 非线性方程组的解 [x1,x2,x3,…]= Solve(e1,e2,e3,…) 给出非线性方程组的解。其中e1,e2,e3,…是符号方程，x1,x2,x3,…是要求的未知量

示例4.9 线性方程组求解 求一通过(1,6)，(-1,8)，(2,11) 三点的二次方 程式 ，并画出此方程式。 解：设此方程式为 p(t)= a t2 + b t + c 则 p(1)=6，p(-1)=8， p(2)=11； 即解下列线性方程式的解 a +b+ c =6 a-b+c = 8 4a+2b+c = 11

Matlab程序设计： A=[1 1 1;1 -1 1;4 2 1]; b=b=[6 8 11]'; x=A\b x = 2 -1 5 %故p(t)= 2t2 -t +5 为通过(1,6), (-1,8), (2,11)三点的二 次方程式。 v=[2 -1 5]; %建一个代表p(t)系数的向量 t=-4:0.05:4; y=polyval(v,t); plot(t,y)

示例4.10 线性方程组求解 剑桥减肥食谱问题：一种在20世纪80年代很流行的食谱， 称为剑桥食谱，是经过多年研究编制出来的。这是由 Alan H. Howard博士领导的科学家团队经过8年对过度 肥胖病人的临床研究，在剑桥大学完成的。这种低热量的 粉状食品精确地平衡了碳水化合物、高质量的蛋白质和脂 肪、配合维生素、矿物质、微量元素和电解质。为得到所 希望的数量和比例的营养，Howard博士在食谱中加入了 多种食品。每种食品供应了多种所需要的成分，然而没有 按正确的比例。例如脱脂牛奶是蛋白质的主要来源但包含 过多的钙，因此大豆粉用来作为蛋白质的来源，它包含较 少量的钙。然而大豆粉包含过多的脂肪，因而加上乳清， 因乳清含脂肪较少，然而乳清又含有过多的碳水化合物…

在这里把问题简化，看看这个问题小规模的情 形。表是该食谱中的3种食物以及100克每种食 物成分含有某些营养素的数量。 每100客食物所含营养（g） 减肥所要求得每日营养量 脱脂牛奶 大豆面粉 乳清 蛋白质 36 51 13 33 碳水化合物 52 34 74 45 脂肪 7 1.1 3

如果用这三种食物作为每天的主要食物，那么它们的用量应各取多少才能全面准确地实现这个营养要求？ 以100克为一个单位，为了保证减肥所要求的每日营养量，设每日需食用的脱脂牛奶x1个单位，大豆面粉x2个单位，乳清x3个单位，则由所给条件得

Matlab程序设计： A=[36,51,13;52,34,74;0,7,1.1]; b=[33;45;3]; U=rref([A,b])； %命令1:最简行阶梯型 x=U(:,end) 解上方程组得，解为：

即为了保证减肥所要求的每日营养量，每日需 食用脱脂牛奶27.72克，大豆面粉39.19克，乳清 23.32克。 另有格式： inv(A)*b %命令2:求逆法 A\b %命令3:左除

示例4.11 非齐次线性方程组的通解 求解方程组 在Matlab编辑器中建立M文件：LX0601.m 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') %求AX=0的基础解系 else X='equation no solve' %判断无解 end 运行后结果显示： R_A = 2 R_B = 3 X = equition no solve 说明该方程组无解。

示例4.12 非齐次线性方程组的通解 网络流问题 当科学家、工程师或者经济学家 研究一些数量在网络中的流动时自然推导出线 性方程组。例如，城市规划和交通工程人员监 控一个网络状的市区道路的交通流量模式；电 气工程师计算流经电路的电流；以及经济学家 分析通过分销商和零售商的网络从制造商到顾 客的产品销售。许多网络中的方程组涉及成百 甚至上千的变量和方程。

一个网络包含一组称为接合点或节点的点集，并 由称为分支的线或弧连接部分或全部的节点。流 的方向在每个分支上有标示，流量（速度）也有 显示或用变量标记。 网络流的基本假设是全部流入网络的总流量等于 全部流出网络的总流量，且全部流入一个节点的 流量等于全部流出此节点的流量。于是，对于每 个节点的流量可以用一个方程来描述。网络分析 的问题就是确定当局部信息（如网络的输入）已 知时，求每一分支的流量。

交通流问题 图给出了某城市部分单行街道在一个下午早些 时候的交通流量（每小时车辆数目）。计算该 网络的车流量。

由网络流量假设，有 对于节点A： 对于节点B： 对于节点C： 对于节点D： 对于节点E： 于是，所给问题可以归结为如下线性方程组 的求解。

Matlab程序设计： A=[-1,1,0,0,0,0;0,-1,1,-1,1,0;0,0,0,0,-1,1;0,0,0,1,0,-1; 1,0,-1,0,0,0]; b=[50;0;-60;50;-40]; [R,s]=rref([A,b]); [m,n]=size(A); x0=zeros(n,1); r=length(s); x0(s,:)=R(1:r,end); disp('非齐次线性方程组的特解为：') x0 disp('对应齐次线性方程组的基础解系为：') x=null(A,'r')

解这个方程组，得

示例4.13 向量组的线性相关性的应用 通过中成药药方配制问题，理解向量组的线性 相关性、最大线性无关组向量的线性表示以及 向量空间等线性代数的知识。 药方配制问题：某中药厂用9种中草药A-I，根据 不同的比例配制成了7种特效药，各用量成分见 表1（单位：克）。

中药 1号成药 2号成药 3号成药 4号成药 5号成药 6号成药 7号成药 A 10 2 14 12 20 38 100 B 25 35 60 55 C 5 3 11 D 7 9 15 47 E 1 33 6 F 50 G 4 17 39 H 16 I 8

试解答： （1）某医院要购买这7种特效药，但药厂的第3 号药和第6号药已经卖完，请问能否用其他特效 药配制出这两种脱销的药品。 （2）现在该医院想用这7种草药配制三种新的特 效药，表2给出了三种新的特效药的成分，请问 能否配制？如何配制？

中药 1号新药 2号新药 3号新药 A 40 162 88 B 62 141 67 C 14 27 8 D 44 102 51 E 53 60 7 F 50 155 80 G 71 118 38 H 41 68 21 I 52 30

解：（1）把每一种特效药看成一个九维列向量: u1, u2, u3, u4, u5 ,u6, u7 分析7个列向量构成向量组的线性相关性。若向量组线性无关，则无法配制脱销的特效药；若向量组线性相关，且能将u3, u6 用其余向量线性表示，则可以配制3号和6号药品 Matlab代码 u1=[10;12;5;7;0;25;9;6;8]; u2=[2;0;3;9;1;5;4;5;2]; u3=[14;12;11;25;2;35;17;16;12]; u4=[12;25;0;5;25;5;25;10;0]; u5=[20;35;5;15;5;35;2;10;0]; u6=[38;60;14;47;33;55;39;35;6]; u7=[100;55;0;35;6;50;25;10;20]; U=[u1,u2,u3,u4,u5,u6,u7] [U0,r]=rref(U)

计算结果为 U = 10 2 14 12 20 38 100 12 0 12 25 35 60 55 5 3 11 0 5 14 0 7 9 25 5 15 47 35 0 1 2 25 5 33 6 25 5 35 5 35 55 50 9 4 17 25 2 39 25 6 5 16 10 10 35 10 8 2 12 0 0 6 20 U0 = 1 0 1 0 0 0 0 0 1 2 0 0 3 0 0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 r = 1 2 4 5 7 故可以配制3号和6号药。

（2）三种新药用v1，v2，v3表示，问题化为v1， v2，v3能否由u1-u7线性表示，若能表示，则可配 制；否则，不能配制。 令v1=[40 62 14 44 53 50 71 41 14]'; v2=[162 141 27 102 60 155 118 68 52]'; v3=[88 67 8 51 7 80 38 21 30]'; U=[u1,u2,u3,u4,u5,u6,u7,v1,v2,v3] [U0,r]=rref(U)

计算结果为 v1 v2 v3 U0 = 1 0 1 0 0 0 0 1 3 0 0 1 2 0 0 3 0 3 4 0 0 0 0 1 0 1 0 2 2 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 r = 1 2 4 5 7 10

由r=1,2,4,5,7,10可以看出一个最大无关组为： u1, u2, u4, u5, u7,v3, 且 v1=u1+3u2+2u4， v2=3u1+4u2+2u4+u7 由于v3在最大无关组,不能被线性表示，所以无法配制。

4.6 矩阵特征值和特征向量 特征值与特征向量是线性代数中非常重要的 概念，在实际的工程应用中占有非常重要的 地位。在本节中要介绍如何利用 去求特征值 与特征向量、矩阵的对角化等问题，培养把实 际问题转化为数学问题来求解的能力。

矩阵的函数运算名称 命令格式 功能简介 矩阵A的特征值和特征向量 E=eig(A); [V,D]=eig(A); [V,D]=eig(A,‘nobalance’) 求矩阵A的全部特征值，构成向量E； 求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量（先对A作相似变换后求矩阵A的特征值和特征向量）； 直接求矩阵A的特征值和特征向量 正交化 [U, T]= schur (A) T = schur (A) U为正交矩阵，使得A = UTU’和U’U=E 生成schur矩阵T。当A为实对称矩阵时，T为特征值对角阵。

示例4.13 矩阵的schur分解 A=[4 0 0;0 3 1;0 1 3]; [U,T]=schur(A) U = 0 0 1.0000 -0.7071 0.7071 0 0.7071 0.7071 0 T = 2 0 0 0 4 0 0 0 4 这里U就是所求的正交矩阵P，T就是对角矩阵Λ。

[V,D]=eig(A) V = 0 0 1.0000 -0.7071 0.7071 0 0.7071 0.7071 0 D = 2 0 0 0 4 0 0 0 4 这里V就是所求的正交矩阵P，D就是对角矩阵Λ。 说明：对于实对称矩阵，用eig和schur分解效果一样。

示例4.14 矩阵对角化的应用 行业就业人数预测：设某中小城市及郊区乡镇共 有30万人从事农、工、商工作，假定这个总人数 在若干年内保持不变，而社会调查表明： （1）在这30万就业人员中，目前约有15万人从 事农业，9万人从事工业，6万人经商。 （2）在务农人员中，每年约有20%改为务工， 10%改为经商。 （3）在务工人员中，每年约有20%改为务农， （4）在经商人员中，每年约有10%改为务农， 10%改为务工。 现欲预测一、二年后从事各业人员的人数，以 及经过多年之后，从事各业人员总数之发展趋势。

解：若用3维向量[xi,yi,zi]‘表示第i年后从事这三种 职业的人员总数，则已知[x0,y0,z0]‘=[5,9,6]‘， 问题是求第一年[x1,y1,z1]'、第二年[x2,y2,z2]‘ 这三种职业的人员总数，以及在n→∞时 [xn,yn,zn]'的发展趋势。 依题意，一年后，从事农、工、商的人员总数 应为

即 进而推得 即n年之后从事各业人员的人数完全由An决定。 事实上，运用实对称矩阵的正交对角化方法， 可以轻松求得An。

Matlab代码 A=[0.7 0.2 0.1;0.2 0.7 0.1;0.1 0.1 0.8]; t0=[5 9 6]'; t1=A*t0; t2=A^2*t0 无限增加时间n，三种职业的人员总数的发展趋 势过程趋向于一个稳态值。

求A的特征值和特征向量，得到 [V,D]=eig(A) V = 0.7071 0.4082 0.5774 -0.7071 0.4082 0.5774 0.0000 -0.8165 0.5774 D = 0.5000 0 0 0 0.7000 0 0 0 1.0000 其中：V是由特征向量组成的矩阵；D是由全部特征 值构成的对角阵

三个特征向量线性无关，由相似对角的理论知， 矩阵A可对角化，即 则 由于A的最大特征值为1，所以 收敛， 因此

tninf=P*Tinf*inv(P)*t0 tninf = 20/3 Matlab代码 syms n Tn=T^n Tn = [ (1/2)^n, 0, 0] [ 0, (7/10)^n, 0] [ 0, 0, 1] Tinf=limit(Tn,n,inf) Tinf = [ 0, 0, 0] [ 0, 0, 1] tninf=P*Tinf*inv(P)*t0 tninf = 20/3 double(tninf) ans = 6.6667

结论：三种职业的人员总数的发展趋势过程趋 向于一个稳态值，其比例为1：1：1，都为6.6667 万人。 画出近20年的人口变化图形 Matlab代码 for x=1:20 y=A^x*t0; plot(x,y,'*') hold on end 由图中可以看出，从第15后，三种职业的人员总 数的发展趋势过程趋向于一个稳态值。

总结：马尔科夫分析模型 实际分析中，往往需要知道经过一段时间后，市 场趋势分析对象可能处于的状态，这就要求建立 一个能反映变化规律的数学模型。马尔科夫市场 趋势分析模型是利用概率建立一种随机型的时序 模型，并用于进行市场趋势分析的方法。 马尔科夫分析法的基本模型为： X(k+1)=X(k)×P 公式中：X(k)表示趋势分析与预测对象在t=k时 刻的状态向量，P表示一步转移矩阵概率， X(k+1)表示趋势分析与预测对象在t=k+1时刻的 状态向量。

必须指出的是，上述模型只适用于具有马尔科 夫性的时间序列，并且各时刻的状态转移概率 保持稳定。若时间序列的状态转移概率随不同 的时刻在变化，不宜用此方法。由于实际的客 观事物很难长期保持同一状态的转移概率，故 此法一般适用于短期的趋势分析与预测。