matlab相关矩阵操作及相关实例

独善其身 · 发表于 2016-4-7 16:21

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

　　矩阵及其基本运算
　　1.1 矩阵的表示
　　1.1.1 数值矩阵的生成
　　1.实数值矩阵输入
　　MATLAB的强大功能之一体现在能直接处理向量或矩阵。当然首要任务是输入待处理的向量或矩阵。
　　不管是任何矩阵(向量)，我们可以直接按行方式输入每个元素：同一行中的元素用逗号(，)或者用空格符来分隔，且空格个数不限;不同的行用分号(;)分隔。所有元素处于一方括号([ ])内;当矩阵是多维(三维以上)，且方括号内的元素是维数较低的矩阵时，会有多重的方括号。如：

　　>> Time = [11 12 1 2 3 4 5 6 7 8 9 10]
　　Time =
　　11 12 1 2 3 4 5 6 7 8 9 10
　　>> X_Data = [2.32 3.43;4.37 5.98]
　　X_Data =
　　2.43 3.43
　　4.37 5.98
　　>> vect_a = [1 2 3 4 5]
　　vect_a =
　　1 2 3 4 5
　　>> Matrix_B = [1 2 3;
　　>> 2 3 4;3 4 5]
　　Matrix_B = 1 2 3
　　2 3 4
　　3 4 5
　　>> Null_M = [ ] %生成一个空矩阵

复制代码

　　2.复数矩阵输入
　　复数矩阵有两种生成方式：
　　第一种方式

　　例1-1

>> a="2".7;b=13/25;
　　>> C=[1,2*a+i*b,b*sqrt(a); sin(pi/4),a+5*b,3.5+1]
　　C=
　　1.0000 5.4000 + 0.5200i 0.8544
　　0.7071 5.3000 4.5000
　　第2种方式
　　例1-2
　　>> R=[1 2 3;4 5 6], M=[11 12 13;14 15 16]
　　R =
　　1 2 3
　　4 5 6
　　M =
　　11 12 13
　　14 15 16
　　>> CN="R"+i*M
　　CN =
　　1.0000 +11.0000i 2.0000 +12.0000i 3.0000 +13.0000i
　　4.0000 +14.0000i 5.0000 +15.0000i 6.0000 +16.0000i

复制代码

　　1.1.2 符号矩阵的生成
　　在MATLAB中输入符号向量或者矩阵的方法和输入数值类型的向量或者矩阵在形式上很相像，只不过要用到符号矩阵定义函数sym，或者是用到符号定义函数syms，先定义一些必要的符号变量，再像定义普通矩阵一样输入符号矩阵。
　　1.用命令sym定义矩阵：
　　这时的函数sym实际是在定义一个符号表达式，这时的符号矩阵中的元素可以是任何的符号或者是表达式，而且长度没有限制，只是将方括号置于用于创建符号表达式的单引号中。如下例：
　　例1-3

>> sym_matrix = sym('[a b c;Jack，Help Me!，NO WAY!]，')
　　sym_matrix =
　　[a b c]
　　[Jack Help Me! NO WAY!]
　　>> sym_digits = sym('[1 2 3;a b c;sin(x)cos(y)tan(z)]')
　　sym_digits =
　　[1 2 3]
　　[a b c]
　　[sin(x)cos(y)tan(z)]

复制代码

　　2.用命令syms定义矩阵
　　先定义矩阵中的每一个元素为一个符号变量，而后像普通矩阵一样输入符号矩阵。
　　例1-4

　　>> syms a b c ;
　　>> M1 = sym('Classical');
　　>> M2 = sym(' Jazz');
　　>> M3 = sym('Blues')
　　>> syms_matrix = [a b c; M1， M2， M3;int2str([2 3 5])]
　　syms_matrix =
　　[ a b c]
　　[Classical Jazz Blues]
　　[ 2 3 5]

复制代码

　　把数值矩阵转化成相应的符号矩阵。
　　数值型和符号型在MATLAB中是不相同的，它们之间不能直接进行转化。MATLAB提供了一个将数值型转化成符号型的命令，即sym。
　　例1-5
　　>> Digit_Matrix = [1/3 sqrt(2) 3.4234;exp(0.23) log(29) 23^(-11.23)]
　　>> Syms_Matrix = sym(Digit_Matrix)
　　结果是：
　　Digit_Matrix =
　　0.3333 1.4142 3.4234
　　1.2586 3.3673 0.0000
　　Syms_Matrix =
　　[ 1/3， sqrt(2)， 17117/5000]
　　[5668230535726899*2^(-52)，7582476122586655*2^(-51)，5174709270083729*2^(-103)]
　　注意：矩阵是用分数形式还是浮点形式表示的，将矩阵转化成符号矩阵后，都将以最接近原值的有理数形式表示或者是函数形式表示。
　　1.1.3 大矩阵的生成
　　对于大型矩阵，一般创建M文件，以便于修改：
　　例1-6 用M文件创建大矩阵，文件名为example.m
　　exm=[ 456 468 873 2 579 55
　　21 687 54 488 8 13
　　65 4567 88 98 21 5
　　456 68 4589 654 5 987
　　5488 10 9 6 33 77]
　　在MATLAB窗口输入：
　　>>example;
　　>>size(exm) %显示exm的大小
　　ans=
　　5 6 %表示exm有5行6列。
　　1.1.4 多维数组的创建
　　函数 cat
　　格式 A="cat"(n,A1,A2,…,Am)
　　说明 n="1和n"=2时分别构造[A1;A2]和[A1，A2]，都是二维数组，而n=3时可以构造出三维数组。
　　例1-7

　　>> A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A2;
　　>> A4=cat(3,A1,A2,A3)
　　A4(:,:,1) =
　　1 2 3
　　4 5 6
　　7 8 9
　　A4(:,:,2) =
　　1 4 7
　　2 5 8
　　3 6 9
　　A4(:,:,3) =
　　0 -2 -4
　　2 0 -2
　　4 2 0

复制代码

　　或用另一种原始方式可以定义：
　　例1-8

　　>> A1=[1,2,3;4,5,6;7,8,9];A2=A1';A3=A1-A2;
　　>> A5(:,:,1)=A1, A5(:,:,2)=A2, A5(:,:,3)=A3
　　A5(:,:,1) =
　　1 2 3
　　4 5 6
　　7 8 9
　　A5(:,:,2) =
　　1 4 7
　　2 5 8
　　3 6 9
　　A5(:,:,3) =
　　0 -2 -4
　　2 0 -2
　　4 2 0

复制代码

　　1.1.5 特殊矩阵的生成
　　命令全零阵
　　函数 zeros
　　格式 B = zeros(n) %生成n×n全零阵
　　B = zeros(m,n) %生成m×n全零阵
　　B = zeros([m n]) %生成m×n全零阵
　　B = zeros(d1,d2,d3…) %生成d1×d2×d3×…全零阵或数组
　　B = zeros([d1 d2 d3…]) %生成d1×d2×d3×…全零阵或数组
　　B = zeros(size(A)) %生成与矩阵A相同大小的全零阵
　　命令单位阵
　　函数 eye
　　格式 Y = eye(n) %生成n×n单位阵
　　Y = eye(m,n) %生成m×n单位阵
　　Y = eye(size(A)) %生成与矩阵A相同大小的单位阵
　　命令全1阵
　　函数 ones
　　格式 Y = ones(n) %生成n×n全1阵
　　Y = ones(m,n) %生成m×n全1阵
　　Y = ones([m n]) %生成m×n全1阵
　　Y = ones(d1,d2,d3…) %生成d1×d2×d3×…全1阵或数组
　　Y = ones([d1 d2 d3…]) %生成d1×d2×d3×…全1阵或数组
　　Y = ones(size(A)) %生成与矩阵A相同大小的全1阵
　　命令均匀分布随机矩阵
　　函数 rand
　　格式 Y = rand(n) %生成n×n随机矩阵，其元素在(0，1)内
　　Y = rand(m,n) %生成m×n随机矩阵
　　Y = rand([m n]) %生成m×n随机矩阵
　　Y = rand(m,n,p,…) %生成m×n×p×…随机矩阵或数组
　　Y = rand([m n p…]) %生成m×n×p×…随机矩阵或数组
　　Y = rand(size(A)) %生成与矩阵A相同大小的随机矩阵
　　rand %无变量输入时只产生一个随机数
　　s = rand('state') %产生包括均匀发生器当前状态的35个元素的向量
　　rand('state', s) %使状态重置为s
　　rand('state', 0) %重置发生器到初始状态
　　rand('state', j) %对整数j重置发生器到第j个状态
　　rand('state', sum (100*clock)) %每次重置到不同状态
　　例1-9 产生一个3×4随机矩阵
　　>> R="rand"(3,4)
　　R =
　　0.9501 0.4860 0.4565 0.4447
　　0.2311 0.8913 0.0185 0.6154
　　0.6068 0.7621 0.8214 0.7919
　　例1-10 产生一个在区间[10, 20]内均匀分布的4阶随机矩阵
　　>> a="10";b=20;
　　>> x="a"+(b-a)*rand(4)
　　x =
　　19.2181 19.3547 10.5789 11.3889
　　17.3821 19.1690 13.5287 12.0277
　　11.7627 14.1027 18.1317 11.9872
　　14.0571 18.9365 10.0986 16.0379
　　命令正态分布随机矩阵
　　函数 randn
　　格式 Y = randn(n) %生成n×n正态分布随机矩阵
　　Y = randn(m,n) %生成m×n正态分布随机矩阵
　　Y = randn([m n]) %生成m×n正态分布随机矩阵
　　Y = randn(m,n,p,…) %生成m×n×p×…正态分布随机矩阵或数组
　　Y = randn([m n p…]) %生成m×n×p×…正态分布随机矩阵或数组
　　Y = randn(size(A)) %生成与矩阵A相同大小的正态分布随机矩阵
　　randn %无变量输入时只产生一个正态分布随机数
　　s = randn('state') %产生包括正态发生器当前状态的2个元素的向量
　　s = randn('state', s) %重置状态为s
　　s = randn('state', 0) %重置发生器为初始状态
　　s = randn('state', j) %对于整数j重置状态到第j状态
　　s = randn('state', sum(100*clock)) %每次重置到不同状态
　　例1-11 产生均值为0.6，方差为0.1的4阶矩阵
　　>> mu="0".6; sigma="0".1;
　　>> x="mu"+sqrt(sigma)*randn(4)
　　x =
　　0.8311 0.7799 0.1335 1.0565
　　0.7827 0.5192 0.5260 0.4890
　　0.6127 0.4806 0.6375 0.7971
　　0.8141 0.5064 0.6996 0.8527
　　命令产生随机排列
　　函数 randperm
　　格式 p = randperm(n) %产生1～n之间整数的随机排列
　　例1-12
　　>> randperm(6)
　　ans =
　　3 2 1 5 4 6
　　命令产生线性等分向量
　　函数 linspace
　　格式 y = linspace(a,b) %在(a, b)上产生100个线性等分点
　　y = linspace(a,b,n) %在(a, b)上产生n个线性等分点
　　命令产生对数等分向量
　　函数 logspace
　　格式 y = logspace(a,b) %在( )之间产生50个对数等分向量
　　y = logspace(a,b,n)
　　y = logspace(a,pi)
　　命令计算矩阵中元素个数
　　n = numel(a) %返回矩阵A的元素的个数
　　命令产生以输入元素为对角线元素的矩阵
　　函数 blkdiag
　　格式 out = blkdiag(a,b,c,d,…) %产生以a,b,c,d,…为对角线元素的矩阵
　　例1-13
　　>> out = blkdiag(1,2,3,4)
　　out =
　　1 0 0 0
　　0 2 0 0
　　0 0 3 0
　　0 0 0 4
　　命令友矩阵
　　函数 compan
　　格式 A = compan(u) %u为多项式系统向量，A为友矩阵，A的第1行元素为 -u (2:n)/u(1)，其中u (2:n)为u的第2到第n个元素，A为特征值就是多项式的特征根。
　　例1-14 求多项式的友矩阵和根
　　>> u=[1 0 -7 6];
　　>> A="compan"(u) %求多项式的友矩阵
　　A =
　　0 7 -6
　　1 0 0
　　0 1 0
　　>> eig(A) %A的特征值就是多项式的根
　　ans =
　　-3.0000
　　2.0000
　　1.0000
　　命令 hadamard矩阵
　　函数 hadamard
　　格式 H = hadamard(n) %返回n阶hadamard矩阵
　　例1-15
　　>> h="hadamard"(4)
　　h =
　　1 1 1 1
　　1 -1 1 -1
　　1 1 -1 -1
　　1 -1 -1 1
　　命令 Hankel方阵
　　函数 hankel
　　格式 H = hankel(c) %第1列元素为c，反三角以下元素为0。
　　H = hankel(c,r) %第1列元素为c，最后一行元素为r，如果c的最后一个元素与r的第一个元素不同，交叉位置元素取为c的最后一个元素。
　　例1-16
　　>> c="1:3",r=7:10
　　c =
　　1 2 3
　　r =
　　7 8 9 10
　　>> h="hankel"(c,r)
　　h =
　　1 2 3 8
　　2 3 8 9
　　3 8 9 10
　　命令 Hilbert矩阵
　　函数 hilb
　　格式 H = hilb(n) %返回n阶Hilbert矩阵，其元素为H(i,j)=1/(i+j-1)。
　　例1-17 产生一个3阶Hilbert矩阵
　　>> format rat %以有理形式输出
　　>> H="hilb"(3)
　　H =
　　1 1/2 1/3
　　1/2 1/3 1/4
　　1/3 1/4 1/5
　　命令逆Hilbert矩阵
　　函数 invhilb
　　格式 H = invhilb(n) %产生n阶逆Hilbert矩阵
　　命令 Magic(魔方)矩阵
　　函数 magic
　　格式 M = magic(n) %产生n 阶魔方矩阵
　　例1-18
　　>> M="magic"(3)
　　M =
　　8 1 6
　　3 5 7
　　4 9 2
　　命令 Pascal矩阵
　　函数 pascal
　　格式 A = pascal(n) %产生n阶Pascal矩阵，它是对称、正定矩阵，它的元素由Pascal三角组成，它的逆矩阵的所有元素都是整数。
　　A = pascal(n,1) %返回由下三角的Cholesky系数组成的Pascal矩阵
　　A = pascal(n,2) %返回Pascal(n,1)的转置和交换的形式
　　例1-19
　　>> A="pascal"(4)
　　A =
　　1 1 1 1
　　1 2 3 4
　　1 3 6 10
　　1 4 10 20
　　>> A="pascal"(3,1)
　　A =
　　1 0 0
　　1 -1 0
　　1 -2 1
　　>> A="pascal"(3,2)
　　A =
　　1 1 1
　　-2 -1 0
　　1 0 0
　　命令托普利兹矩阵
　　函数 toeplitz
　　格式 T = toeplitz(c,r) %生成一个非对称的托普利兹矩阵，将c作为第1列，将r作为第1 行，其余元素与左上角相邻元素相等。
　　T = toeplitz(r) %用向量r生成一个对称的托普利兹矩阵
　　例1-20
　　>> c=[1 2 3 4 5];
　　>> r=[1.5 2.5 3.5 4.5 5.5];
　　>> T="toeplitz"(c,r)
　　T =
　　1 5/2 7/2 9/2 11/2
　　2 1 5/2 7/2 9/2
　　3 2 1 5/2 7/2
　　4 3 2 1 5/2
　　5 4 3 2 1
　　命令 Wilkinson特征值测试阵
　　函数 wilkinson
　　格式 W = wilkinson(n) %返回n阶Wilkinson特征值测试阵
　　例1-21

