第五讲 线性代数中的数值计算问题
【引 例 】求下列三阶线性代数方程组的近似解 MATLAB程序为: A=[2 -5 4;1 5 -2;-1 2 4]; b=[5;6;5]; x=A\b
在MATLAB命令窗口,先输入下列命令构造系数矩阵A和右端向量b: 1 5 -2 -1 2 4 b=[5;6;5] b = 5 6 5 然后只需输入命令x=A\b即可求得解x: x=A\b x = 2.7674 1.1860 1.3488
一、 特殊矩阵的实现
1.零矩阵:所有元素值为零的矩阵称为零矩阵。零矩阵可 以用zeros函数实现。zeros是MATLAB内部函数,使用格式 如下: 一、 特殊矩阵的实现 常见的特殊矩阵有零矩阵、幺矩阵、单位矩阵、三角形矩阵等,这类特殊矩阵在线性代数中具有通用性;还有一类特殊矩阵在专门学科中有用,如有名的希尔伯特(Hilbert)矩阵、范德蒙(Vandermonde) 矩阵等。 1.零矩阵:所有元素值为零的矩阵称为零矩阵。零矩阵可 以用zeros函数实现。zeros是MATLAB内部函数,使用格式 如下: zeros(m):产生m m阶零矩阵; zeros(m,n):产生m n阶零矩阵,当m=n时等同于zeros(m); zeros(size(A)):产生与矩阵A同样大小的零矩阵。
2.幺矩阵:所有元素值为1的矩阵称为幺矩阵。幺矩阵可以用ones函数实现。它的调用格式与zeros函数一样。 一、 特殊矩阵的实现 2.幺矩阵:所有元素值为1的矩阵称为幺矩阵。幺矩阵可以用ones函数实现。它的调用格式与zeros函数一样。 【例1】 试用ones分别建立32阶幺矩阵、和与前例矩阵A同样大小的幺矩阵。 用ones(3,2) 建立一个3 2阶幺阵: ones(3,2) % 一个32阶幺阵 ans =1 1 1 1
3.单位矩阵:主对角线的元素值为1、其余元素值为0的 矩阵称为单位矩阵。它可以用MATLAB内部函数eye建立, 使用格式与zeros相同。 一、 特殊矩阵的实现 3.单位矩阵:主对角线的元素值为1、其余元素值为0的 矩阵称为单位矩阵。它可以用MATLAB内部函数eye建立, 使用格式与zeros相同。 4.数量矩阵:主对角线的元素值为一常数d、其余元素值为0的矩阵称为数量矩阵。显然,当d=1时,即为单位矩阵,故数量矩阵可以用eye(m)*d或eye(m,n)*d建立。 5.对角阵:对角线的元素值为常数、其余元素值为0的 矩阵称为对角阵。我们可以通过MATLAB内部函数diag, 利用一个向量构成对角阵;或从矩阵中提取某对角线构 成一个向量。使用 格式为diag(V)和diag(V,k)两种。
【例2】已知向量v,试建立以向量v作为主对角线的对角阵A;建立分别以向量v作为主对角线两侧的对角线的对角阵B和C。 一、 特殊矩阵的实现 6.用一个向量V构成一个对角阵 设V为具有m个元素的向量,diag(V)将产生一个mm阶对角阵,其主对角线的元素值即为向量的元素值;diag(V,k)将产生一个nn(n=m+|k|,k为一整数)阶对角阵,其第k条对角线的元素值即为向量的元素值。注意:当k>0,则该对角线位于主对角线的上方第k条;当k<0,该对角线位于主对角线的下方第|k|条;当k=0,则等同于diag(V)。用diag建立的对角阵必是方阵。 【例2】已知向量v,试建立以向量v作为主对角线的对角阵A;建立分别以向量v作为主对角线两侧的对角线的对角阵B和C。 MATLAB程序如下:
一、 特殊矩阵的实现 v =[1;2;3]; % 建立一个已知的向量A A=diag(v) A= 1 0 0 0 2 0 0 0 3 % 按各种对角线情况构成相应的对角阵A、B和C v =[1;2;3]; % 建立一个已知的向量A A=diag(v) A= 1 0 0 0 2 0 0 0 3 B=diag(v,1) B = 0 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 C=diag(v,-1) C = 0 0 0 0 1 0 0 0 0 2 0 0
一、 特殊矩阵的实现 7.从矩阵中提取某对角线 我们也可以用diag从矩阵中提取某对角线构成一个向量。设A为m n阶矩阵,diag(A)将从矩阵A中提取其主对角线产生一个具有min(m,n)个元素的向量。diag(A,k)的功能是: 当k>0,则将从矩阵A中提取位于主对角线的上方第k条对角线构成一个具有n-k个元素的向量;当k<0,则将从矩阵A中提取位于主对角线的下方第|k|条对角线构成一个具有m+k个元素的向量;当k=0,则等同于diag(A)。
【例3】 已知矩阵A,试从矩阵A分别提取主对角线及它两侧的对角线构成向量B、C和D。 MATLAB程序如下: 一、 特殊矩阵的实现 【例3】 已知矩阵A,试从矩阵A分别提取主对角线及它两侧的对角线构成向量B、C和D。 MATLAB程序如下: A=[1 2 3;4 5 6]; % 建立一个已知的23阶矩阵A % 按各种对角线情况构成向量B、C和D B=diag(A) B = 1 5 C=diag(A,1) C = 2 6 D=diag(A,-1) D = 4
8.上三角阵:使用格式为triu(A)、triu(A,k) 一、 特殊矩阵的实现 8.上三角阵:使用格式为triu(A)、triu(A,k) 设A为mn阶矩阵,triu(A)将从矩阵A中提取主对角线之上的上三角部分构成一个m n阶上三角阵;triu(A,k)将从矩阵A中提取主对角线第|k|条对角线之上的上三角部分构成一个mn阶上三角阵。注意:这里的k与diag(A,k)的用法类似,当k>0,则该对角线位于主对角线的上方第k条;当k<0,该对角线位于主对角线的下方第|k|条;当k=0,则等同于triu (A)
一、 特殊矩阵的实现 【例4】试分别用triu(A)、triu(A,1)和、triu(A,-1)从矩阵A提取相应的上三角部分构成上三角阵B、C和D。 MATLAB程序如下: A=[1 2 3;4 5 6;7 8 9;9 8 7]; % 一个已知的43阶矩阵A % 构成各种情况的上三角阵B、C和D B=triu(A) B = 1 2 3 0 5 6 0 0 9 0 0 0 C=triu(A,1) D=triu(A,-1) 9.下三角阵:使用格式为tril(A)、tril(A,k) tril的功能是从矩阵A中提取下三角部分构成下三角阵。用法与triu相同。
10.空矩阵 A=[0.1 0.2 0.3;0.4 0.5 0.6]; B=find(A>1.0) B = [ ] B=[ ] 一、 特殊矩阵的实现 10.空矩阵 在MATLAB里,把行数、列数为零的矩阵定义为空矩阵。空矩阵在数学意义上讲是空的,但在MATLAB里确是很有用的。例如 A=[0.1 0.2 0.3;0.4 0.5 0.6]; B=find(A>1.0) B = [ ] 这里[ ]是空矩阵的符号,B=find(A>1.0)表示列出矩阵A中值大于1.0的元素的序号。当不能满足括号中的条件时,返回空矩阵。另外,也可以将空矩阵赋给一个变量,如: B=[ ]
二、矩阵的特征值 与特征向量
对于NN阶方阵A,所谓A的特征值问题是:求数λ和N维非零向量x(通常为复数),使之满足下式: 二、矩阵的特征值与特征向量 对于NN阶方阵A,所谓A的特征值问题是:求数λ和N维非零向量x(通常为复数),使之满足下式: A. x=λ x 则称λ为矩阵A的一个特征值(特征根),而非零向量x为矩阵A的特征值λ所对应的特征向量。 对一般的N N阶方阵A,其特征值通常为复数,若A为实对称矩阵,则A的特征值为实数。
二、矩阵的特征值与特征向量 MATLAB提供的内部函数eig可以用来计算特征值与特征向量。eig函数的使用格式有五种,其中常见的有E=eig(A)、[V,D]=eig(A)和[V,D]=eig(A,’nobalance’)三种,另外两种格式用来计算矩阵的广义特征值与特征向量:E=eig(A,B)和[V,D]=eig(A,B)。
(1) E=eig(A):由eig(A)返回方阵A的N个特征值,构成向量E; 二、矩阵的特征值与特征向量 (1) E=eig(A):由eig(A)返回方阵A的N个特征值,构成向量E; (2) [V,D]=eig(A):由eig(A)返回方阵A的N个特征值,构成NN阶对角阵D,其对角线上的N个元素即为相应的特征值,同时将返回相应的特征向量赋予NN阶方阵V的对应列,且A、V、D满足AV=V D; (3) [V,D]=eig(A,’nobalance’):本格式的功能与格式(2)一样,只是格式(2)是先对A作相似变换(balance),然后再求其特征值与相应的特征向量;而本格式则事先不作相似变换;
(4) E=eig(A,B):由eig(A,B)返回NN阶方阵A和B的N个广义特征值,构成向量E。 二、矩阵的特征值与特征向量 (4) E=eig(A,B):由eig(A,B)返回NN阶方阵A和B的N个广义特征值,构成向量E。 (5) [V,D]=eig(A,B):由eig(A,B)返回方阵A和B的N个广义特征值,构成N N阶对角阵D,其对角线上的N个元素即为相应的广义特征值,同时将返回相应的特征向量构成NN阶满秩矩阵,且 满足AV=B V D。
【例5】试用格式(1)求下列对称矩阵A的特征值;用格式(2)求A的特征值和相应的特征向量,且验证之。 二、矩阵的特征值与特征向量 【例5】试用格式(1)求下列对称矩阵A的特征值;用格式(2)求A的特征值和相应的特征向量,且验证之。 A =[ 1.0000 1.0000 0.5000 1.0000 1.0000 0.2500 0.5000 0.2500 2.0000 ]; 执行eig(A)将直接获得对称矩阵A的三个实特征值:
而下列命令则将其三个实特征值作为向量赋予变量E: E=eig(A) E = -0.0166 二、矩阵的特征值与特征向量 eig(A) ans = -0.0166 1.4801 2.5365 而下列命令则将其三个实特征值作为向量赋予变量E: E=eig(A) E = -0.0166
三、行列式的值
【例6】利用随机函数产生一个三阶方阵A,然后计算方阵之行列式的值。 三、行列式的值 MATLAB提供的内部函数det用来计算矩阵的行列式的值。设矩阵A为一方阵(必须是方阵),求矩阵A的行列式值的格式为:det(A)。注意:本函数同样能计算通过构造出的稀疏矩阵的行列式的值。关于如何构造稀疏矩阵,将在本章最后一节介绍。 【例6】利用随机函数产生一个三阶方阵A,然后计算方阵之行列式的值。 A=rand(3) A = 0.9501 0.4860 0.4565 0.2311 0.8913 0.0185 0.6068 0.7621 0.8214 det(A) ans = 0.4289
四、 矩阵求逆及其 线性代数方程组求解
则称A为B的逆矩阵,或称B为A的逆矩阵。这时A,B都称为可逆矩阵(或非奇异矩阵、或满秩矩阵),否则称为不可逆矩阵(或奇异矩阵、或降秩矩阵)。 四、矩阵求逆及其线性代数方程组求解 1 . 矩阵求逆 若方阵A,B满足等式 A*B = B*A = I (I为单位矩阵) 则称A为B的逆矩阵,或称B为A的逆矩阵。这时A,B都称为可逆矩阵(或非奇异矩阵、或满秩矩阵),否则称为不可逆矩阵(或奇异矩阵、或降秩矩阵)。
【例7】试用inv函数求方阵A的逆阵A-1赋值给B,且验证A与A-1是互逆的。 四、矩阵求逆及其线性代数方程组求解 【例7】试用inv函数求方阵A的逆阵A-1赋值给B,且验证A与A-1是互逆的。 A=[1 -1 1;5 -4 3;2 1 1]; B=inv(A) B = -1.4000 0.4000 0.2000 0.2000 -0.2000 0.4000 2.6000 -0.6000 0.2000 A*B ans = 1.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 1.0000 B*A ans = 1.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 1.0000
利用求系数矩阵A的逆阵A-1,我们可以得到矩阵求逆解法。对于线性代数方程组Ax=b,等号两侧各左乘A-1,有: 四、矩阵求逆及其线性代数方程组求解 2. 矩阵求逆解法 利用求系数矩阵A的逆阵A-1,我们可以得到矩阵求逆解法。对于线性代数方程组Ax=b,等号两侧各左乘A-1,有: A-1Ax=A-1b 由于A-1A=I,故得: x=A-1b
【例8】试用矩阵求逆解法求解例6.20中矩阵A为系数矩阵的线性代数方程组Ax=b的解。 A=[1 -1 1;5 -4 3;2 1 1]; 四、矩阵求逆及其线性代数方程组求解 【例8】试用矩阵求逆解法求解例6.20中矩阵A为系数矩阵的线性代数方程组Ax=b的解。 A=[1 -1 1;5 -4 3;2 1 1]; b=[2;-3;1]; x=inv(A)*b x = -3.8000 1.4000 7.2000
对于线性代数方程组Ax=b,我们可以运用左除运算符“\”象解一元一次方程那样简单地求解: x=A\b 四、矩阵求逆及其线性代数方程组求解 3. 直接解法 对于线性代数方程组Ax=b,我们可以运用左除运算符“\”象解一元一次方程那样简单地求解: x=A\b 当系数矩阵A为N*N的方阵时,MATLAB会自行用高斯消去法求解线性代数方程组。若右端项b为N*1的列向量,则x=A\b可获得方程组的数值解x(N*1的列向量);若右端项b为N*M的矩阵,则x=A\b可同时获得同一系数矩阵A、M个方程组数值解x(为N*M的矩阵),即x(:,j)=A\b(:,j),j=1,2,…M。
四、矩阵求逆及其线性代数方程组求解
解法1:分别解方程组 (1)Ax=b1;(2)Ay=b2 四、矩阵求逆及其线性代数方程组求解 解法1:分别解方程组 (1)Ax=b1;(2)Ay=b2 A=[1 -1 1;5 -4 3;2 1 1]; b1=[2;-3;1]; b2=[3;4;-5]; x=A\b1 x = -3.8000 1.4000 7.2000 y=A\b2 -3.6000 -2.2000 4.4000 得两个线性代数方程组的解: (1) x1= -3.8, x2= 1.4, x3= 7.2; (2) y1= -3.8, y2= 1.4, y3= 7.2
解法2:将两个方程组连在一起求解:Az=b 四、矩阵求逆及其线性代数方程组求解 解法2:将两个方程组连在一起求解:Az=b b=[2 3;-3 4;1 -5] z=A\b z = -3.8000 -3.6000 1.4000 -2.2000 7.2000 4.4000 很明显,这里的解z的两个列向量便是前面分别求得的两组解x和y
习题 1、解方程组Ax=b,分别用求逆解法与直接解法求其解。 2、编一个m程序,求N阶方阵A的行列式的值。