　　>> W="wilkinson"(4)
　　W =
　　3/2 1 0 0
　　1 1/2 1 0
　　0 1 1/2 1
　　0 0 1 3/2
　　>> W="wilkinson"(7)
　　W =
　　3 1 0 0 0 0 0
　　1 2 1 0 0 0 0
　　0 1 1 1 0 0 0
　　0 0 1 0 1 0 0
　　0 0 0 1 1 1 0
　　0 0 0 0 1 2 1
　　0 0 0 0 0 1 3

复制代码

　　1.2 矩阵运算
　　1.2.1 加、减运算
　　运算符：“+”和“-”分别为加、减运算符。
　　运算规则：对应元素相加、减，即按线性代数中矩阵的“十”，“一”运算进行。
　　例1-22
　　>>A=[1, 1, 1; 1, 2, 3; 1, 3, 6]
　　>>B=[8, 1, 6; 3, 5, 7; 4, 9, 2]
　　>>A+B=A+B
　　>>A-B=A-B
　　结果显示：A+B=
　　9 2 7
　　4 7 10
　　5 12 8
　　A-B=
　　-7 0 -5
　　-2 -3 -4
　　-3 -6 4
　　1.2.2 乘法
　　运算符：*
　　运算规则：按线性代数中矩阵乘法运算进行，即放在前面的矩阵的各行元素，分别与放在后面的矩阵的各列元素对应相乘并相加。
　　1.两个矩阵相乘
　　例1-23
　　>>X= [2 3 4 5;
　　1 2 2 1];
　　>>Y=[0 1 1;
　　1 1 0;
　　0 0 1;
　　1 0 0];
　　Z=X*Y
　　结果显示为：
　　Z=
　　8 5 6
　　3 3 3
　　2.矩阵的数乘：数乘矩阵
　　上例中：a=2*X
　　则显示：a =
　　4 6 8 10
　　2 4 4 2
　　向量的点乘(内积)：维数相同的两个向量的点乘。
　　数组乘法：
　　A.*B表示A与B对应元素相乘。
　　3.向量点积
　　函数 dot
　　格式 C = dot(A,B) %若A、B为向量，则返回向量A与B的点积，A与B长度相同;若为矩阵，则A与B有相同的维数。
　　C = dot(A,B,dim) %在dim维数中给出A与B的点积
　　例 >>X=[-1 0 2];
　　>>Y=[-2 -1 1];
　　>>Z=dot(X, Y)
　　则显示：Z =
　　4
　　还可用另一种算法：
　　sum(X.*Y)
　　ans=
　　4
　　4.向量叉乘
　　在数学上，两向量的叉乘是一个过两相交向量的交点且垂直于两向量所在平面的向量。在Matlab中，用函数cross实现。
　　函数 cross
　　格式 C = cross(A,B) %若A、B为向量，则返回A与B的叉乘，即C=A×B，A、B必须是3个元素的向量;若A、B为矩阵，则返回一个3×n矩阵，其中的列是A与B对应列的叉积，A、B都是3×n矩阵。
　　C = cross(A,B,dim) %在dim维数中给出向量A与B的叉积。A和B必须具有相同的维数，size(A,dim)和size(B,dim)必须是3。
　　例1-24 计算垂直于向量(1, 2, 3)和(4, 5, 6)的向量。
　　>>a=[1 2 3];
　　>>b=[4 5 6];
　　>>c=cross(a,b)
　　结果显示：
　　c=
　　-3 6 -3
　　可得垂直于向量(1, 2, 3)和(4, 5, 6)的向量为±(-3, 6, -3)
　　5.混合积
　　混合积由以上两函数实现：
　　例1-25 计算向量a=(1, 2, 3)、b=(4, 5, 6)和c=(-3, 6, -3) 的混合积
　　解：
　　>>a=[1 2 3]; b=[4 5 6]; c=[-3 6 -3];
　　>>x=dot(a, cross(b, c))
　　结果显示：x =
　　54
　　注意：先叉乘后点乘，顺序不可颠倒。
　　6.矩阵的卷积和多项式乘法
　　函数 conv
　　格式 w = conv(u,v) %u、v为向量，其长度可不相同。
　　说明长度为m的向量序列u和长度为n的向量序列v的卷积(Convolution)定义为：式中：w向量序列的长度为(m+n-1)，当m=n时，

　w(1) = u(1)*v(1)
　　w(2) = u(1)*v(2)+u(2)*v(1)
　　w(3) = u(1)*v(3)+u(2)*v(2)+u(3)*v(1)
　　…
　　w(n) = u(1)*v(n)+u(2)*v(n-1)+ … +u(n)*v(1)
　　…
　　w(2*n-1) = u(n)*v(n)
　　例1-26 展开多项式
　　解： >> w="conv"([1,2,2],conv([1,4],[1,1]))
　　w =
　　1 7 16 18 8
　　>> P="poly2str"(w,'s') %将w表示成多项式
　　P =
　　s^4 + 7 s^3 + 16 s^2 + 18 s + 8

复制代码

　　7.反褶积(解卷)和多项式除法运算
　　函数 deconv
　　格式 [q,r] = deconv(v,u) %多项式v除以多项式u，返回商多项式q和余多项式r。
　　注意：v、u、q、r都是按降幂排列的多项式系数向量。
　　例1-27 卷积
　　>>u = [1 2 3 4]
　　>>v = [10 20 30]
　　>>c = conv(u,v)
　　c =
　　10 40 100 160 170 120
　　则反褶积为
　　>>[q,r] = deconv(c,u)
　　q =
　　10 20 30
　　r =
　　0 0 0 0 0 0
　　8.张量积
　　函数 kron
　　格式 C="kron" (A,B) %A为m×n矩阵，B为p×q矩阵，则C为mp×nq矩阵。
　　说明 A与B的张量积定义为： A B与B A均为mp×nq矩阵，但一般地A B B A。
　　例1-28 求A B。
　　>> A=[1 2;3 4];B=[1 2 3;4 5 6;7 8 9];
　　>> C="kron"(A,B)
　　C =
　　1 2 3 2 4 6
　　4 5 6 8 10 12
　　7 8 9 14 16 18
　　3 6 9 4 8 12
　　12 15 18 16 20 24
　　21 24 27 28 32 36
　　1.2.3 集合运算
　　1.两个集合的交集
　　函数 intersect
　　格式 c = intersect(a,b) %返回向量a、b的公共部分，即c= a∩b。
　　c = intersect(A,B,'rows') %A、B为相同列数的矩阵，返回元素相同的行。
　　[c,ia,ib] = intersect(a,b) %c为a、b的公共元素，ia表示公共元素在a中的位置，ib表示公共元素在b中位置。
　　例1-29

　　>> A=[1 2 3 4;1 2 4 6;6 7 1 4]
　　A =
　　1 2 3 4
　　1 2 4 6
　　6 7 1 4
　　>> B=[1 2 3 8;1 1 4 6;6 7 1 4]
　　B =
　　1 2 3 8
　　1 1 4 6
　　6 7 1 4
　　>> C="intersect"(A,B,'rows')
　　C =
　　6 7 1 4
　　例1-30
　　>> A = [1 9 6 20]; B = [1 2 3 4 6 10 20];
　　>> [c,ia,ib] = intersect(A,B)
　　c =
　　1 6 20
　　ia =
　　1 3 4
　　ib =
　　1 5 7

复制代码

　　2.检测集合中的元素
　　函数 ismember
　　格式 k = ismember(a,S) %当a中元素属于S时，k取1，否则，k取0。
　　k = ismember(A,S,'rows') %A、S有相同的列，返回行相同k取1，不相同取0的列向量。
　　例1-31
　　>> S=[0 2 4 6 8 10 12 14 16 18 20];
　　>> a=[1 2 3 4 5 6];
　　>> k="ismember"(a,S)
　　k =
　　0 1 0 1 0 1 %1表示相同元素的位置
　　例1-32
　　>> A=[1 2 3 4;1 2 4 6;6 7 1 4]
　　>> B=[1 2 3 8;1 1 4 6;6 7 1 4]
　　>> k="ismember"(A,B,'rows')
　　k =
　　0
　　0
　　1 %1表示元素相同的行
　　3.两集合的差
　　函数 setdiff
　　格式 c = setdiff(a,b) %返回属于a但不属于b的不同元素的集合，C = a-b。
　　c = setdiff(A,B,'rows') %返回属于A但不属于B的不同行
　　[c,i] = setdiff(…) %c与前面一致，i表示c中元素在A中的位置。
　　例1-33
　　>> A = [1 7 9 6 20]; B = [1 2 3 4 6 10 20];
　　>> c="setdiff"(A,B)
　　c =
　　7 9
　　例1-34
　　>> A=[1 2 3 4;1 2 4 6;6 7 1 4]
　　>> B=[1 2 3 8;1 1 4 6;6 7 1 4]
　　>> c="setdiff"(A,B,'rows')
　　c =
　　1 2 3 4
　　1 2 4 6
　　4.两个集合交集的非(异或)
　　函数 setxor
　　格式 c = setxor(a,b) %返回集合a、b交集的非
　　c = setxor(A,B,'rows') %返回矩阵A、B交集的非，A、B有相同列数。
　　[c,ia,ib] = setxor(…) %ia、ib表示c中元素分别在a (或A)、b(或B)中位置
　　例1-35

　　>> A=[1 2 3 4];
　　>> B=[2 4 5 8];
　　>> C="setxor"(A,B)
　　C =
　　1 3 5 8
　　例1-36
　　>> A=[1 2 3 4;1 2 4 6;6 7 1 4]
　　A =
　　1 2 3 4
　　1 2 4 6
　　6 7 1 4
　　>> B=[1 2 3 8;1 1 4 6;6 7 1 4]
　　B =
　　1 2 3 8
　　1 1 4 6
　　6 7 1 4
　　>> [C,ia,ib]=setxor(A,B,'rows')
　　C =
　　1 1 4 6
　　1 2 3 4
　　1 2 3 8
　　1 2 4 6
　　ia =
　　1
　　2
　　ib =
　　2
　　1

复制代码

　　5.两集合的并集
　　函数 union
　　格式 c = union(a,b) %返回a、b的并集，即c = a∪b。
　　c = union(A,B,'rows') %返回矩阵A、B不同行向量构成的大矩阵，其中相同行向量只取其一。
　　[c,ia,ib] = union(…) %ia、ib分别表示c中行向量在原矩阵(向量)中的位置
　　例1-37
　　>> A=[1 2 3 4];
　　>> B=[2 4 5 8];
　　>> c="union"(A,B)
　　则结果为
　　c =
　　1 2 3 4 5 8
　　例1-38

　>> A=[1 2 3 4;1 2 4 6]
　　A =
　　1 2 3 4
　　1 2 4 6
　　>> B=[1 2 3 8;1 1 4 6]
　　B =
　　1 2 3 8
　　1 1 4 6
　　>> [c,ia,ib]=union(A,B,'rows')
　　c =
　　1 1 4 6
　　1 2 3 4
　　1 2 3 8
　　1 2 4 6
　　ia =
　　1
　　2
　　ib =
　　2
　　1

复制代码

独善其身 · 发表于 2016-4-7 16:22

　　6.取集合的单值元素

　　函数 unique

　　格式 b = unique (a) %取集合a的不重复元素构成的向量

　　b = unique (A,'rows') %返回A、B不同行元素组成的矩阵

　　[b,i,j] = unique (…) %i、j体现b中元素在原向量(矩阵)中的位置

　　例1-39

　

　>> A=[1 1 2 2 4 4 6 4 6]
　　A =
　　1 1 2 2 4 4 6 4 6
　　>> [c,i,j]=unique(A)
　　c =
　　1 2 4 6
　　i =
　　2 4 8 9
　　j =
　　1 1 2 2 3 3 4 3 4
　　例1-40
　　>> A=[1 2 2 4;1 1 4 6;1 1 4 6]
　　A =
　　1 2 2 4
　　1 1 4 6
　　1 1 4 6
　　>> [c,i,j]=unique(A,'rows')
　　c =
　　1 1 4 6
　　1 2 2 4
　　i =
　　3
　　1
　　j =
　　2
　　1
　　1

　　1.2.4 除法运算

　　Matlab提供了两种除法运算：左除(\)和右除(/)。一般情况下，x=a\b是方程a*x =b的解，而x=b/a是方程x*a=b的解。

　　例：a=[1 2 3; 4 2 6; 7 4 9]

　　b=[4; 1; 2];

　　x=a\b

　　则显示：x=

　　-1.5000

　　2.0000

　　0.5000

　　如果a为非奇异矩阵，则a\b和b/a可通过a的逆矩阵与b阵得到：

　　a\b = inv(a)*b

　　b/a = b*inv(a)

　　数组除法：

　　A./B表示A中元素与B中元素对应相除。

　　1.2.5 矩阵乘方

　　运算符：^

　　运算规则：

　　(1)当A为方阵，P为大于0的整数时，A^P表示A的P次方，即A自乘P次;P为小于0的整数时，A^P表示A-1的P次方。

　　(2)当A为方阵，p为非整数时，则其中V为A的特征向量，为特征值对角矩阵。如果有重根，以上指令不成立。

　　(3)标量的矩阵乘方PA，标量的矩阵乘方定义为式中V，D取自特征值分解AV=AD。

　　(4)标量的数组乘方P.^A，标量的数组乘方定义为数组乘方：A.^P：表示A的每个元素的P次乘方。

　　1.2.6 矩阵函数

　　命令方阵指数

　　函数 expm

　　格式 Y = expm(A) %使用Pade近似算法计算eA，这是一个内部函数，A为方阵。

　　Y="expm1"(A) %使用一个M文件和内部函数相同的算法计算eA

　　Y="expm2"(A) %使用泰勒级数计算eA

　　Y="expm3"(A) %使用特征值和特征向量计算eA

　　命令矩阵的对数

　　函数 logm

　　格式 Y = logm(X) %计算矩阵X的对数，它是expm(X)的反函数。

　　[Y,esterr] = logm(X) %esterr为相对残差的估计值：norm(expm(Y)-X)/norm(X)

　　例1-41

　　>> A=[1 1 0;0 0 2;0 0 -1];

　　>> Y="expm"(A)

　　Y =

　　2.7183 1.7183 1.0862

　　0 1.0000 1.2642

　　0 0 0.3679

　　>> A="logm"(Y)

　　A =

　　1.0000 1.0000 0.0000

　　0 0 2.0000

　　0 0 -1.0000

　　命令方阵的函数

　　函数 funm

　　格式 F = funm(A,fun) %A为方阵，计算由fun指定的A的矩阵函数，fun可以是任意基本函数，如sin、cos等等，例如：funm(A, ’exp’)=expm(A)。

　　[F,esterr] = funm(A,fun) %esterr为结果所产生的相对误差的估计值。

　　命令矩阵的方根

　　函数 sqrtm

　　格式 X = sqrtm(A) %矩阵A的平方根A1/2，相当于X*X=A，求X。若A的特征值有非负实部，则X是唯一的;若A的特征值有负的实部，则X为复矩阵;若A为奇异矩阵，则X不存在。

　　[X,resnorm] = sqrtm(A) % resnorm为结果产生的相对误差

　　[X,alpha,condest] = sqrtm(A) % alpha为稳定因子，condest为结果的条件数的估计值。

　　命令矩阵A的多项式

　　函数 polyvalm

　　格式 polyvalm(P, A) %P为多项式系数向量，方阵A为多项式变量，返回多项式值。

　　1.2.7 矩阵转置

　　运算符：′

　　运算规则：若矩阵A的元素为实数，则与线性代数中矩阵的转置相同。

　　若A为复数矩阵，则A转置后的元素由A对应元素的共轭复数构成。

　　若仅希望转置，则用如下命令：A.′。

　　1.2.8 方阵的行列式

　　函数 det

　　格式 d = det(X) %返回方阵X的多项式的值

　　例1-42

　　>> A=[1 2 3;4 5 6;7 8 9]

　　A =

　　1 2 3

　　4 5 6

　　7 8 9

　　>> D="det"(A)

　　D =

　　0

　　1.2.9 逆与伪逆

　　命令逆

　　函数 inv

　　格式 Y="inv"(X) %求方阵X的逆矩阵。若X为奇异阵或近似奇异阵，将给出警告信息。

　　例1-43 求的逆矩阵

　　方法一

　　>>A=[1 2 3; 2 2 1; 3 4 3];

　　>>Y=inv(A)或Y=A^(-1)

　　则结果显示为

　　Y =

　　1.0000 3.0000 -2.0000

　　-1.5000 -3.0000 2.5000

　　1.0000 1.0000 -1.0000

　　方法二：由增广矩阵进行初等行变换

　　>>B=[1, 2, 3, 1, 0, 0; 2, 2, 1, 0, 1, 0; 3, 4, 3, 0, 0, 1];

　　>>C=rref(B) %化行最简形

　　>>X=C(:, 4:6) %取矩阵C中的A^(-1)部分

　　显示结果如下：

　　C =

　　1.0000 0 0 1.0000 3.0000 -2.0000

　　0 1.0000 0 -1.5000 -3.0000 2.5000

　　0 0 1.0000 1.0000 1.0000 -1.0000

　　X =

　　1.0000 3.0000 -2.0000

　　-1.5000 -3.0000 2.5000

　　1.0000 1.0000 -1.0000

　　例1-44

　　>> A=[2 1 -1;2 1 2;1 -1 1];

　　>> format rat %用有理格式输出

　　>> D="inv"(A)

　　D =

　　1/3 0 1/3

　　0 1/3 -2/3

　　-1/3 1/3 0

　　命令伪逆

　　函数 pinv

　　格式 B = pinv(A) %求矩阵A的伪逆

　　B = pinv(A, tol) %tol为误差：max(size(A))*norm(A)*eps

　　说明当矩阵为长方阵时，方程AX=I和XA=I至少有一个无解，这时A的伪逆能在某种程度上代表矩阵的逆，若A为非奇异矩阵，则pinv(A) = inv(A)。

　　例1-45

　　>> A="magic"(5); %产生5阶魔方阵。

　　>> A="A"(:,1:4) %取5阶魔方阵的前4列元素构成矩阵A。

　　A =

　　17 24 1 8

　　23 5 7 14

　　4 6 13 20

　　10 12 19 21

　　11 18 25 2

　　>> X="pinv"(A) %计算A的伪逆

　　X =

　　-0.0041 0.0527 -0.0222 -0.0132 0.0069

　　0.0437 -0.0363 0.0040 0.0033 0.0038

　　-0.0305 0.0027 -0.0004 0.0068 0.0355

　　0.0060 -0.0041 0.0314 0.0211 -0.0315

　　1.2.10 矩阵的迹

　　函数 trace

　　格式 b="trace" (A) %返回矩阵A的迹，即A的对角线元素之和。

　　1.2.11 矩阵和向量的范数

　　命令向量的范数

　　函数 norm

　　格式 n = norm(X) %X为向量，求欧几里德范数，即。

　　n = norm(X,inf) %求 -范数，即。

　　n = norm(X,1) %求1-范数，即。

　　n = norm(X,-inf) %求向量X的元素的绝对值的最小值，即。

　　n = norm(X, p) %求p-范数，即，所以norm(X,2) = norm(X)。

　　命令矩阵的范数

　　函数 norm

　　格式 n = norm(A) %A为矩阵，求欧几里德范数，等于A的最大奇异值。

　　n = norm(A,1) %求A的列范数，等于A的列向量的1-范数的最大值。

　　n = norm(A,2) %求A的欧几里德范数，和norm(A)相同。

　　n = norm(A,inf) %求行范数，等于A的行向量的1-范数的最大值

　　即：max(sum(abs(A')))。

　　n = norm(A, 'fro' ) %求矩阵A的Frobenius范数，

　　即sqrt(sum(diag(A'*A)))，不能用矩阵p-范数的定义来求。

　　命令范数的估计值

　　函数 normest

　　格式 nrm = normest(A) %矩阵A的2-范数(欧几里德范数)的估计值，相对误差小于106。

　　nrm = normest(A,tol) %tol为指定相对误差

　　[nrm,count] = normest(…) %count给出计算估计值的迭代次数

　　1.2.12 条件数

　　命令矩阵的条件数

　　函数 cond

　　格式 c = cond(X) %求X的2-范数的条件数，即X的最大奇异值和最小奇异值的商。

　　c = cond(X,p) %求p-范数的条件数，p的值可以是1、2、inf或者’fro’。

　　说明线性方程组AX=b的条件数是一个大于或者等于1的实数，用来衡量关于数据中的扰动，也就是A/或b对解X的灵敏度。一个差条件的方程组的条件数很大。条件数的定义为：

　　命令 1-范数的条件数估计

　　函数 condest

　　格式 c = condest (A) %方阵A的1-范数的条件数的下界估值。

　　[c,v] = condest (A) %v为向量，满足，即norm(A*v,1) =norm(A,1)*norm(v,1)/c。

　　[c,v] = condest (A,t) %求上面的c和v，同时显示出关于计算的步骤信息。如果t=1，则计算的每步都显示出来;如果t=-1，则给出商c/rcond(A)。

　　命令矩阵可逆的条件数估值

　　函数 rcond

　　格式 c = rcond(A) %对于差条件矩阵A来说，给出一个接近于0的数;对于好条件矩阵A，则给出一个接近于1的数。

　　命令特征值的条件数

　　函数 condeig

　　格式 c = condeig(A) %返回矩阵A的特征值的条件数

　　[V,D,c] = condeig(A) %D为A的特征值对角阵，V为A的特征向量。

　　1.2.13 矩阵的秩

　　函数 rank

　　格式 k = rank (A) %求矩阵A的秩

　　k = rank (A,tol) %tol为给定误差

　　1.2.14 特殊运算

　　1.矩阵对角线元素的抽取

　　函数 diag

　　格式 X = diag(v,k) %以向量v的元素作为矩阵X的第k条对角线元素，当k=0时，v为X的主对角线;当k>0时，v为上方第k条对角线;当k<0时，v为下方第k条对角线。

　　X = diag(v) %以v为主对角线元素，其余元素为0构成X。

　　v = diag(X,k) %抽取X的第k条对角线元素构成向量v。k=0：抽取主对角线元素;k>0：抽取上方第k条对角线元素;k<0抽取下方第k条对角线元素。

　　v = diag(X) %抽取主对角线元素构成向量v。

　　例1-46

　

　>> v=[1 2 3];
　　>> x="diag"(v,-1)
　　x =
　　0 0 0 0
　　1 0 0 0
　　0 2 0 0
　　0 0 3 0
　　>> A=[1 2 3;4 5 6;7 8 9]
　　A =
　　1 2 3
　　4 5 6
　　7 8 9
　　>> v="diag"(A,1)
　　v =
　　2
　　6

　　2.上三角阵和下三角阵的抽取

　　函数 tril %取下三角部分

　　格式 L = tril(X) %抽取X的主对角线的下三角部分构成矩阵L

　　L = tril(X,k) %抽取X的第k条对角线的下三角部分;k=0为主对角线;k>0为主对角线以上;k<0为主对角线以下。

　　函数 triu %取上三角部分

　　格式 U = triu(X) %抽取X的主对角线的上三角部分构成矩阵U

　　U = triu(X,k) %抽取X的第k条对角线的上三角部分;k=0为主对角线;k>0为主对角线以上;k<0为主对角线以下。

　　例1-47

　

　>> A="ones"(4) %产生4阶全1阵
　　A =
　　1 1 1 1
　　1 1 1 1
　　1 1 1 1
　　1 1 1 1
　　>> L="tril"(A,1) %取下三角部分
　　L =
　　1 1 0 0
　　1 1 1 0
　　1 1 1 1
　　1 1 1 1
　　>> U="triu"(A,-1) %取上三角部分
　　U =
　　1 1 1 1
　　1 1 1 1
　　0 1 1 1
　　0 0 1 1

　　3.矩阵的变维

　　矩阵的变维有两种方法，即用“：”和函数“reshape”，前者主要针对2个已知维数矩阵之间的变维操作;而后者是对于一个矩阵的操作。

　　(1)“：”变维

　　例1-48

　

　> A=[1 2 3 4 5 6;6 7 8 9 0 1]
　　A =
　　1 2 3 4 5 6
　　6 7 8 9 0 1
　　>> B="ones"(3,4)
　　B =
　　1 1 1 1
　　1 1 1 1
　　1 1 1 1
　　>> B(:)=A(:)
　　B =
　　1 7 4 0
　　6 3 9 6
　　2 8 5 1

　　(2)Reshape函数变维

　　格式 B = reshape(A,m,n) %返回以矩阵A的元素构成的m×n矩阵B

　　B = reshape(A,m,n,p,…) %将矩阵A变维为m×n×p×…

　　B = reshape(A,[m n p…]) %同上

　　B = reshape(A,siz) %由siz决定变维的大小，元素个数与A中元素个数

　　相同。

　　例1-49 矩阵变维

　　>> a=[1:12];

　　>> b="reshape"(a,2,6)

　　b =

　　1 3 5 7 9 11

　　2 4 6 8 10 12

　　4.矩阵的变向

　　(1)矩阵旋转

　　函数 rot90

　　格式 B = rot90 (A) %将矩阵A逆时针方向旋转90°

　　B = rot90 (A,k) %将矩阵A逆时针方向旋转(k×90°)，k可取正负整数。

　　例1-50

　　>> A=[1 2 3;4 5 6;7 8 9]

　　A =

　　1 2 3

　　4 5 6

　　7 8 9

　　>> Y1=rot90(A),Y2=rot90(A,-1)

　　Y1 = %逆时针方向旋转

　　3 6 9

　　2 5 8

　　1 4 7

　　Y2 = %顺时针方向旋转

　　7 4 1

　　8 5 2

　　9 6 3

　　(2)矩阵的左右翻转

　　函数 fliplr

　　格式 B = fliplr(A) %将矩阵A左右翻转

　　(3)矩阵的上下翻转

　　函数 flipud

　　格式 B = flipud(A) %将矩阵A上下翻转

　　例1-51

　　>> A=[1 2 3;4 5 6]

　　A =

　　1 2 3

　　4 5 6

　　>> B1=fliplr(A),B2=flipud(A)

　　B1 =

　　3 2 1

　　6 5 4

　　B2 =

　　4 5 6

　　1 2 3

　　(4)按指定维数翻转矩阵

　　函数 flipdim

　　格式 B = flipdim(A,dim) % flipdim(A,1) = flipud(A)，并且flipdim(A,2)=fliplr(A)。

　　例1-52

　　>> A=[1 2 3;4 5 6]

　　A =

　　1 2 3

　　4 5 6

　　>> B1=flipdim(A,1),B2=flipdim(A,2)

　　B1 =

　　4 5 6

　　1 2 3

　　B2 =

　　3 2 1

　　6 5 4

　　(5)复制和平铺矩阵

　　函数 repmat

　　格式 B = repmat(A,m,n) %将矩阵A复制m×n块，即B由m×n块A平铺而成。

　　B = repmat(A,[m n]) %与上面一致

　　B = repmat(A,[m n p…]) %B由m×n×p×…个A块平铺而成

　　repmat(A,m,n) %当A是一个数a时，该命令产生一个全由a组成的m×n矩阵。

　　例1-53

　　>> A=[1 2;5 6]
　　A =
　　1 2
　　5 6
　　>> B="repmat"(A,3,4)
　　B =
　　1 2 1 2 1 2 1 2
　　5 6 5 6 5 6 5 6
　　1 2 1 2 1 2 1 2
　　5 6 5 6 5 6 5 6
　　1 2 1 2 1 2 1 2
　　5 6 5 6 5 6 5 6

　　5.矩阵的比较关系

　　矩阵的比较关系是针对于两个矩阵对应元素的，所以在使用关系运算时，首先应该保证两个矩阵的维数一致或其中一个矩阵为标量。关系运算是对两个矩阵的对应运算进行比较，若关系满足，则将结果矩阵中该位置元素置为1，否则置0。

　　MATLAB的各种比较关系运算有见表1-1。

　　表1-1

　　运算符含义运算符含义

　　> 大于关系 < 小于关系

　　= = 等于关系 >= 大于或等于关系

　　<= 小于或等于关系 ~ = 不等于关系

　　例1-54

　　>> A=[1 2 3 4;5 6 7 8];B=[0 2 1 4;0 7 7 2];
　　>> C1=A==B, C2=A>=B, C3=A~=B
　　C1 =
　　0 1 0 1
　　0 0 1 0
　　C2 =
　　1 1 1 1
　　1 0 1 1
　　C3 =
　　1 0 1 0
　　1 1 0 1

　　6.矩阵元素的数据变换

　　对于小数构成的矩阵A来说，如果我们想对它取整数，有以下几种方法：

　　(1)按-∞方向取整

　　函数 floor

　　格式 floor(A) %将A中元素按-∞方向取整，即取不足整数。

　　(2)按+∞方向取整

　　函数 ceil

　　格式 ceil(A) %将A中元素按+∞方向取整，即取过剩整数。

　　(3)四舍五入取整

　　函数 round

　　格式 round (A) %将A中元素按最近的整数取整，即四舍五入取整。

　　(4)按离0近的方向取整

　　函数 fix

　　格式 fix (A) %将A中元素按离0近的方向取整

　　例1-55

　

　>> A="-1".5+4*rand(3)
　　A =
　　2.3005 0.4439 0.3259
　　-0.5754 2.0652 -1.4260
　　0.9274 1.5484 1.7856
　　>> B1=floor(A),B2=ceil(A),B3=round(A),B4=fix(A)
　　B1 =
　　2 0 0
　　-1 2 -2
　　0 1 1
　　B2 =
　　3 1 1
　　0 3 -1
　　1 2 2
　　B3 =
　　2 0 0
　　-1 2 -1
　　1 2 2
　　B4 =
　　2 0 0
　　0 2 -1
　　0 1 1

　　(5)矩阵的有理数形式

　　函数 rat

　　格式 [n,d]=rat (A) %将A表示为两个整数矩阵相除，即A=n./d。

　　例1-56 对于上例中的A

　　>> [n,d]=rat(A)

　　n =

　　444 95 131

　　-225 2059 -472

　　166 48 1491

　　d =

　　193 214 402

　　391 997 331

　　179 31 835

　　(6)矩阵元素的余数

　　函数 rem

　　格式 C = rem (A, x) %表示A矩阵除以模数x后的余数。若x=0，则定义rem(A, 0)=NaN，若x≠0，则整数部分由fix(A./x)表示，余数C=A-x.*fix (A./x)。允许模x为小数。

　　7.矩阵逻辑运算

　　设矩阵A和B都是m×n矩阵或其中之一为标量，在MATLAB中定义了如下的逻辑运算：

　　(1)矩阵的与运算

　　格式 A&B或and(A, B)

　　说明 A与B对应元素进行与运算，若两个数均非0，则结果元素的值为1，否则为0。

　　(2)或运算

　　格式 A|B或or(A, B)

　　说明 A与B对应元素进行或运算，若两个数均为0，则结果元素的值为0，否则为1。

　　(3)非运算

　　格式 ~A或not (A)

　　说明若A的元素为0，则结果元素为1，否则为0。

　　(4)异或运算

　　格式 xor (A,B)

　　说明 A与B对应元素进行异或运算，若相应的两个数中一个为0，一个非0，则结果为0，否则为1。

　　例1-57

　　>> A=[0 2 3 4;1 3 5 0],B=[1 0 5 3;1 5 0 5]
　　A =
　　0 2 3 4
　　1 3 5 0
　　B =
　　1 0 5 3
　　1 5 0 5
　　>> C1=A&B,C2=A|B,C3=~A,C4=xor(A,B)
　　C1 =
　　0 0 1 1
　　1 1 0 0
　　C2 =
　　1 1 1 1
　　1 1 1 1
　　C3 =
　　1 0 0 0
　　0 0 0 1
　　C4 =
　　1 1 0 0
　　0 0 1 1

　　1.2.15 符号矩阵运算

　　1.符号矩阵的四则运算

　　Matlab 6.x 抛弃了在4.2版中为符号矩阵设计的复杂函数形式，把符号矩阵的四则运算简化为与数值矩阵完全相同的运算方式，其运算符为：加(+)、减(-)、乘(×)、除(/、\)等或：符号矩阵的和(symadd)、差(symsub)、乘 (symmul)。

　　例1-58 >>

　　>> ;

　　>>C=B-A

　　>>D=a\b

　　则显示：

　　C=

　　x-1/x 1-1/(x+1)

　　x+2-1/(x+2) -1/(x+3)

　　D=

　　-6*x-2*x^3-7*x^2 1/2*x^3+x+3/2*x^2

　　6+2*x^3+10*x^2+14*x -2*x^2-3/2*x-1/2*x^3

　　2.其他基本运算

　　符号矩阵的其他一些基本运算包括转置(')、行列式(det)、逆(inv)、秩(rank)、幂(^)和指数(exp和expm)等都与数值矩阵相同

　　3.将数值矩阵转化为符号矩阵

　　函数 sym

　　格式 B="sym"(A) %将A转化为符号矩阵B

　　例1-59

　　>> A=[2/3,sqrt(2),0.222;1.4,1/0.23,log(3)]

　　A =

　　0.6667 1.4142 0.2220

　　1.4000 4.3478 1.0986

　　>> B="sym"(A)

　　B =

　　[ 2/3, sqrt(2), 111/500]

　　[ 7/5, 100/23, 4947709893870346*2^(-52)]

　　4.符号矩阵的索引与修改

　　符号矩阵的索引与修改同数值矩阵的索引与修改完全相同，即用矩阵的坐标括号表达式实现。

　　例1-60 对上例中的矩阵B

　　>> B(2,3) %矩阵的索引

　　ans =

　　4947709893870346*2^(-52)

　　>> B(2,3)='log(7)' %矩阵的修改

　　B =

　　[ 2/3, sqrt(2), 111/500]

　　[ 7/5, 100/23, log(7)]

　　5.符号矩阵的简化

　　符号工具箱中提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作函数。

　　(1)因式分解

　　函数 factor

　　格式 factor(s) %符号表达式s的因式分解函数

　　说明 S为符号矩阵或符号表达式，常用于多项式的因式分解。

　　例1-61 将x 9-1分解因式

　　在Matlab命令窗口键入：

　　syms x

　　factor(x^9-1)

　　则显示：ans =

　　(x-1)*(x^2+x+1)*(x6+x^3+1)

　　例1-62 问“入”取何值时，齐次方程组有非0解?

　　解：在Matlab编辑器中建立M文件：

　　syms k

　　A=[1-k -2 4;2 3-k 1;1 1 1-k];

　　D=det(A)

　　factor(D)

　　其结果显示如下：

　　D =

　　-6*k+5*k^2-k^3

　　ans =

　　-k*(k-2)*(-3+k)

　　从而得到：当k=0、k=2或k=3时，原方程组有非0解。

　　(2)符号矩阵的展开

　　函数 expand

　　格式：expand(s) %符号表达式s的展开函数

　　说明：s为符号矩阵或表达式。常用在多项式的因式分解中，也常用于三角函数，指数函数和对数函数的展开中。

　　例1-63 将(x+1)3、sin(x+y)展开

　　在Matlab编辑器中建立M文件：

　　syms x y

　　p=expand((x+1)^3)

　　q=expand(sin(x+y))

　　则结果显示为

　　p =

　　x^3+3*x^2+3*x+1

　　q =

　　sin(x)*cos(y)+cos(x)*sin(y)

　　(3)同类式合并

　　函数 Collect

　　格式 Collect(s,v) %将s中的变量v的同幂项系数合并

　　Collect(s) % s是矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合并。

　　(4)符号简化

　　函数 simple或simplify %寻找符号矩阵或符号表达式的最简型

　　格式 simple (s) % s是矩阵或表达式

　　[R,how]=simple (s) %R为返回的最简形，how为简化过程中使用的主要方法。

　　说明 Simple(s)将表达式s的长度化到最短。若还想让表达式更加精美，可使用函数Pretty。

　　格式 Pretty(s) %使表达式s更加精美

　　例1-64 计算行列式的值。

　　在Matlab编辑器中建立M文件：

　　syms a b c d

　　A=[1 1 1 1;a b c d;a^2 b^2 c^2 d^2;a^4 b^4 c^4 d^4];

　　d1=det(A)

　　d2=simple(d1) %化简表达式d1

　　pretty(d2) %让表达式d2符合人们的书写习惯

　　则显示结果如下：

　　d1 =

　　b*c^2*d^4-b*d^2*c^4-b^2*c*d^4+b^2*d*c^4+b^4*c*d^2-b^4*d*c^2-a*c^2*d^4+a*d^2*c^4+a*b^2*d^4-a*b^2*c^4-a*b^4*d^2+a*b^4*c^2+a^2*c*d^4-a^2*d*c^4-a^2*b*d^4+a^2*b*c^4+a^2*b^4*d-a^2*b^4*c-a^4*c*d^2+a^4*d*c^2+a^4*b*d^2-a^4*b*c^2-a^4*b^2*d+a^4*b^2*c

　　d2 =

　　(-d+c)*(b-d)*(b-c)*(-d+a)*(a-c)*(a-b)*(a+c+d+b)

　　(-d+c)(b-d)(b-c)(-d+a)(a-c)(a-b)(a+c+d+b)

　　1.2.16 矩阵元素个数的确定

　　函数 numel

　　格式 n = numel(a) %计算矩阵A中元素的个数

　　例1-65

　　>> A=[1 2 3 4;5 6 7 8];

　　>> n="numel"(A)

　　n =

　　8

　　1.3 矩阵分解

　　1.3.1 Cholesky分解

　　函数 chol

　　格式 R = chol(X) %如果X为n阶对称正定矩阵，则存在一个实的非奇异上三角阵R，满足R'*R = X;若X非正定，则产生错误信息。

　　[R,p] = chol(X) %不产生任何错误信息，若X为正定阵，则p=0，R与上相同;若X非正定，则p为正整数，R是有序的上三角阵。

　　例1-66

　　>> X="pascal"(4) %产生4阶pascal矩阵

　　X =

　　1 1 1 1

　　1 2 3 4

　　1 3 6 10

　　1 4 10 20

　　>> [R,p]=chol(X)

　　R =

　　1 1 1 1

　　0 1 2 3

　　0 0 1 3

　　0 0 0 1

　　p =

　　0

　　1.3.2 LU分解

　　矩阵的三角分解又称LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。

　　函数 lu

　　格式 [L,U] = lu(X) %U为上三角阵，L为下三角阵或其变换形式，满足LU=X。

　　[L,U,P] = lu(X) %U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，满足LU=PX。

　　例1-67

　

　>> A=[1 2 3;4 5 6;7 8 9];
　　>> [L,U]=lu(A)
　　L =
　　0.1429 1.0000 0
　　0.5714 0.5000 1.0000
　　1.0000 0 0
　　U =
　　7.0000 8.0000 9.0000
　　0 0.8571 1.7143
　　0 0 0.0000
　　>> [L,U,P]=lu(A)
　　L =
　　1.0000 0 0
　　0.1429 1.0000 0
　　0.5714 0.5000 1.0000
　　U =
　　7.0000 8.0000 9.0000
　　0 0.8571 1.7143
　　0 0 0.0000
　　P =
　　0 0 1
　　1 0 0
　　0 1 0

　　1.3.3 QR分解

　　将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。

　　函数 qr

　　格式 [Q,R] = qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。

　　[Q,R,E] = qr(A) %求得正交矩阵Q和上三角阵R，E为单位矩阵的变换形式，R的对角线元素按大小降序排列，满足AE=QR。

　　[Q,R] = qr(A,0) %产生矩阵A的“经济大小”分解

　　[Q,R,E] = qr(A,0) %E的作用是使得R的对角线元素降序，且Q*R=A(:, E)。

　　R = qr(A) %稀疏矩阵A的分解，只产生一个上三角阵R，满足R'*R = A'*A，这种方法计算A'*A时减少了内在数字信息的损耗。

　　[C,R] = qr(A,b) %用于稀疏最小二乘问题：minimize||Ax-b||的两步解：[C,R] = qr(A,b)，x = R\c。

　　R = qr(A,0) %针对稀疏矩阵A的经济型分解

　　[C,R] = qr(A,b,0) %针对稀疏最小二乘问题的经济型分解

　　例1-68

　　>>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.4461

　　R =

　　-12.8841 -14.5916 -16.2992

　　0 -1.0413 -2.0826

　　0 0 0.0000

　　0 0 0

　　函数 qrdelete

　　格式 [Q,R] = qrdelete(Q,R,j) %返回将矩阵A的第j列移去后的新矩阵的qr分解

　　例1-69

　　>> A=[-149 -50 -154;537 180 546;-27 -9 -25];

　　>> [Q,R]=qr(A)

　　Q =

　　-0.2671 -0.7088 0.6529

　　0.9625 -0.1621 0.2176

　　-0.0484 0.6865 0.7255

　　R =

　　557.9418 187.0321 567.8424

　　0 0.0741 3.4577

　　0 0 0.1451

　　>> [Q,R]=qrdelete(Q,R,3) %将A的第3列去掉后进行qr分解。

　　Q =

　　-0.2671 -0.7088 0.6529

　　0.9625 -0.1621 0.2176

　　-0.0484 0.6865 0.7255

　　R =

　　557.9418 187.0321

　　0 0.0741

　　0 0

　　函数 qrinsert

　　格式 [Q,R] = qrinsert(Q,R,j,x) %在矩阵A中第j列插入向量x后的新矩阵进行qr分解。若j大于A的列数，表示在A的最后插入列x。

　　例1-70

　　>> A=[-149 -50 -154;537 180 546;-27 -9 -25];

　　>> x=[35 10 7]';

　　>> [Q,R]=qrinsert(Q,R,4,x)

　　Q =

　　-0.2671 -0.7088 0.6529

　　0.9625 -0.1621 0.2176

　　-0.0484 0.6865 0.7255

　　R =

　　557.9418 187.0321 567.8424 -0.0609

　　0 0.0741 3.4577 -21.6229

　　0 0 0.1451 30.1073

　　1.3.4 Schur分解

　　函数 schur

　　格式 T = schur(A) %产生schur矩阵T，即T的主对角线元素为特征值的三角阵。

　　T = schur(A,flag) %若A有复特征根，则flag='complex'，否则flag='real'。

　　[U,T] = schur(A,…) %返回正交矩阵U和schur矩阵T，满足A = U*T*U'。

　　例1-71

　　>> H = [ -149 -50 -154; 537 180 546; -27 -9 -25 ];

　　>> [U,T]=schur(H)

　　U =

　　0.3162 -0.6529 0.6882

　　-0.9487 -0.2176 0.2294

　　0.0000 0.7255 0.6882

　　T =

　　1.0000 -7.1119 -815.8706

　　0 2.0000 -55.0236

　　0 0 3.0000

　　1.3.5 实Schur分解转化成复Schur分解

　　函数 rsf2csf

　　格式 [U,T] = rsf2csf (U,T) %将实舒尔形式转化成复舒尔形式

　　例1-72

　　>> A=[1 1 1 3;1 2 1 1;1 1 3 1;-2 1 1 4];

　　>> [u,t]=schur (A)

　　u =

　　-0.4916 -0.4900 -0.6331 -0.3428

　　-0.4980 0.2403 -0.2325 0.8001

　　-0.6751 0.4288 0.4230 -0.4260

　　-0.2337 -0.7200 0.6052 0.2466

　　t =

　　4.8121 1.1972 -2.2273 -1.0067

　　0 1.9202 -3.0485 -1.8381

　　0 0.7129 1.9202 0.2566

　　0 0 0 1.3474

　　>> [U,T]=rsf2csf (u,t)

　　U =

　　-0.4916 -0.2756 - 0.4411i 0.2133 + 0.5699i -0.3428

　　-0.4980 -0.1012 + 0.2163i -0.1046 + 0.2093i 0.8001

　　-0.6751 0.1842 + 0.3860i -0.1867 - 0.3808i -0.4260

　　-0.2337 0.2635 - 0.6481i 0.3134 - 0.5448i 0.2466

　　T =

　　4.8121 -0.9697 + 1.0778i -0.5212 + 2.0051i -1.0067

　　0 1.9202 + 1.4742i 2.3355 0.1117 + 1.6547i

　　0 0 1.9202 - 1.4742i 0.8002 + 0.2310i

　　0 0 0 1.3474

　　1.3.6 特征值分解

　　函数 eig

　　格式 d = eig(A) %求矩阵A的特征值d，以向量形式存放d。

　　d = eig(A,B) %A、B为方阵，求广义特征值d，以向量形式存放d。

　　[V,D] = eig(A) %计算A的特征值对角阵D和特征向量V，使AV=VD成立。

　　[V,D] = eig(A,'nobalance') %当矩阵A中有与截断误差数量级相差不远的值时，该指令可能更精确。'nobalance'起误差调节作用。

　　[V,D] = eig(A,B) %计算广义特征值向量阵V和广义特征值阵D，满足AV=BVD。

　　[V,D] = eig(A,B,flag) % 由flag指定算法计算特征值D和特征向量V，flag的可能值为：'chol' 表示对B使用Cholesky分解算法，这里A为对称Hermitian矩阵，B为正定阵。'qz' 表示使用QZ算法，这里A、B为非对称或非Hermitian矩阵。

独善其身 · 发表于 2016-4-7 16:23

　　1.3.7 奇异值分解

　　函数 svd

　　格式 s = svd (X) %返回矩阵X的奇异值向量

　　[U,S,V] = svd (X) %返回一个与X同大小的对角矩阵S，两个酉矩阵U和V，且满足= U*S*V'。若A为m×n阵，则U为m×m阵，V为n×n阵。奇异值在S的对角线上，非负且按降序排列。

　　[U,S,V] = svd (X,0) %得到一个“有效大小”的分解，只计算出矩阵U的前n列，矩阵S的大小为n×n。

　　例1-73

　　>> A=[1 2;3 4;5 6;7 8];

　　>> [U,S,V]=svd(A)

　　U =

　　-0.1525 -0.8226 -0.3945 -0.3800

　　-0.3499 -0.4214 0.2428 0.8007

　　-0.5474 -0.0201 0.6979 -0.4614

　　-0.7448 0.3812 -0.5462 0.0407

　　S =

　　14.2691 0

　　0 0.6268

　　0 0

　　0 0

　　V =

　　-0.6414 0.7672

　　-0.7672 -0.6414

　　>> [U,S,V]=svd(A,0)

　　U =

　　-0.1525 -0.8226

　　-0.3499 -0.4214

　　-0.5474 -0.0201

　　-0.7448 0.3812

　　S =

　　14.2691 0

　　0 0.6268

　　V =

　　-0.6414 0.7672

　　-0.7672 -0.6414

　　1.3.8 广义奇异值分解

　　函数 gsvd

　　格式 [U,V,X,C,S] = gsvd(A,B) %返回酉矩阵U和V、一个普通方阵X、非负对角矩阵C和S，满足A = U*C*X'，B = V*S*X'，C'*C + S'*S = I (I为单位矩阵);A和B的列数必须相同，行数可以不同。

　　[U,V,X,C,S] = gsvd(A,B,0) %含义与前面相似

　　sigma = gsvd (A,B) %返回广义奇异值sigma

　　例1-74

　　>> A="reshape"(1:12,3,4) %产生3行4列矩阵，元素由1，2，…，12构成。

　　A =

　　1 4 7 10

　　2 5 8 11

　　3 6 9 12

　　>> B="magic"(4) %产生4阶魔方阵

　　B =

　　16 2 3 13

　　5 11 10 8

　　9 7 6 12

　　4 14 15 1

　　>> [U,V,X,C,S]=gsvd(A,B)

　　U =

　　0.4082 0.7071 0.5774

　　-0.8165 0.0000 0.5774

　　0.4082 -0.7071 0.5774

　　V =

　　0.2607 -0.7950 -0.5000 0.2236

　　-0.4029 0.3710 -0.5000 0.6708

　　-0.5452 -0.0530 -0.5000 -0.6708

　　0.6874 0.4770 -0.5000 -0.2236

　　X =

　　0 -9.4340 -17.0587 3.4641

　　1.8962 8.7980 -17.0587 8.6603

　　3.7924 8.1620 -17.0587 13.8564

　　-5.6885 -7.5260 -17.0587 19.0526

　　C =

　　0 0.0000 0 0

　　0 0 0.0829 0

　　0 0 0 1.0000

　　S =

　　1.0000 0 0 0

　　0 1.0000 0 0

　　0 0 0.9966 0

　　0 0 0 0.0000

　　1.3.9 特征值问题的QZ分解

　　函数 qz

　　格式 [AA,BB,Q,Z,V] = qz(A,B) %A、B为方阵，产生上三角阵AA和BB，正交矩阵Q、Z或其列变换形式，V为特征向量阵。且满足：Q*A*Z= AA 和Q*B*Z = BB。

　　[AA,BB,Q,Z,V] = qz(A,B,flag) %产生由flag决定的分解结果，flag取值为：'complex'：表示复数分解(默认)，取值为'real'：表示实数分解。

　　1.3.10 海森伯格形式的分解

　　如果矩阵H的第一子对角线下元素都是0，则H为海森伯格(Hessenberg)矩阵。如果矩阵是对称矩阵，则它的海森伯格形式是对角三角阵。MATLAB可以通过相似变换将矩阵变换成这种形式。

　　函数 hess

　　格式 H = hess(A) %返回矩阵A的海森伯格形式

　　[P,H] = hess(A) %P为酉矩阵，满足：A = PHP' 且P'P = eye(size(A))。

　　例1-75

　　>> A=[-149 -50 -154;537 180 546;-27 -9 -25];

　　>> [P,H]=hess(A)

　　P =

　　1.0000 0 0

　　0 -0.9987 0.0502

　　0 0.0502 0.9987

　　H =

　　-149.0000 42.2037 -156.3165

　　-537.6783 152.5511 -554.9272

　　0 0.0728 2.4489

　　H的第一子对角元素是H(3,1)=0。

　　1.4 线性方程的组的求解

　　我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：

　　若系数矩阵的秩r=n(n为方程组中未知变量的个数)，则有唯一解;

　　若系数矩阵的秩r

　　线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解;其特解的求法属于解的第一类问题，通解部分属第二类问题。

　　1.4.1 求线性方程组的唯一解或特解(第一类问题)

　　这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 —— 直接法;另一类是解大型稀疏矩阵 —— 迭代法。

　　1.利用矩阵除法求线性方程组的特解(或一个解)

　　方程：AX=b

　　解法：X=A\b

　　例1-76 求方程组的解。

　　解：

　　>>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) %求秩

　　X=A\B %求解

　　运行后结果如下

　　R_A =

　　5

　　X =

　　2.2662

　　-1.7218

　　1.0571

　　-0.5940

　　0.3188

　　这就是方程组的解。

　　或用函数rref求解：

　　>> C=[A,B] %由系数矩阵和常数列构成增广矩阵C

　　>> 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的最后一列元素就是所求之解。

　　例1-77 求方程组的一个特解。

　　解：

　　>>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];

　　>>B=[1 4 0]';

　　>>X=A\B %由于系数矩阵不满秩，该解法可能存在误差。

　　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.2500

　　0 1.0000 -1.5000 -1.7500 -0.2500

　　0 0 0 0 0

　　由此得解向量X=[1.2500 – 0.2500 0 0]’(一个特解)。

　　2.利用矩阵的LU、QR和cholesky分解求方程组的解

　　(1)LU分解：

　　LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式(行交换)和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。

　　则：A*X=b 变成L*U*X=b

　　所以X=U\(L\b) 这样可以大大提高运算速度。

　　命令 [L，U]=lu (A)

　　例1-78 求方程组的一个特解。

　　解：

　　>>A=[4 2 -1;3 -1 2;11 3 0];

　　>>B=[2 10 8]';

　　>>D=det(A)

　　>>[L,U]=lu(A)

　　>>X=U\(L\B)

　　显示结果如下：

　　D =

　　0

　　L =

　　0.3636 -0.5000 1.0000

　　0.2727 1.0000 0

　　1.0000 0 0

　　U =

　　11.0000 3.0000 0

　　0 -1.8182 2.0000

　　0 0 0.0000

　　Warning: Matrix is close to singular or badly scaled.

　　Results may be inaccurate. RCOND = 2.018587e-017.

　　> In D:\Matlab\pujun\lx0720.m at line 4

　　X =

　　1.0e+016 *

　　-0.4053

　　1.4862

　　1.3511

　　说明结果中的警告是由于系数行列式为零产生的。可以通过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"\(Q\b)

　　上例中 [Q, R]=qr(A)

　　X=R\(Q\B)

　　说明这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。

　　1.4.2 求线性齐次方程组的通解

　　在Matlab中，函数null用来求解零空间，即满足A?X=0的解空间，实际上是求出解空间的一组基(基础解系)。

　　格式 z = null % z的列向量为方程组的正交规范基，满足。

　　% z的列向量是方程AX=0的有理基

　　例1-79 求解方程组的通解：

　　解：

　　>>A=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3];

　　>>format rat %指定有理式格式输出

　　>>B=null(A,'r') %求解空间的有理基

　　运行后显示结果如下：

　　B =

　　2 5/3

　　-2 -4/3

　　1 0

　　0 1

　　或通过行最简行得到基：

　　>> B="rref"(A)

　　B =

　　1.0000 0 -2.0000 -1.6667

　　0 1.0000 2.0000 1.3333

　　0 0 0 0

　　即可写出其基础解系(与上面结果一致)。

　　写出通解：

　　syms k1 k2

　　X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解

　　pretty(X) %让通解表达式更加精美

　　运行后结果如下：

　　X =

　　[ 2*k1+5/3*k2]

　　[ -2*k1-4/3*k2]

　　[ k1]

　　[ k2]

　　% 下面是其简化形式

　　[2k1 + 5/3k2 ]

　　[ ]

　　[-2k1 - 4/3k2]

　　[ ]

　　[ k1 ]

　　[ ]

　　[ k2 ]

　　1.4.3 求非齐次线性方程组的通解

　　非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。

　　因此，步骤为：

　　第一步：判断AX=b是否有解，若有解则进行第二步

　　第二步：求AX=b的一个特解

　　第三步：求AX=0的通解

　　第四步：AX=b的通解= AX="0的通解"+AX=b的一个特解。

　　例1-80 求解方程组

　　解：在Matlab中建立M文件如下：

　　A=[1 -2 3 -1;3 -1 5 -3;2 1 2 -2];

　　b=[1 2 3]';

　　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

　　X="A"\b %求特解

　　C="null"(A,'r') %求AX=0的基础解系

　　else X='equition no solve' %判断无解

　　end

　　运行后结果显示：

　　R_A =

　　2

　　R_B =

　　3

　　X =

　　equition no solve

　　说明该方程组无解

　　例1-81 求解方程组的通解：

　　解法一：在Matlab编辑器中建立M文件如下：

　　A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];

　　b=[1 4 0]';

　　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

　　X="A"\b

　　C="null"(A,'r')

　　else X='Equation has no solves'

　　end

　　运行后结果显示为：

　　R_A =

　　2

　　R_B =

　　2

　　Warning: Rank deficient, rank = 2 tol = 8.8373e-015.

　　> In D:\Matlab\pujun\lx0723.m at line 11

　　X =

　　0

　　0

　　-8/15

　　3/5

　　C =

　　3/2 -3/4

　　3/2 7/4

　　1 0

　　0 1

　　所以原方程组的通解为X=k1 +k2 +

　　解法二：用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/4

　　0 0 0 0 0

　　对应齐次方程组的基础解系为：，非齐次方程组的特解为：所以，原方程组的通解为：X=k1 +k2 + 。

　　1.4.4 线性方程组的LQ解法

　　函数 symmlq

　　格式 x = symmlq(A,b) %求线性方程组AX=b的解X。A必须为n阶对称方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息;如果收敛失败，将给出警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。

　　symmlq(A,b,tol) %指定误差tol，默认值是1e-6。

　　symmlq(A,b,tol,maxit) %maxit指定最大迭代次数

　　symmlq(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子

　　symmlq(A,b,tol,maxit,M1,M2) %M=M1×M2

　　symmlq(A,b,tol,maxit,M1,M2,x0) %x0为初始估计值，默认值为0。

　　[x,flag] = symmlq(A,b,…) %flag的取值为：0表示在指定迭代次数之内按要求精度收敛;1表示在指定迭代次数内不收敛;2表示M为坏条件的预处理因子;3表示两次连续迭代完全相同;4表示标量参数太小或太大;5表示预处理因子不是对称正定的。

　　[x,flag,relres] = symmlq(A,b,…) % relres表示相对误差norm(b-A*x)/norm(b)

　　[x,flag,relres,iter] = symmlq(A,b,…) %iter表示计算x的迭代次数

　　[x,flag,relres,iter,resvec] = symmlq(A,b,…) %resvec表示每次迭代的残差：norm(b-A*x0)

　　[x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,…) %resveccg表示每次迭代共轭梯度残差的范数

　　1.4.5 双共轭梯度法解方程组

　　函数 bicg

　　格式 x = bicg(A,b) %求线性方程组AX=b的解X。A必须为n阶方阵，b为n元列向量。A可以是由afun定义并返回A*X的函数。如果收敛，将显示结果信息;如果收敛失败，将给出警告信息并显示相对残差norm(b-A*x)/norm(b)和计算终止的迭代次数。

　　bicg(A,b,tol) %指定误差tol，默认值是1e-6。

　　bicg(A,b,tol,maxit) %maxit指定最大迭代次数

　　bicg(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子

　　bicg(A,b,tol,maxit,M1,M2) %M=M1×M2

　　bicg(A,b,tol,maxit,M1,M2,x0) %x0为初始估计值，默认值为0。

　　[x,flag] = bicg(A,b,…) %flag的取值为：0表示在指定迭代次数之内按要求精度收敛;1表示在指定迭代次数内不收敛;2表示M为坏条件的预处理因子;3表示两次连续迭代完全相同;4表示标量参数太小或太大。

　　[x,flag,relres] = bicg(A,b,…) % relres表示相对误差norm(b-A*x)/norm(b)

　　[x,flag,relres,iter] = bicg(A,b,…) %iter表示计算x的迭代次数

　　[x,flag,relres,iter,resvec] = bicg(A,b,…) %resvec表示每次迭代的残差：norm(b-A*x0)

　　例1-83 调用MATLAB6.0数据文件west0479。

　　>> load west0479

　　>> A="west0479"; %将数据取为系数矩阵A。

　　>> b="sum" (A,2); %将A的各行求和，构成一列向量。

　　>> X="A"\b; %用“\”求AX=b的解。

　　>> norm(b-A*X)/norm(b) %计算解的相对误差。

　　ans =

　　1.2454e-017

　　>> [x,flag,relres,iter,resvec] = bicg(A,b) %用bicg函数求解。

　　x = (全为0，由于太长，不显示出来)

　　flag =

　　1 %表示在默认迭代次数(20次)内不收敛。

　　relres = %相对残差relres = norm(b-A*x)/norm(b) =norm(b)/norm(b) = 1。

　　1

　　iter = %表明解法不当，使得初始估计值0向量比后来所有迭代值都好。

　　0

　　resvec = (略) %每次迭代的残差。

　　>> semilogy(0:20,resvec/norm(b),'-o') %作每次迭代的相对残差图形，结果如下图。

　　>> xlabel('iteration number') %x轴为迭代次数。

　　>> ylabel('relative residual') %y轴为相对残差。

　　1.4.6 稳定双共轭梯度方法解方程组

　　函数 bicgstab

　　格式 x =bicgstab(A,b)

　　bicgstab(A,b,tol)

　　bicgstab(A,b,tol,maxit)

　　bicgstab(A,b,tol,maxit,M)

　　bicgstab(A,b,tol,maxit,M1,M2)

　　bicgstab(A,b,tol,maxit,M1,M2,x0)

　　[x,flag] = bicgstab(A,b,…)

　　[x,flag,relres] = bicgstab(A,b,…)

　　[x,flag,relres,iter] = bicgstab(A,b,…)

　　[x,flag,relres,iter,resvec] = bicgstab(A,b,…)

　　稳定双共轭梯度法解方程组，调用方式和返回的结果形式和命令bicg一样。

　　例1-84

　　>>load west0479;

　　>>A=west0479;

　　>>b=sum(A,2);

　　>>[x,flag]=bicgstab(A,b)

　　显示结果是x的值全为0，flag=1。表示在默认误差和默认迭代次数(20次)下不收敛。

　　若改为：

　　>>[L1,U1] = luinc(A,1e-5);

　　>>[x1,flag1] = bicgstab(A,b,1e-6,20,L1,U1)

　　即指定误差，并用A的不完全LU分解因子L和U作为预处理因子M=L*U，其结果是x1的值全为0，flag=2表示预处理因子为坏条件的预处理因子。

　　若改为

　　>>[L2,U2]=luinc(A,1e-6); %稀疏矩阵的不完全LU分解。

　　>>[x2,flag2,relres2,iter2,resvec2]=bicgstab(A,b,1e-15,10,L2,U2)

　　%指定最大迭代次数为10次，预处理因子M=L*U。

　　>>semilogy(0:0.5:iter2,resvec2/norm(b),'-o') %每次迭代的相对残差图形。

　　>>xlabel('iteration number')

　　>>ylabel('relative residual')

　　结果为

　　x2=(其值全为1，略)

　　flag2 =

　　0 %表示收敛

　　relres2 =

　　2.8534e-016 %收敛时的相对误差

　　iter2 =

　　6 %计算终止时的迭代次数

　　resvec2 = %每次迭代的残差

　　1.4.7 复共轭梯度平方法解方程组

　　函数 cgs

　　格式 x = cgs(A,b)

　　cgs(A,b,tol)

　　cgs(A,b,tol,maxit)

　　cgs(A,b,tol,maxit,M)

　　cgs(A,b,tol,maxit,M1,M2)

　　cgs(A,b,tol,maxit,M1,M2,x0)

　　[x,flag] = cgs(A,b,…)

　　[x,flag,relres] = cgs(A,b,…)

　　[x,flag,relres,iter] = cgs(A,b,…)

　　[x,flag,relres,iter,resvec] = cgs(A,b,…)

　　调用方式和返回的结果形式与命令bicg一样。

　　1.4.8 共轭梯度的LSQR方法

　　函数 lsqr

　　格式 x = lsqr(A,b)

　　lsqr(A,b,tol)

　　lsqr(A,b,tol,maxit)

　　lsqr(A,b,tol,maxit,M)

　　lsqr(A,b,tol,maxit,M1,M2)

　　lsqr(A,b,tol,maxit,M1,M2,x0)

　　[x,flag] = lsqr(A,b,…)

　　[x,flag,relres] = lsqr(A,b,…)

　　[x,flag,relres,iter] = lsqr(A,b,…)

　　[x,flag,relres,iter,resvec] = lsqr(A,b,…)

　　调用方式和返回的结果形式与命令bicg一样。

　　例1-85

　　>>n = 100;

　　>>on = ones(n,1);

　　>>A = spdiags([-2*on 4*on -on],-1:1,n,n); %产生一个对角矩阵

　　>>b = sum(A,2);

　　>>tol = 1e-8; %指定精度

　　>>maxit = 15; %指定最大迭代次数

　　>>M1 = spdiags([on/(-2) on],-1:0,n,n);

　　>>M2 = spdiags([4*on -on],0:1,n,n); %M1*M2=M，即产生预处理因子

　　>>[x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,[])

　　结果显示

　　x的值全为1

　　flag =

　　0 %表示收敛

　　relres =

　　3.5241e-009 %表示相对残差

　　iter =

　　12 %计算终止时的迭代次数

　　1.4.9 广义最小残差法

　　函数 qmres

　　格式 x = gmres(A,b)

　　gmres(A,b,restart)

　　gmres(A,b,restart,tol)

　　gmres(A,b,restart,tol,maxit)

　　gmres(A,b,restart,tol,maxit,M)

　　gmres(A,b,restart,tol,maxit,M1,M2)

　　gmres(A,b,restart,tol,maxit,M1,M2,x0)

　　[x,flag] = gmres(A,b,…)

　　[x,flag,relres] = gmres(A,b,…)

　　[x,flag,relres,iter] = gmres(A,b,…)

　　[x,flag,relres,iter,resvec] = gmres(A,b,…)

　　除参数restart(缺省时为系数方阵A的阶数n)可以给出外，调用方式和返回的结果形式与命令bicg一样。

　　例1-86

　　>>load west0479;

　　>>A = west0479;

　　>>b = sum(A,2);

　　>>[x,flag] = gmres(A,b,5)

　　结果显示flag=1，表示在默认精度和默认迭代次数下不收敛。

　　若最后一行改为

　　>>[L1,U1] = luinc(A,1e-5);

　　>>[x1,flag1] = gmres(A,b,5,1e-6,5,L1,U1)

　　结果flag1=2，说明该方法失败，原因是U1的对角线上有0元素，构成坏条件的预处理因子。

　　若改为：

　　[L2,U2] = luinc(A,1e-6);

　　tol = 1e-15;

　　[x4,flag4,relres4,iter4,resvec4] = gmres(A,b,4,tol,5,L2,U2)

　　[x6,flag6,relres6,iter6,resvec6] = gmres(A,b,6,tol,3,L2,U2)

　　[x8,flag8,relres8,iter8,resvec8] = gmres(A,b,8,tol,3,L2,U2)

　　结果flag4=0，flag6=0，flag8=0表示参数restart分别取为4，6，8时收敛。

　　1.4.10 最小残差法解方程组

　　函数 minres

　　格式 x = minres(A,b)

　　minres(A,b,tol)

　　minres(A,b,tol,maxit)

　　minres(A,b,tol.maxit,M)

　　minres(A,b,tol,maxit,M1,M2)

　　minres(A,b,tol,maxit,M1,M2,x0)

　　[x,flag] = minres(A,b,…)

　　[x,flag,relres] = minres(A,b,…)

　　[x,flag,relres,iter] = minres(A,b,…)

　　[x,flag,relres,iter,resvec] = minres(A,b,…)

　　[x,flag,relres,iter,resvec,resveccg] = minres(A,b,…)

　　这里A为对称矩阵，这种方法是寻找最小残差来求x。调用方式和返回的结果形式与命令bicg一样。

　　例1-87

　　>>n = 100; on = ones(n,1);

　　>>A = spdiags([-2*on 4*on -2*on],-1:1,n,n);

　　>>b = sum(A,2);

　　>>tol = 1e-10;

　　>>maxit = 50;

　　>>M1 = spdiags(4*on,0,n,n);

　　>> [x,flag,relres,iter,resvec,resveccg] = minres(A,b,tol,maxit,M1,[],[])

　　结果显示：

　　flag =

　　0

　　relres =

　　4.6537e-014

　　iter =

　　49

　　resvec =(略)

　　resveccg =(略)

独善其身 · 发表于 2016-4-7 16:23

　　1.4.11 预处理共轭梯度方法

　　函数 pcg

　　格式 x = pcg(A,b)

　　pcg(A,b,tol)

　　pcg(A,b,tol,maxit)

　　pcg(A,b,tol,maxit,M)

　　pcg(A,b,tol,maxit,M1,M2)

　　pcg(A,b,tol,maxit,M1,M2,x0)

　　pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)

　　[x,flag] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)

　　[x,flag,relres] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)

　　[x,flag,relres,iter] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)

　　[x,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)

　　调用方式和返回的结果形式与命令bicg一样，这里A为对称正定矩阵。

　　1.4.12 准最小残差法解方程组

　　函数 qmr

　　格式 x = qmr(A,b)

　　qmr(A,b,tol)

　　qmr(A,b,tol,maxit)

　　qmr(A,b,tol,maxit,M)

　　qmr(A,b,tol,maxit,M1,M2)

　　qmr(A,b,tol,maxit,M1,M2,x0)

　　qmr(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,…)

　　[x,flag] = qmr(A,b,…)

　　[x,flag,relres] = qmr(A,b,…)

　　[x,flag,relres,iter] = qmr(A,b,…)

　　[x,flag,relres,iter,resvec] = qmr(A,b,…)

　　调用方式和返回的结果形式与命令bicg一样，这里A为方阵。

　　例1-88

　　>>load west0479;

　　>>A = west0479;

　　>>b = sum(A,2);

　　>>[x,flag] = qmr(A,b)

　　结果显示flag=1，表示在缺省情况下不收敛

　　若改为

　　>>[L1,U1] = luinc(A,1e-5);

　　>>[x1,flag1] = qmr(A,b,1e-6,20,L1,U1)

　　结果显示 flag="2"，表示是坏条件的预处理因子，不收敛。

　　若改为

　　>>[L2,U2] = luinc(A,1e-6);

　　>>[x2,flag2,relres2,iter2,resvec2] = qmr(A,b,1e-15,10,L2,U2)

　　>>semilogy(0:iter2,resvec2/norm(b),'-o') %每次迭代的相对残差图形

　　>>xlabel('iteration number')

　　>>ylabel('relative residual')

　　结果为

　　x=(全为1，略)

　　flag2 =

　　0 %表示收敛

　　relres2 = %表示相对残差

　　2.8715e-016

　　iter2 = %计算终止时的迭代次数

　　8

　　resvec2 = %每次迭代的残差

　　1.0e+005 *

　　7.0557

　　7.1773

　　3.4032

　　1.7220

　　0.0000

　　0.0000

　　0.0000

　　0.0000

　　0.0000

　　1.5 特征值与二次型

　　工程技术中的一些问题，如振动问题和稳定性问题，常归结为求一个方阵的特征值和特征向量。

　　1.5.1 特征值与特征向量的求法

　　设A为n阶方阵，如果数“ ”和n维列向量x使得关系式成立，则称为方阵A的特征值，非零向量x称为A对应于特征值“ ”的特征向量。

　　详见1.3.5和1.3.6节：特征值分解问题。

　　例1-89 求矩阵的特征值和特征向量

　　解：

　　>>A=[-2 1 1;0 2 0;-4 1 3];

　　>>[V,D]=eig(A)

　　结果显示：

　　V =

　　-0.7071 -0.2425 0.3015

　　0 0 0.9045

　　-0.7071 -0.9701 0.3015

　　D =

　　-1 0 0

　　0 2 0

　　0 0 2

　　即：特征值-1对应特征向量(-0.7071 0 -0.7071)T

　　特征值2对应特征向量(-0.2425 0 -0.9701)T和(-0.3015 0.9045 -0.3015)T

　　例1-90 求矩阵的特征值和特征向量。

　　解：

　　>>A=[-1 1 0;-4 3 0;1 0 2];

　　>>[V,D]=eig(A)

　　结果显示为

　　V =

　　0 0.4082 -0.4082

　　0 0.8165 -0.8165

　　1.0000 -0.4082 0.4082

　　D =

　　2 0 0

　　0 1 0

　　0 0 1

　　说明当特征值为1 (二重根)时，对应特征向量都是k (0.4082 0.8165 -0.4082)T，k为任意常数。

　　1.5.2 提高特征值的计算精度

　　函数 balance

　　格式 [T,B] = balance(A) %求相似变换矩阵T和平衡矩阵B，满足。

　　B = balance(A) %求平衡矩阵B

　　1.5.3 复对角矩阵转化为实对角矩阵

　　函数 cdf2rdf

　　格式 [V,D] = cdf2rdf (v,d) %将复对角阵d变为实对角阵D，在对角线上，用2×2实数块代替共轭复数对。

　　例1-91

　　>> A=[1 2 3;0 4 5;0 -5 4];

　　>> [v,d]=eig(A)

　　v =

　　1.0000 -0.0191 - 0.4002i -0.0191 + 0.4002i

　　0 0 - 0.6479i 0 + 0.6479i

　　0 0.6479 0.6479

　　d =

　　1.0000 0 0

　　0 4.0000 + 5.0000i 0

　　0 0 4.0000 - 5.0000i

　　>> [V,D]=cdf2rdf(v,d)

　　V =

　　1.0000 -0.0191 -0.4002

　　0 0 -0.6479

　　0 0.6479 0

　　D =

　　1.0000 0 0

　　0 4.0000 5.0000

　　0 -5.0000 4.0000

　　1.5.4 正交基

　　命令 orth

　　格式 B="orth"(A) %将矩阵A正交规范化，B的列与A的列具有相同的空间，B的列向量是正交向量，且满足：B'*B = eye(rank(A))。

　　例1-92 将矩阵正交规范化。

　　解：

　　>>A=[4 0 0; 0 3 1; 0 1 3];

　　>>B=orth(A)

　　>>Q=B'*B

　　则显示结果为

　　P =

　　1.0000 0 0

　　0 0.7071 -0.7071

　　0 0.7071 0.7071

　　Q =

　　1.0000 0 0

　　0 1.0000 0

　　0 0 1.0000

　　1.5.5 二次型

　　例1-93 求一个正交变换X=PY，把二次型

　　化成标准形。

　　解：先写出二次型的实对称矩阵

　　在Matlab编辑器中建立M文件如下：

　　A=[0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0];

　　[P,D]=schur(A)

　　syms y1 y2 y3 y4

　　y=[y1;y2;y3;y4];

　　X=vpa(P,2)*y %vpa表示可变精度计算，这里取2位精度

　　f=[y1 y2 y3 y4]*D*y

　　运行后结果显示如下：

　　P =

　　780/989 780/3691 1/2 -390/1351

　　780/3691 780/989 -1/2 390/1351

　　780/1351 -780/1351 -1/2 390/1351

　　0 0 1/2 1170/1351

　　D =

　　1 0 0 0

　　0 1 0 0

　　0 0 -3 0

　　0 0 0 1

　　X =

　　[ .79*y1+.21*y2+.50*y3-.29*y4]

　　[ .21*y1+.79*y2-.50*y3+.29*y4]

　　[ .56*y1-.56*y2-.50*y3+.29*y4]

　　[ .50*y3+.85*y4]

　　f =

　　y1^2+y2^2-3*y3^2+y4^2

　　1.6 秩与线性相关性

　　1.6.1 矩阵和向量组的秩以及向量组的线性相关性

　　矩阵A的秩是矩阵A中最高阶非零子式的阶数;向量组的秩通常由该向量组构成的矩阵来计算。

　　函数 rank

　　格式 k = rank(A) %返回矩阵A的行(或列)向量中线性无关个数

　　k = rank(A,tol) %tol为给定误差

　　例1-94 求向量组的秩，并判断其线性相关性。

　　>>A=[1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4];

　　>>k=rank(A)

　　结果为

　　k =

　　3

　　由于秩为3 < 向量个数，因此向量组线性相关。

　　1.6.2 求行阶梯矩阵及向量组的基

　　行阶梯使用初等行变换，矩阵的初等行变换有三条：

　　1.交换两行 (第i、第j两行交换)

　　2.第i行的K倍

　　3.第i行的K倍加到第j行上去

　　通过这三条变换可以将矩阵化成行最简形，从而找出列向量组的一个最大无关组，Matlab将矩阵化成行最简形的命令是rref或rrefmovie。

　　函数 rref或rrefmovie

　　格式 R = rref(A) %用高斯—约当消元法和行主元法求A的行最简行矩阵R

　　[R,jb] = rref(A) %jb是一个向量，其含义为：r = length(jb)为A的秩;A(:, jb)为A的列向量基;jb中元素表示基向量所在的列。

　　[R,jb] = rref(A,tol) %tol为指定的精度

　　rrefmovie(A) %给出每一步化简的过程

　　例1-95 求向量组a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)的一个最大无关组。

　　>> a1=[1 -2 2 3]';

　　>>a2=[-2 4 -1 3]';

　　>>a3=[-1 2 0 3]';

　　>>a4=[0 6 2 3]';

　　>>a5=[2 -6 3 4]';

　　A=[a1 a2 a3 a4 a5]

　　A =

　　1 -2 -1 0 2

　　-2 4 2 6 -6

　　2 -1 0 2 3

　　3 3 3 3 4

　　>> [R,jb]=rref(A)

　　R =

　　1.0000 0 0.3333 0 1.7778

　　0 1.0000 0.6667 0 -0.1111

　　0 0 0 1.0000 -0.3333

　　0 0 0 0 0

　　jb =

　　1 2 4

　　>> A(:,jb)

　　ans =

　　1 -2 0

　　-2 4 6

　　2 -1 2

　　3 3 3

　　即：a1 a2 a4为向量组的一个基。

　　1.7 稀疏矩阵技术

　　1.7.1 稀疏矩阵的创建

　　函数 sparse

　　格式 S = sparse(A) %将矩阵A转化为稀疏矩阵形式，即由A的非零元素和下标构成稀疏矩阵S。若A本身为稀疏矩阵，则返回A本身。

　　S = sparse(m,n) %生成一个m×n的所有元素都是0的稀疏矩阵

　　S = sparse(i,j,s) %生成一个由长度相同的向量i，j和s定义的稀疏矩阵S，其中i，j是整数向量，定义稀疏矩阵的元素位置(i,j)，s是一个标量或与i，j长度相同的向量，表示在(i,j)位置上的元素。

　　S = sparse(i,j,s,m,n) %生成一个m×n的稀疏矩阵，(i,j)对应位置元素为si，m = max(i)且n =max(j)。

　　S = sparse(i,j,s,m,n,nzmax) %生成一个m×n的含有nzmax个非零元素的稀疏矩阵S，nzmax的值必须大于或者等于向量i和j的长度。

　　例1-96

　　>> S="sparse"(1:10,1:10,1:10)

　　S =

　　(1,1) 1

　　(2,2) 2

　　(3,3) 3

　　(4,4) 4

　　(5,5) 5

　　(6,6) 6

　　(7,7) 7

　　(8,8) 8

　　(9,9) 9

　　(10,10) 10

　　>> S="sparse"(1:10,1:10,5)

　　S =

　　(1,1) 5

　　(2,2) 5

　　(3,3) 5

　　(4,4) 5

　　(5,5) 5

　　(6,6) 5

　　(7,7) 5

　　(8,8) 5

　　(9,9) 5

　　(10,10) 5

　　1.7.2 将稀疏矩阵转化为满矩阵

　　函数 full

　　格式 A="full"(S) %S为稀疏矩阵，A为满矩阵。

　　例1-97

　　>> S="sparse"(1:5,1:5,4:8)

　　S =

　　(1,1) 4

　　(2,2) 5

　　(3,3) 6

　　(4,4) 7

　　(5,5) 8

　　>> A="full"(S)

　　A =

　　4 0 0 0 0

　　0 5 0 0 0

　　0 0 6 0 0

　　0 0 0 7 0

　　0 0 0 0 8

　　1.7.3 稀疏矩阵非零元素的索引

　　函数 find

　　格式 k = find(x) %按行检索X中非零元素的点，若没有非零元素，将返回空矩阵。

　　[i,j] = find(X) %检索X中非零元素的行标i和列标j

　　[i,j,v] = find(X) %检索X中非零元素的行标i和列标j以及对应的元素值v

　　例1-98 上例中

　　>> [i,j,v]=find(S)

　　则显示为：

　　I =[ 1 2 3 4 5]’

　　j =[ 1 2 3 4 5]’

　　v =[4 5 6 7 8]’

　　1.7.4 外部数据转化为稀疏矩阵

　　函数 spconvert

　　格式 S="spconvert"(D) %D是只有3列或4列的矩阵

　　说明：先运用load函数把外部数据(.mat文件或.dat文件)装载于MATLAB内存空间中的变量T;T数组的行维为nnz或nnz+1，列维为3(对实数而言)或列维为4(对复数而言);T数组的每一行(以[i,j,Sre,Sim]形式)指定一个稀疏矩阵元素。

　　例1-99

　　>> D=[1 2 3;2 5 4;3 4 6;3 6 7]

　　D =

　　1 2 3

　　2 5 4

　　3 4 6

　　3 6 7

　　>> S="spconvert"(D)

　　S =

　　(1,2) 3

　　(3,4) 6

　　(2,5) 4

　　(3,6) 7

　　>> D=[1 2 3 4; 2 5 4 0;3 4 6 9;3 6 7 4];

　　D =

　　1 2 3 4

　　2 5 4 0

　　3 4 6 9

　　3 6 7 4

　　>> S="spconvert"(D)

　　S =

　　(1,2) 3.0000 + 4.0000i

　　(3,4) 6.0000 + 9.0000i

　　(2,5) 4.0000

　　(3,6) 7.0000 + 4.0000i

　　1.7.5 基本稀疏矩阵

　　1.带状(对角)稀疏矩阵

　　函数 spdiags

　　格式 [B,d] = spdiags(A) %从矩阵A中提取所有非零对角元素，这些元素保存在矩阵B中，向量d表示非零元素的对角线位置。

　　B = spdiags(A,d) %从A中提取由d指定的对角线元素，并存放在B中。

　　A = spdiags(B,d,A) %用B中的列替换A中由d指定的对角线元素，输出稀疏矩阵。

　　A = spdiags(B,d,m,n) %产生一个m×n稀疏矩阵A，其元素是B中的列元素放

　　在由d指定的对角线位置上。

　　例1-100

　　>>A = [11 0 13 0

　　0 22 0 24

　　0 0 33 0

　　41 0 0 44

　　0 52 0 0

　　0 0 63 0

　　0 0 0 74];

　　>>[B,d] = spdiags(A)

　　B =

　　41 11 0

　　52 22 0

　　63 33 13

　　74 44 24

　　d =

　　-3 %表示B的第1列元素在A中主对角线下方第3条对角线上

　　0 %表示B的第2列在A的主对角线上

　　2 %表示B的第3列在A的主对角线上方第2条对角线上

　　例1-101

　　>> B=[1 2 3 4

　　5 6 7 8

　　9 10 11 12

　　13 14 15 16];

　　>> d=[-2 0 1 3];

　　>> A="spdiags"(B,d,4,4);

　　>> full(A)

　　ans =

　　2 7 0 16

　　0 6 11 0

　　1 0 10 15

　　0 5 0 14

　　2.单位稀疏矩阵

　　函数 speye

　　格式 S = speye(m,n) %生成m×n的单位稀疏矩阵

　　S = speye(n) %生成n×n的单位稀疏矩阵

　　3.稀疏均匀分布随机矩阵

　　函数 sprand

　　格式 R = sprand(S) %生成与S具有相同稀疏结构的均匀分布随机矩阵

　　R = sprand(m,n,density) %生成一个m×n的服从均匀分布的随机稀疏矩阵，非零元素的分布密度是density。

　　R = sprand(m,n,density,rc) %生成一个近似的条件数为1/rc、大小为m×n的均匀分布的随机稀疏矩阵。

　　4.稀疏正态分布随机矩阵

　　函数 sprandn

　　格式 R = sprandn(S) %生成与S具有相同稀疏结构的正态分布随机矩阵。

　　R = sprandn(m,n,density) %生成一个m×n的服从正态分布的随机稀疏矩阵，非零元素的分布密度是density。

　　R = sprandn(m,n,density,rc) %生成一个近似的条件数为1/rc、大小为m×n的均匀分布的随机稀疏矩阵。

　　5.稀疏对称随机矩阵

　　函数 sprandsym

　　格式 R = sprandsym(S) %生成稀疏对称随机矩阵，其下三角和对角线与S具有相同的结构，其元素服从均值为0、方差为1的标准正态分布。

　　R = sprandsym(n,density) %生成n×n的稀疏对称随机矩阵，矩阵元素服从正态分布，分布密度为density。

　　R = sprandsym(n,density,rc) %生成近似条件数为1/rc的稀疏对称随机矩阵

　　R = sprandsym(n,density,rc,kind) %生成一个正定矩阵，参数kind取值为kind=1表示矩阵由一正定对角矩阵经随机Jacobi旋转得到，其条件数正好为1/rc;kind=2表示矩阵为外积的换位和，其条件数近似等于1/rc;kind=3表示生成一个与矩阵S结构相同的稀疏随机矩阵，条件数近似为1/rc ，density被忽略。

　　1.7.6 稀疏矩阵的运算

　　1.稀疏矩阵非零元素的个数

　　函数 nnz

　　格式 n = nnz(X) %返回矩阵X中非零元素的个数

　　2.稀疏矩阵的非零元素

　　函数 nonzeros

　　格式 s = nonzeros(A) %返回矩阵A中非零元素按列顺序构成的列向量

　　例1-102

　　>>A =[ 2 7 0 16

　　0 6 11 0

　　1 0 10 15

　　0 5 0 14];

　　>> s="nonzeros"(A)

　　结果为：

　　s =[ 2 1 7 6 5 11 10 16 15 14]’

　　3.稀疏矩阵非零元素的内存分配

　　函数 nzmax

　　格式 n = nzmax(S) %返回非零元素分配的内存总数n

　　4.稀疏矩阵的存贮空间

　　函数 spalloc

　　格式 S = spalloc(m,n,nzmax) %产生一个m×n阶只有nzmax个非零元素的稀疏矩阵，这样可以有效减少存贮空间和提高运算速度。

　　5.稀疏矩阵的非零元素应用

　　函数 spfun

　　格式 f = spfun('function',S) %用S中非零元素对函数'function'求值，如果'function'不是对稀疏矩阵定义的，同样可以求值。

　　例1-103 4阶稀疏矩阵对角矩阵

　　S =

　　(1,1) 1

　　(2,2) 2

　　(3,3) 3

　　(4,4) 4

　　f= spfun('exp',S) %即指数e的非零元素方幂。

　　结果为

　　f =

　　(1,1) 2.7183

　　(2,2) 7.3891

　　(3,3) 20.0855

　　(4,4) 54.5982

　　6.把稀疏矩阵的非零元素全换为1

　　函数 spones

　　格式 R = spones(S) %将稀疏矩阵S中的非零元素全换为1

　　1.7.7 画稀疏矩阵非零元素的分布图形

　　函数 spy

　　格式 spy(S) %画出稀疏矩阵S中非零元素的分布图形。S也可以是满矩阵。

　　spy(S,markersize) % markersize为整数，指定点阵大小。

　　spy(S,'LineSpec') %'LineSpec'指定绘图标记和颜色

　　spy(S,'LineSpec',markersize) %参数与上面相同

　　例1-104

　　>> load west0479

　　>> A="west0479";

　　>> spy(A,'ro',3)

　　1.7.8 矩阵变换

　　1.列近似最小度排序

　　函数 colamd

　　格式 p = colamd(S) %返回稀疏矩阵S的列的近似最小度排序向量p

　　2.列最小度排序

　　函数 colmmd

　　格式 p = colmmd(S) %返回稀疏矩阵S的列的最小度排序向量p，按p排列后的矩阵为S(: , p)。

　　例1-105 比较稀疏矩阵S与排序后的矩阵S(: , p)

　　>> load west0479;

　　>> S="west0479";

　　>> p="colmmd"(S);

　　>> subplot(2,2,1),spy(S)

　　>> subplot(2,2,2),spy(S(:,p))

　　>> subplot(2,2,3),spy(lu(S))

　　>> subplot(2,2,4),spy(lu(S(:,p)))

　　3.非零元素的列变换

　　函数 colperm

　　格式 j = colperm(S) %返回一个稀疏矩阵S的列变换的向量。列按非0元素升序排列。有时这是LU分解前有用的变换：LU(S(:,j))。如果S是一个对称矩阵，对行和列进行排序，有利于Cholesky分解：chol(S(j,j))

　　4.Dulmage-Mendelsohn分解

　　函数 dmperm

　　格式 p = dmperm (A) %返回A的行排列向量p，这样，如果A满列秩，就使得A(p,:)是具有非0对角线元素的方阵。

　　[p,q,r] = dmperm(A) %A为方阵，p为行排列向量，q为列排列向量，使得A(p,q)

　　是上三角块形式，r为索引向量。

　　[p,q,r,s] = dmperm(A) %A不是方阵，p,q,r含义与前面相同，s也是索引向量。

　　例1-106

　　>> A=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0]

　　A =

　　11 0 13 0

　　41 0 0 44

　　0 22 0 24

　　0 0 63 0

　　>> [p,q,r]=dmperm(A)

　　p =

　　3 2 1 4

　　q =

　　2 4 1 3

　　r =

　　1 2 3 4 5

　　>> A(p,q)

　　ans =

　　22 24 0 0

　　0 44 41 0

　　0 0 11 13

　　0 0 0 63

　　5.整数的随机排列

　　函数 randperm

　　格式 p = randperm (n) %对正整数1，2，3，…，n的随机排列，可以用来创建随机变换矩阵。

　　例1-107

　　>> p="randperm"(6)

　　p =

　　3 4 6 5 1 2

　　6.稀疏对称近似最小度排列

　　函数 symamd

　　格式 p = symamd(S) %S为对称正定矩阵，返回排列向量p。

　　7.稀疏逆Cuthill-McKee排序

　　函数 symrcm

　　格式 r = symrcm (S) %返回S的对称逆Cuthill-McKee排序r，使S的非0元素集中在主对角线附近。

　　8.稀疏对称最小度排列

　　函数 symmmd

　　格式 p = symmmd(S) %返回S的对称最小度排列向量p，S为对称正定矩阵。

　　例1-108

　　>> B = bucky+4*speye(60);

　　>>r = symrcm(B);

　　>>p = symmmd(B);

　　>>R = B(r,r);

　　>>S = B(p,p);

　　>>subplot(2,2,1), spy(R), title('B(r,r)')

　　>>subplot(2,2,2), spy(S), title('B(s,s)')

　　>>subplot(2,2,3), spy(chol(R)), title('chol(B(r,r))')

　　>>subplot(2,2,4), spy(chol(S)), title('chol(B(s,s))')

　　1.7.9 稀疏矩阵的近似欧几里得范数和条件数

　　命令矩阵A条件数的1-范数估计值

　　函数 condest

　　格式 c = condest(A) %计算方阵A的1-范数中条件数的下界值c

　　[c,v] = condest(A) %方阵A的1-范数中条件数的下界值c和向量v，使得，即：

　　norm(A*v,1) = norm(A,1)*norm(v,1)/c。

　　命令 2-范数估计值

　　函数 normest

　　格式 nrm = normest(S) %返回矩阵S的2-范数的估计值，相对误差为10-6。

　　nrm = normest(S,tol) %tol为指定的相对误差，而不用默认误差10-6。

　　[nrm,count] = normest(…) %count为给出的计算范数迭代的次数

　　其条件数的计算与前面1.2.12相同。

　　1.7.10 稀疏矩阵的分解

　　命令不完全的LU分解

　　函数 luinc

　　格式 [L,U] = luinc(X,'0') %X为稀疏方阵;L为下三角矩阵的置换形式;U为上三角矩阵;'0'是一种分解标准。

　　[L,U,P] = luinc(X,'0') %L为下三角阵，其主对角线元素为1;U为上三角阵，p为单位矩阵的置换形式。

　　[L,U] = luinc(X,options) %options取值为：droptol表示指定的舍入误差;milu表示改变分解以便于从上三角分解因子中抽取被去掉的列元素。ugiag为1表示用droptol值代替上三角因子中的对角线上的零元素，默认值为0。thresh为中心临界值。

　　[L,U] = luinc(X,droptol) %droptol表示指定不完全分解的舍入误差

　　[L,U,P] = luinc(X,options)

　　[L,U,P] = luinc(X,droptol)

　　例1-109

　　>> S=[11 0 13 0;41 0 0 44;0 22 0 24;0 0 63 0]

　　S =

　　11 0 13 0

　　41 0 0 44

　　0 22 0 24

　　0 0 63 0

　　>> S="sparse"(S)

　　S =

　　(1,1) 11

　　(2,1) 41

　　(3,2) 22

　　(1,3) 13

　　(4,3) 63

　　(2,4) 44

　　(3,4) 24

　　>> luinc(S,'0')

　　ans =

　　(1,1) 41.0000

　　(4,1) 0.2683

　　(2,2) 22.0000

　　(3,3) 63.0000

　　(4,3) 0.2063

　　(1,4) 44.0000

　　(2,4) 24.0000

　　>> [L,U,p]=luinc(S,'0')

　　L =

　　(1,1) 1.0000

　　(4,1) 0.2683

　　(2,2) 1.0000

　　(3,3) 1.0000

　　(4,3) 0.2063

　　(4,4) 1.0000

　　U =

　　(1,1) 41

　　(2,2) 22

　　(3,3) 63

　　(1,4) 44

　　(2,4) 24

　　p =

　　(4,1) 1

　　(1,2) 1

　　(2,3) 1

　　(3,4) 1

　　命令稀疏矩阵的不完全Cholesky分解

　　函数 cholinc

　　格式 R = (X,droptol) %稀疏矩阵X的不完全Cholesky分解，droptol为指定误差。

　　R = cholinc(X,options) %options取值为：droptol表示舍入误差;michol表示如果michol=1，则从对角线上抽取出被去掉的元素。rdiag表示用droptol值代替上三角分解因子中的对角线上的零元素。

　　R = cholinc(X,'0') %'0'是一种分解标准

　　[R,p] = cholinc(X,'0') %不产生任何出错信息，如果R存在，则p=0;如果R不存在，则p为正整数。

　　R = cholinc(X,'inf') %进行Cholesky无穷分解

　　1.7.11 稀疏矩阵的特征值分解

　　函数 eigs

　　格式 d = eigs(A) %求稀疏矩阵A的6个最大特征值d，d以向量形式存放。

　　d = eigs(A,B) %求稀疏矩阵的广义特征值问题。满足AV=BVD，其中D为特征值对角阵，V为特征向量矩阵，B必须是对称正定阵或Hermitian正定阵。

　　d = eigs(A,k) %返回k个最大特征值

　　d = eigs(A,B,k) %返回k个最大特征值

　　d = eigs(A,k,sigma) %sigma取值：'lm' 表示最大数量的特征值;'sm' 最小数量特征值;对实对称问题：'la'表示最大特征值;'sa'为最小特征值;对非对称和复数问题：'lr' 表示最大实部;'sr' 表示最小实部;'li' 表示最大虚部;'si'表示最小虚部。

　　d = eigs(A,B,k,sigma) %同上

　　d = eigs(A,k,sigma,options) % options为指定参数：参见eigs帮助文件。

　　d = eigs(A,B,k,sigma,options) %同上。以下的参数k、sigma、options相同。

　　d = eigs(Afun,n) %用函数Afun代替A，n为A的阶数，D为特征值。

　　d = eigs(Afun,n,B)

　　d = eigs(Afun,n,k)

　　d = eigs(Afun,n,B,k)

　　d = eigs(Afun,n,k,sigma)

　　d = eigs(Afun,n,B,k,sigma)

　　d = eigs(Afun,n,k,sigma,options)

　　d = eigs(Afun,n,B,k,sigma,options)

　　[V,D] = eigs(A,…) %D为6个最大特征值对角阵，V的列向量为对应特征向量。

　　[V,D] = eigs(Afun,n,…)

　　[V,D,flag] = eigs(A,…) %flag表示特征值的收敛性，若flag=0，则所有特征值都收敛，否则，不是所有都收敛。

　　[V,D,flag] = eigs(Afun,n,…)

　　1.7.12 稀疏矩阵的线性方程组

　　参见1.4节的方程组求解。

账号		自动登录	找回密码
密码			我要加入

[综合讨论] matlab相关矩阵操作及相关实例

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

相关帖子