第五章 求解线性方程组
? 矩阵
? 直接解法
? 迭代法
? 符号解法
? 稀疏矩阵技术
? 特征值与特征向量
5.1 矩阵
5.1.1特殊矩阵的输入
? 数值矩阵的输入
– 零矩阵、幺矩阵及单位矩阵
生成 n?n方阵,
A=zeros(n),B=ones(n),C=eye(n)
生成 m?n矩阵,
A=zeros(m,n),B=ones(m,n),C=eye(m,n)
生成和矩阵 B同样位数的矩阵,
A=zeros(size(B))
– 随机元素矩阵
若矩阵随机元素满足 [0,1]区间上的均匀分布
生成 n?m阶标准均匀分布为随机数矩阵,
A=rand(n,m)
生成 n?n阶标准均匀分布为随机数方阵,
A=rand(n)
– 对角元素矩阵
已知向量生成对角矩阵,
A=diag(V)
已知矩阵提取对角元素列向量,
V＝ diag(A)
生成主对角线上第 k条对角线为 V的矩阵,
A=diag(V,k)
? 例,diag( )函数的不同调用格式
>> C=[1 2 3]; V=diag(C) % 生成对角矩阵
V =
1 0 0
0 2 0
0 0 3
>> V1=diag(V)' % 将列向量通过转置变换成行向量
V1 =
1 2 3
>> C=[1 2 3]; V=diag(C,2) % 主对角线上第 k条对角线为 C的矩阵
V =
0 0 1 0 0
0 0 0 2 0
0 0 0 0 3
0 0 0 0 0
0 0 0 0 0
生成三对角矩阵,
>> V=diag([1 2 3 4])+diag([2 3 4],1)+diag([5 4 3],-1)
V =
1 2 0 0
5 2 3 0
0 4 3 4
0 0 3 4
– Hilbert矩阵及逆 Hilbert矩阵
生成 n阶的 Hilbert矩阵,
A=hilb(n)
求取逆 Hilbert矩阵,
B=invhilb(n)
? 符号矩阵的输入
数值矩阵 A转换成符号矩阵,
B=sys(A)
例,
>> A=hilb(3)
A =
1.0000 0.5000 0.3333
0.5000 0.3333 0.2500
0.3333 0.2500 0.2000
>> B=sym(A)
B =
[ 1,1/2,1/3]
[ 1/2,1/3,1/4]
[ 1/3,1/4,1/5]
5.1.2 矩阵基本概念与性质
? 行列式
格式, d=det(A)
例：求行列式
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1]; det(A)
ans =
0
? 例,
>> tic,A=sym(hilb(20)); det(A),toc
ans =
1/23774547167685345090916442434276164401754
1983775348649303318533123441975931064458
5187585766816573773440565759867265558971
7656384197107933033865823241498112410235
5448916615471780963525779783680000000000
0000000000000000000000000
elapsed_time =
2.3140
高阶的 Hilbert矩阵是接近奇异的矩阵。
? 矩阵的迹
格式,t=trace(A)
? 矩阵的秩
格式,r=rank(A) ％用默认的精度求数值秩
r=rank(A,) ％给定精度下求数值秩 ?
? 例
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1]; rank(A)
ans =
3
该矩阵的秩为 3，小于矩阵的阶次，故为非满秩矩阵。
? 例
>> H=hilb(20); rank(H) ％数值方法
ans =
13
>> H=sym(hilb(20)); rank(H) % 解析方法，原矩阵为非奇异矩阵
ans =
20
? 矩阵范数
? 矩阵的范数定义,
格式,
N=norm(A) ％求解默认的 2范数
N=norm(A，选项 ) ％选项可为 1,2,inf等
? 例：求一向量、矩阵的范数
>> a=[16 2 3 13];
>> [norm(a),norm(a,2),norm(a,1),norm(a,Inf)]
ans =
2.092844953645635e+001 2.092844953645635e+001
3.400000000000000e+001 1.600000000000000e+001
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> [norm(A),norm(A,2),norm(A,1),norm(A,Inf)]
ans =
34 34 34 34
符号运算工具箱未提供 norm( )函数，需先用 double( )
函数转换成双精度数值矩阵，再调用 norm( )函数。
? 特征多项式
格式,C=poly(A)
例,>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> poly(A) ％直接求取
ans =
1.000000000000000e+000 -3.399999999999999e+001
-7.999999999999986e+001 2.719999999999999e+003
-2.819840539024018e-012
>> A=sym(A); poly(A) ％运用符号工具箱
ans =
x^4-34*x^3-80*x^2+2720*x
? 符号多项式与数值多项式的转换
格式,
f=poly2sym(P) 或 f=poly2sym(P,x)
格式,P=sym2poly(f)
? 例,
>> P=[1 2 3 4 5 6]; % 先由系数按降幂顺序排列表示多
项式
>> f=poly2sym(P,'v') % 以 v 为算子表示多项式
f =
v^5+2*v^4+3*v^3+4*v^2+5*v+6
>> P=sym2poly(f)
P =
1 2 3 4 5 6
? 矩阵的逆矩阵
格式,C=inv(A)
例,
>> format long; H=hilb(4); H1=inv(H)
H1 =
1.0e+003 *
0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999
-0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984
0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961
-0.13999999999999 1.67999999999984 -4.19999999999961 2.79999999999974
检验,
>> H*H1
ans =
1.00000000000001 0.00000000000023 -0.00000000000045 0.00000000000023
0.00000000000001 1.00000000000011 -0.00000000000011 0.00000000000011
0.00000000000001 0 1.00000000000011 0
0.00000000000000 0.00000000000011 -0.00000000000011 1.00000000000011
计算误差范数,
>> norm(H*inv(H)-eye(size(H)))
ans =
6.235798190375727e-013
>> H2=invhilb(4); norm(H*H2-eye(size(H)))
ans =
5.684341886080802e-014
>> H=hilb(10); H1=inv(H); norm(H*H1-eye(size(H)))
ans =
0.00264500826202
>> H2=invhilb(10); norm(H*H2-eye(size(H)))
ans =
1.612897415528547e-005
>> H=hilb(13); H1=inv(H); norm(H*H1-eye(size(H)))
Warning,Matrix is close to singular or badly scaled,
Results may be inaccurate,RCOND = 2.339949e-018,
ans =
53.23696008570294
>> H2=invhilb(13); norm(H*H2-eye(size(H)))
ans =
11.37062973181391
对接近于奇异矩阵，高阶一般不建议用 inv( )，可用符号工具箱。
>> H=sym(hilb(7)); inv(H)
ans =
[ 49,-1176,8820,-29400,48510,-38808,12012]
[-1176,37632,-317520,1128960,-1940400,1596672,-504504]
[8820,-317520,2857680,-10584000,18711000,-15717240,5045040]
[-29400,1128960,-10584000,40320000,-72765000,62092800,-20180160]
[48510,-1940400,18711000,-72765000,133402500,-115259760,37837800]
[-38808,1596672,-15717240,62092800,-115259760,100590336,-33297264]
[12012,-504504,5045040,-20180160,37837800,-33297264,11099088]
>> H=sym(hilb(30)); norm(double(H*inv(H)-eye(size(H))))
ans =
0
? 例：奇异阵求逆
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> format long; B = inv(A)
Warning,Matrix is close to singular or badly scaled,
Results may be inaccurate,RCOND = 1.306145e-017,
B =
1.0e+014 *
0.93824992236885 2.81474976710656 -2.81474976710656 -0.93824992236885
2.81474976710656 8.44424930131968 -8.44424930131968 -2.81474976710656
-2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656
-0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885
>> norm(A*B-eye(size(A))) ％检验
ans =
1.64081513306419
>> A=sym(A); inv(A) ％奇异矩阵不存在一个相应的逆矩阵，用符号工具
箱的函数也不行
Error using ==> sym/inv
Error,(in inverse) singular matrix
5.2 线性方程组直接解法
5.2.1线性方程组直接求解－矩阵除法
? 关于线性方程组的直接解法，如 Gauss消去法、
选主元消去法、平方根法、追赶法等等，在
MATLAB中，只需用“／”或,\”就解决问
题。它内部实际包含着许许多多的自适应算
法，如对超定方程用最小二乘法，对欠定方
程时它将给出范数最小的一个解，解三对角
阵方程组时用追赶法等等。
格式,
x=A\b
? 例：解方程组
>> A=[.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129;
.3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927];
>> b=[0.4043 0.1550 0.4240 -0.2557]';
>> x=A\b; x'
ans =
-0.1819 -1.6630 2.2172 -0.4467
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
0,4 0 9 6 0,1 2 3 4 0,3 6 7 8 0,2 9 4 3 0,4 0 4 3
0,2 2 4 6 0,3 8 7 2 0,4 0 1 5 0,1 1 2 9 0,1 1 5 0
0,3 6 4 5 0,1 9 2 0 0,3 7 8 1 0,0 6 4 3 0,4 2 4 0
0,1 7 8 4 0,4 0 0 2 0,2 7 8 6 0,3 9 2 7 0,2 5 5 7
x x x x
x x x x
x x x x
x x x x
? ? ? ??
?
? ? ? ??
?
? ? ? ?
?
? ? ? ? ? ?
?
5.2.2线性方程组直接求解－判定求解
? 例,
>> A=[1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2]; B=[5 1; 4 2; 3 3; 2 4];
>> C=[A B]; rank(A),rank(C)
ans =
4
ans =
4
>> x=inv(A)*B
x =
-1.8000 2.4000
1.8667 -1.2667
3.8667 -3.2667
-2.1333 2.7333
检验
>> norm(A*x-B)
ans =
7.4738e-015
精确解
>> x1=inv(sym(A))*B
x1 =
[ -9/5,12/5]
[ 28/15,-19/15]
[ 58/15,-49/15]
[ -32/15,41/15]
检验
>> norm(double(A*x1-B))
ans =
0
原方程组对应的齐次方程组的解
求取 A矩阵的化零矩阵,
格式,Z=null(A)
求取 A矩阵的化零矩阵的规范形式,
格式,Z=null(A,‘ r ’)
? 例,
判断可解性
>> A=[1 2 3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2]; B=[1;3;2;6];
>> C=[A B]; [rank(A),rank(C)]
ans =
2 2
>> Z=null(A,'r') % 解出规范化的化零空间
Z =
2.0000 3.0000
-2.5000 -3.5000 ％
1.0000 0
0 1.0000
1
23
34
4
23
2.5 3.5
10
01
x
xx
xx
x
?? ??
?? ??
?? ??
?? ???
??
?? ?? ??
?? ??
?? ????
>> x0=pinv(A)*B % 得出一个特解
x0 =
0.9542
0.7328 %全部解
-0.0763
-0.2977
验证得出的解
>> a1=randn(1); a2=rand(1); % 取不同分布的随
机数
>> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(A*x-B)
ans =
4.4409e-015
1
2
12
3
4
2 3 0.9542
2.5 3.5 0.7328
1 0 - 0.0763
0 1 - 0.2977
x
x
aa
x
x
?? ? ? ? ? ? ?
?? ? ? ? ? ? ?
??
?? ? ? ? ? ? ?? ? ?
?? ? ? ? ? ? ?
?? ? ? ? ? ? ?
?? ? ? ? ? ? ???
? 解析解
>> Z=null(sym(A))
Z =
[ 2,3]
[ -5/2,-7/2]
[ 1,0]
[ 0,1]
>> x0=sym(pinv(A)*B)
x0 =
[ 125/131]
[ 96/131]
[ -10/131] ％
[ -39/131]
? 验证得出的解
>> a1=randn(1); a2=rand(1); % 取不同分布的随机数
>> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(double(A*x-B))
ans =
0
? 通解
>> syms a1 a2;
>> x=a1*Z(:,1)+a2*Z(:,2)+x0
x =
[ 2*a1+3*a2+125/131]
[ -5/2*a1-7/2*a2+96/131]
[ a1-10/131]
[ a2-39/131]
摩尔－彭罗斯广义逆求解出的方程最小二乘解
不满足原始代数方程。
5.2.3 线性方程组的直接求解分析
? LU分解
? 格式
[l,u,p]=lu(A)
L是一个单位下三角矩阵,u是一个上三角矩
阵,
p是代表选主元的置换矩阵。
故,Ax=y => PAx=Py
=> LUx=Py => PA=LU
[l,u]=lu(A)
其中 l等于 P-1 L,u等于 U，所以 (P-1 L)U=A
? 例：对 A进行 LU分解
>> A=[1 2 3; 2 4 1; 4 6 7];
>> [l,u,p]=lu(A)
l =
1.0000 0 0
0.5000 1.0000 0
0.2500 0.5000 1.0000
u =
4.0000 6.0000 7.0000
0 1.0000 -2.5000
0 0 2.5000
p =
0 0 1
0 1 0
1 0 0
1 2 3
2 4 1
467
A
??
???
??
??
>> [l,u]=lu(A) ％ l＝ P-1 L
l =
0.2500 0.5000 1.0000
0.5000 1.0000 0
1.0000 0 0
u =
4.0000 6.0000 7.0000
0 1.0000 -2.5000
0 0 2.5000
? Cholesky分解
若矩阵 A为 n阶对称正定阵,则存在唯
一的对角元素为正的三角阵 D，使得
格式,
D=chol(A)
例：进行 Cholesky分解。
>> A=[16 4 8; 4 5 -4; 8 -4 22];
>> D=chol(A)
D =
4 1 2
0 2 -3
0 0 3
1 6 4 8
4 5 4
8 4 2 2
A
??
??
??
??
???
? 三个变换
在线性方程组的迭代求解中，要用到系数
矩阵 A的上三角矩阵、对角阵和下三角矩阵。
此三个变换在 MATLAB中可由以下函数实现。
– 上三角变换,
格式 triu(A,1)
– 对角变换,
格式 diag(A)
– 下三角变换,
格式 tril(A,-1)
例：对此矩阵做三种变换。
1 2 2
1 1 1
2 2 1
A
???
??
? ??
??
>> A=[1 2 -2;1 1 1;2 2 1]; ％
>> triu(A,1)
ans =
0 2 -2
0 0 1
0 0 0
>> tril(A,-1)
ans =
0 0 0
1 0 0
2 2 0
>> b=diag(A); b'
ans =
1 1 1
1 2 2
1 1 1
2 2 1
A
???
??
? ??
??
5.3 迭代解法的几种形式
5.3.1 Jacobi迭代法
? 方程组 Ax=b
A可写成 A=D-L-U
其中,D=diag[a11,a22,…,a nn],-L,-U分别为 A的
严格下、上三角部分（不包括对角线元素）,
由 Ax=b ＝ > x=Bx+f
由此可构造迭代法,
x(k+1)=Bx(k)+f
其中,B=D-1(L+U)=I-D-1A,f=D-1b,
function y=jacobi(a,b,x0)
D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1);
B=D\(L+U); f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=B*x0+f;
n=n+1;
end
n
? 例：用 Jacobi方法求解,
设 x(0)=0,精度为 10-6。
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> jacobi(a,b,[0;0;0])
n =
11
ans =
0.9958
0.9579
0.7916
12
1 2 3
23
10 9
10 2 7
2 10 6
xx
x x x
xx
???
?
? ? ? ??
? ? ? ?
?
5.3.2 Gauss-Seidel迭代法
由原方程构造迭代方程
x(k+1)=G x(k)+f
其中,G=(D-L)-1 U,f=(D-L)-1 b
D=diag[a11,a22,…,a nn],
-L,-U分别为 A的严格下、上三角部
分（不包括对角线元素）,
function y=seidel(a,b,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
G=(D-L)\U ;f=(D-L)\b;
y=G*x0+f; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=G*x0+f;
n=n+1;
end
n
? 例：对上例用 Gauss-Seidel迭代法求解
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> seidel(a,b,[0;0;0])
n =
7
ans =
0.9958
0.9579
0.7916
例：分别用 Jacobi和 G-S
法迭代求解,看是否收敛。
1
2
3
1 2 2 9
1 1 1 7
2 2 1 6
x
x
x
?? ? ? ? ? ?
? ? ? ? ? ?
?? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
>> a=[1 2 -2; 1 1 1; 2 2 1]; b=[9; 7; 6];
>> jacobi(a,b,[0;0;0])
n =
4
ans =
-27
26
8
>> seidel(a,b,[0;0;0])
n =
1011
ans =
1.0e+305 *
-Inf
Inf
-1.7556
5.3.3 SOR迭代法
在很多情况下,J法和 G-S法收敛较
慢，所以考虑对 G-S法进行改进。于是引
入一种新的迭代法－逐次超松弛迭代法
（ Succesise Over-Relaxation),记为 SQR法。
迭代公式为,
X(k+1)= (D-wL)-1((1-w)D+wU)x(k)
+ w(D-wL)-1 b
其中,w最佳值在 [1,2）之间，不易计
算得到，因此 w通常有经验给出。
function y=sor(a,b,w,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
M=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w;
y=M*x0+f; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=M*x0+f;
n=n+1;
end
n
例：上例中，当 w=1.103时，用 SOR法求解
原方程。
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> sor(a,b,1.103,[0;0;0])
n =
8
ans =
0.9958
0.9579
0.7916
5.3.4 两步迭代法
当线性方程系数矩阵为对称正定时，
可用一种特殊的迭代法来解决，其迭代
公式为,
（ D-L） x(k+1/2) =U x(k) +b
（ D-U） x(k+1)=Lx(k+1/2) +b
=>
x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 b
x(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 b
function y=twostp(a,b,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
G1==(D-L)\U; f1=(D-L)\b;
G2==(D-U)\L; f1=(D-U)\b;
y=G1*x0+f1; y=G2*y+f2; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=G1*x0+f1; y=G2*y+f2;
n=n+1;
end
n
? 例：求解方程组
>> a=[10 -1 2 0; -1 11 -1 3; 2 -1 10 3; 0 3 -1 8];
>> b=[6; 25; -11; 15];
>> twostp(a,b,[0; 0; 0; 0])
n =
7
ans =
1.0791
1.9824
-1.4044
0.9560
1
2
3
4
10 1 2 0 6
1 11 1 3 25
2 1 10 1 11
0 3 1 8 15
x
x
x
x
? ??? ? ? ?
??? ? ? ?
?? ??? ? ? ?
?
??? ? ? ?? ? ?
??? ? ? ?? ? ? ?
?? ? ? ???
5.4 线性方程组的符号解法
在 MATLAB的 Symbolic Toolbox中提供了
线性方程的符号求解函数，如
linsolve(A,b)
等同于 X = sym(A)\sym(b),
solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN ')
? 例,
>> A=sym('[10,-1,0;-1,10,-2;0,-2,10]')；
>> b=('[9; 7; 6]')；
>> linsolve(A,b)
ans =
[ 473/475]
[ 91/95]
[ 376/475]
>> vpa(ans)
ans =
[,99578947368421052631578947368421]
[,95789473684210526315789473684211]
[,79157894736842105263157894736842]
例,
>> [x,y] = solve('x^2 + x*y + y = 3','x^2
- 4*x + 3 = 0','x','y')
x =
[ 1]
[ 3]
y =
[ 1]
[ -3/2]
5.5 稀疏矩阵技术
? 稀疏矩阵的建立,
– 格式 S=sparse(i,j,s,m,n)
生成一 mxn阶的稀疏矩阵，以向量 i和 j为坐标的位置上对应
元素值为 s。
– 例,
>> n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n)
a1 =
(1,1) 4
(2,2) 4
(3,3) 4
(4,4) 4
(5,5) 4
– 例,
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n)
a2 =
(2,1) 1
(3,2) 1
(4,3) 1
(5,4) 1
>> full(a2)
– ans =
0 0 0 0 0
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
? 例,n=5,建立主对角线上元素为 4，两条次对角
线为 1的三对角阵。
>> n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n);
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n);
>> a=a1+a2+a2'
a =
(1,1) 4
(2,1) 1
(1,2) 1
(2,2) 4
(3,2) 1
(2,3) 1
(3,3) 4
(4,3) 1
(3,4) 1
(4,4) 4
(5,4) 1
(4,5) 1
(5,5) 4
>> full(a)
ans =
4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4
– 格式 A=spdiags(B,d,m,n)
生成一 mxn阶的稀疏矩阵，使得 B的列放在由 d指定的
位置。
– 例,
>> n=5
>> b=spdiags([ones(n,1),4*ones(n,1),ones(n,1)],…
[-1,0,1],n,n);
>> full(b)
ans =
4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4
– 格式,spconvert(dd)
对于无规律的稀疏矩阵，可使用此命令
由外部数据转化为稀疏矩阵。
– 调用形式为：先用 load函数加载以行表示对
应位置和元素值的,dat文本文件，再用此命
令转化为稀疏矩阵。
– 例：无规律稀疏矩阵的建立。
首先编制文本文件 sp.dat如下,
5 1 5.00
3 5 8.00
4 4 2.00
5 5 0
>> load sp.dat
>> spconvert(sp)
ans =
(5,1) 5
(4,4) 2
(3,5) 8
>> full(ans)
ans =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 8
0 0 0 2 0
5 0 0 0 0
? 稀疏矩阵的计算,
同满矩阵比较，稀疏矩阵在算法上有很大的不
同。具体表现在存储空间减少，计算时间减少。
例：比较求解下面方程组 n＝ 1000时两种方法的差别。
1
2
4 1 1
1 4 1
1
1 4 1
n
x
x
x
??? ? ? ?
??? ? ? ?
??? ? ? ??
??? ? ? ?
??? ? ? ?? ? ? ???
? ? ? ???
O
MO O M
>> n=1000;
>> a1=sparse(1:n,1:n,4*ones(1,n),n,n);
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n);
>> a=a1+a2+a2';
>> b=ones(1000,1);
>> tic; x=a\b; t1=toc
t1 =
0.4800
>> a=full(a);
>> tic; x=a\b; t2=toc
t2 =
1.3220
5.6 矩阵的特征值问题
5.6.1一般矩阵的特征值与特征向量
? 格式,d=eig (A)
只求解特征值。
? 格式,[V,D]=eig (A)
求解特征值和特征向量。
? 例,
直接求解,
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> eig(A)
ans =
34.0000
8.9443
-8.9443
0.0000
精确解,
>> eig(sym(A))
ans =
[ 0]
[ 34]
[ 4*5^(1/2)]
[ -4*5^(1/2)]
高精度数值解,
>> vpa(ans,70)
ans =
[ 0]
[ 34.]
[ 8.9442719099991587856366946749251049417624734384461028
97083588981642084]
[ -8.94427190999915878563669467492510494176247343844610
28 97083588981642084]
? 同时求出特征值与特征向量,
直接求解,
>> [v,d] = eig(A)
v =
-0.5000 -0.8236 0.3764 -0.2236
-0.5000 0.4236 0.0236 -0.6708
-0.5000 0.0236 0.4236 0.6708
-0.5000 0.3764 -0.8236 0.2236
d =
34.0000 0 0 0
0 8.9443 0 0
0 0 -8.9443 0
0 0 0 0.0000
解析解,
>> [v,d]=eig(sym(A))
v =
[ -1,1,-8*5^(1/2)-17,8*5^(1/2)-17 ]
[ -3,1,4*5^(1/2)+9,-4*5^(1/2)+9 ]
[ 3,1,1,1 ]
[ 1,1,4*5^(1/2)+7,-4*5^(1/2)+7 ]
d =
[ 0,0,0,0]
[ 0,34,0,0]
[ 0,0,4*5^(1/2),0]
[ 0,0,0,-4*5^(1/2)]
5.6.2 矩阵的广义特征向量问题
若 B=I，则化成普通矩阵特征值问题。
（避免有重特征根时，特征向量矩阵为奇异
矩阵。）
? 格式,d=eig (A)
求解广义特征值。
? 格式,[V,D]=eig (A)
求解广义特征值和特征向量。
? 例,
直接求解,
>> A=[5,7,6,5; 7,10,8,7; 6,8,10,9; 5,7,9,10];
>> B=[2,6,-1,-2; 5,-1,2,3; -3,-4,1,10; 5,-2,-3,8];
>> [V,D] = eig(A,B)
V =
0.3697 -0.3741 + 0.6259i -0.3741 - 0.6259i 1.0000
0.9948 -0.0674 - 0.2531i -0.0674 + 0.2531i -0.6090
0.7979 0.9239 + 0.0264i 0.9239 - 0.0264i -0.2316
1.0000 -0.6599 - 0.3263i -0.6599 + 0.3263i 0.1319
D =
4.7564 0 0 0
0 0.0471 + 0.1750i 0 0
0 0 0.0471 - 0.1750i 0
0 0 0 -0.0037
检验,
>> norm(A*V-B*V*D)
ans =
1.3897e-014
符号运算工具箱中的 eig( )函数不支持广
义特征值的运算。
? 矩阵
? 直接解法
? 迭代法
? 符号解法
? 稀疏矩阵技术
? 特征值与特征向量
5.1 矩阵
5.1.1特殊矩阵的输入
? 数值矩阵的输入
– 零矩阵、幺矩阵及单位矩阵
生成 n?n方阵,
A=zeros(n),B=ones(n),C=eye(n)
生成 m?n矩阵,
A=zeros(m,n),B=ones(m,n),C=eye(m,n)
生成和矩阵 B同样位数的矩阵,
A=zeros(size(B))
– 随机元素矩阵
若矩阵随机元素满足 [0,1]区间上的均匀分布
生成 n?m阶标准均匀分布为随机数矩阵,
A=rand(n,m)
生成 n?n阶标准均匀分布为随机数方阵,
A=rand(n)
– 对角元素矩阵
已知向量生成对角矩阵,
A=diag(V)
已知矩阵提取对角元素列向量,
V＝ diag(A)
生成主对角线上第 k条对角线为 V的矩阵,
A=diag(V,k)
? 例,diag( )函数的不同调用格式
>> C=[1 2 3]; V=diag(C) % 生成对角矩阵
V =
1 0 0
0 2 0
0 0 3
>> V1=diag(V)' % 将列向量通过转置变换成行向量
V1 =
1 2 3
>> C=[1 2 3]; V=diag(C,2) % 主对角线上第 k条对角线为 C的矩阵
V =
0 0 1 0 0
0 0 0 2 0
0 0 0 0 3
0 0 0 0 0
0 0 0 0 0
生成三对角矩阵,
>> V=diag([1 2 3 4])+diag([2 3 4],1)+diag([5 4 3],-1)
V =
1 2 0 0
5 2 3 0
0 4 3 4
0 0 3 4
– Hilbert矩阵及逆 Hilbert矩阵
生成 n阶的 Hilbert矩阵,
A=hilb(n)
求取逆 Hilbert矩阵,
B=invhilb(n)
? 符号矩阵的输入
数值矩阵 A转换成符号矩阵,
B=sys(A)
例,
>> A=hilb(3)
A =
1.0000 0.5000 0.3333
0.5000 0.3333 0.2500
0.3333 0.2500 0.2000
>> B=sym(A)
B =
[ 1,1/2,1/3]
[ 1/2,1/3,1/4]
[ 1/3,1/4,1/5]
5.1.2 矩阵基本概念与性质
? 行列式
格式, d=det(A)
例：求行列式
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1]; det(A)
ans =
0
? 例,
>> tic,A=sym(hilb(20)); det(A),toc
ans =
1/23774547167685345090916442434276164401754
1983775348649303318533123441975931064458
5187585766816573773440565759867265558971
7656384197107933033865823241498112410235
5448916615471780963525779783680000000000
0000000000000000000000000
elapsed_time =
2.3140
高阶的 Hilbert矩阵是接近奇异的矩阵。
? 矩阵的迹
格式,t=trace(A)
? 矩阵的秩
格式,r=rank(A) ％用默认的精度求数值秩
r=rank(A,) ％给定精度下求数值秩 ?
? 例
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1]; rank(A)
ans =
3
该矩阵的秩为 3，小于矩阵的阶次，故为非满秩矩阵。
? 例
>> H=hilb(20); rank(H) ％数值方法
ans =
13
>> H=sym(hilb(20)); rank(H) % 解析方法，原矩阵为非奇异矩阵
ans =
20
? 矩阵范数
? 矩阵的范数定义,
格式,
N=norm(A) ％求解默认的 2范数
N=norm(A，选项 ) ％选项可为 1,2,inf等
? 例：求一向量、矩阵的范数
>> a=[16 2 3 13];
>> [norm(a),norm(a,2),norm(a,1),norm(a,Inf)]
ans =
2.092844953645635e+001 2.092844953645635e+001
3.400000000000000e+001 1.600000000000000e+001
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> [norm(A),norm(A,2),norm(A,1),norm(A,Inf)]
ans =
34 34 34 34
符号运算工具箱未提供 norm( )函数，需先用 double( )
函数转换成双精度数值矩阵，再调用 norm( )函数。
? 特征多项式
格式,C=poly(A)
例,>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> poly(A) ％直接求取
ans =
1.000000000000000e+000 -3.399999999999999e+001
-7.999999999999986e+001 2.719999999999999e+003
-2.819840539024018e-012
>> A=sym(A); poly(A) ％运用符号工具箱
ans =
x^4-34*x^3-80*x^2+2720*x
? 符号多项式与数值多项式的转换
格式,
f=poly2sym(P) 或 f=poly2sym(P,x)
格式,P=sym2poly(f)
? 例,
>> P=[1 2 3 4 5 6]; % 先由系数按降幂顺序排列表示多
项式
>> f=poly2sym(P,'v') % 以 v 为算子表示多项式
f =
v^5+2*v^4+3*v^3+4*v^2+5*v+6
>> P=sym2poly(f)
P =
1 2 3 4 5 6
? 矩阵的逆矩阵
格式,C=inv(A)
例,
>> format long; H=hilb(4); H1=inv(H)
H1 =
1.0e+003 *
0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999
-0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984
0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961
-0.13999999999999 1.67999999999984 -4.19999999999961 2.79999999999974
检验,
>> H*H1
ans =
1.00000000000001 0.00000000000023 -0.00000000000045 0.00000000000023
0.00000000000001 1.00000000000011 -0.00000000000011 0.00000000000011
0.00000000000001 0 1.00000000000011 0
0.00000000000000 0.00000000000011 -0.00000000000011 1.00000000000011
计算误差范数,
>> norm(H*inv(H)-eye(size(H)))
ans =
6.235798190375727e-013
>> H2=invhilb(4); norm(H*H2-eye(size(H)))
ans =
5.684341886080802e-014
>> H=hilb(10); H1=inv(H); norm(H*H1-eye(size(H)))
ans =
0.00264500826202
>> H2=invhilb(10); norm(H*H2-eye(size(H)))
ans =
1.612897415528547e-005
>> H=hilb(13); H1=inv(H); norm(H*H1-eye(size(H)))
Warning,Matrix is close to singular or badly scaled,
Results may be inaccurate,RCOND = 2.339949e-018,
ans =
53.23696008570294
>> H2=invhilb(13); norm(H*H2-eye(size(H)))
ans =
11.37062973181391
对接近于奇异矩阵，高阶一般不建议用 inv( )，可用符号工具箱。
>> H=sym(hilb(7)); inv(H)
ans =
[ 49,-1176,8820,-29400,48510,-38808,12012]
[-1176,37632,-317520,1128960,-1940400,1596672,-504504]
[8820,-317520,2857680,-10584000,18711000,-15717240,5045040]
[-29400,1128960,-10584000,40320000,-72765000,62092800,-20180160]
[48510,-1940400,18711000,-72765000,133402500,-115259760,37837800]
[-38808,1596672,-15717240,62092800,-115259760,100590336,-33297264]
[12012,-504504,5045040,-20180160,37837800,-33297264,11099088]
>> H=sym(hilb(30)); norm(double(H*inv(H)-eye(size(H))))
ans =
0
? 例：奇异阵求逆
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> format long; B = inv(A)
Warning,Matrix is close to singular or badly scaled,
Results may be inaccurate,RCOND = 1.306145e-017,
B =
1.0e+014 *
0.93824992236885 2.81474976710656 -2.81474976710656 -0.93824992236885
2.81474976710656 8.44424930131968 -8.44424930131968 -2.81474976710656
-2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656
-0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885
>> norm(A*B-eye(size(A))) ％检验
ans =
1.64081513306419
>> A=sym(A); inv(A) ％奇异矩阵不存在一个相应的逆矩阵，用符号工具
箱的函数也不行
Error using ==> sym/inv
Error,(in inverse) singular matrix
5.2 线性方程组直接解法
5.2.1线性方程组直接求解－矩阵除法
? 关于线性方程组的直接解法，如 Gauss消去法、
选主元消去法、平方根法、追赶法等等，在
MATLAB中，只需用“／”或,\”就解决问
题。它内部实际包含着许许多多的自适应算
法，如对超定方程用最小二乘法，对欠定方
程时它将给出范数最小的一个解，解三对角
阵方程组时用追赶法等等。
格式,
x=A\b
? 例：解方程组
>> A=[.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129;
.3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927];
>> b=[0.4043 0.1550 0.4240 -0.2557]';
>> x=A\b; x'
ans =
-0.1819 -1.6630 2.2172 -0.4467
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4
0,4 0 9 6 0,1 2 3 4 0,3 6 7 8 0,2 9 4 3 0,4 0 4 3
0,2 2 4 6 0,3 8 7 2 0,4 0 1 5 0,1 1 2 9 0,1 1 5 0
0,3 6 4 5 0,1 9 2 0 0,3 7 8 1 0,0 6 4 3 0,4 2 4 0
0,1 7 8 4 0,4 0 0 2 0,2 7 8 6 0,3 9 2 7 0,2 5 5 7
x x x x
x x x x
x x x x
x x x x
? ? ? ??
?
? ? ? ??
?
? ? ? ?
?
? ? ? ? ? ?
?
5.2.2线性方程组直接求解－判定求解
? 例,
>> A=[1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2]; B=[5 1; 4 2; 3 3; 2 4];
>> C=[A B]; rank(A),rank(C)
ans =
4
ans =
4
>> x=inv(A)*B
x =
-1.8000 2.4000
1.8667 -1.2667
3.8667 -3.2667
-2.1333 2.7333
检验
>> norm(A*x-B)
ans =
7.4738e-015
精确解
>> x1=inv(sym(A))*B
x1 =
[ -9/5,12/5]
[ 28/15,-19/15]
[ 58/15,-49/15]
[ -32/15,41/15]
检验
>> norm(double(A*x1-B))
ans =
0
原方程组对应的齐次方程组的解
求取 A矩阵的化零矩阵,
格式,Z=null(A)
求取 A矩阵的化零矩阵的规范形式,
格式,Z=null(A,‘ r ’)
? 例,
判断可解性
>> A=[1 2 3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2]; B=[1;3;2;6];
>> C=[A B]; [rank(A),rank(C)]
ans =
2 2
>> Z=null(A,'r') % 解出规范化的化零空间
Z =
2.0000 3.0000
-2.5000 -3.5000 ％
1.0000 0
0 1.0000
1
23
34
4
23
2.5 3.5
10
01
x
xx
xx
x
?? ??
?? ??
?? ??
?? ???
??
?? ?? ??
?? ??
?? ????
>> x0=pinv(A)*B % 得出一个特解
x0 =
0.9542
0.7328 %全部解
-0.0763
-0.2977
验证得出的解
>> a1=randn(1); a2=rand(1); % 取不同分布的随
机数
>> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(A*x-B)
ans =
4.4409e-015
1
2
12
3
4
2 3 0.9542
2.5 3.5 0.7328
1 0 - 0.0763
0 1 - 0.2977
x
x
aa
x
x
?? ? ? ? ? ? ?
?? ? ? ? ? ? ?
??
?? ? ? ? ? ? ?? ? ?
?? ? ? ? ? ? ?
?? ? ? ? ? ? ?
?? ? ? ? ? ? ???
? 解析解
>> Z=null(sym(A))
Z =
[ 2,3]
[ -5/2,-7/2]
[ 1,0]
[ 0,1]
>> x0=sym(pinv(A)*B)
x0 =
[ 125/131]
[ 96/131]
[ -10/131] ％
[ -39/131]
? 验证得出的解
>> a1=randn(1); a2=rand(1); % 取不同分布的随机数
>> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(double(A*x-B))
ans =
0
? 通解
>> syms a1 a2;
>> x=a1*Z(:,1)+a2*Z(:,2)+x0
x =
[ 2*a1+3*a2+125/131]
[ -5/2*a1-7/2*a2+96/131]
[ a1-10/131]
[ a2-39/131]
摩尔－彭罗斯广义逆求解出的方程最小二乘解
不满足原始代数方程。
5.2.3 线性方程组的直接求解分析
? LU分解
? 格式
[l,u,p]=lu(A)
L是一个单位下三角矩阵,u是一个上三角矩
阵,
p是代表选主元的置换矩阵。
故,Ax=y => PAx=Py
=> LUx=Py => PA=LU
[l,u]=lu(A)
其中 l等于 P-1 L,u等于 U，所以 (P-1 L)U=A
? 例：对 A进行 LU分解
>> A=[1 2 3; 2 4 1; 4 6 7];
>> [l,u,p]=lu(A)
l =
1.0000 0 0
0.5000 1.0000 0
0.2500 0.5000 1.0000
u =
4.0000 6.0000 7.0000
0 1.0000 -2.5000
0 0 2.5000
p =
0 0 1
0 1 0
1 0 0
1 2 3
2 4 1
467
A
??
???
??
??
>> [l,u]=lu(A) ％ l＝ P-1 L
l =
0.2500 0.5000 1.0000
0.5000 1.0000 0
1.0000 0 0
u =
4.0000 6.0000 7.0000
0 1.0000 -2.5000
0 0 2.5000
? Cholesky分解
若矩阵 A为 n阶对称正定阵,则存在唯
一的对角元素为正的三角阵 D，使得
格式,
D=chol(A)
例：进行 Cholesky分解。
>> A=[16 4 8; 4 5 -4; 8 -4 22];
>> D=chol(A)
D =
4 1 2
0 2 -3
0 0 3
1 6 4 8
4 5 4
8 4 2 2
A
??
??
??
??
???
? 三个变换
在线性方程组的迭代求解中，要用到系数
矩阵 A的上三角矩阵、对角阵和下三角矩阵。
此三个变换在 MATLAB中可由以下函数实现。
– 上三角变换,
格式 triu(A,1)
– 对角变换,
格式 diag(A)
– 下三角变换,
格式 tril(A,-1)
例：对此矩阵做三种变换。
1 2 2
1 1 1
2 2 1
A
???
??
? ??
??
>> A=[1 2 -2;1 1 1;2 2 1]; ％
>> triu(A,1)
ans =
0 2 -2
0 0 1
0 0 0
>> tril(A,-1)
ans =
0 0 0
1 0 0
2 2 0
>> b=diag(A); b'
ans =
1 1 1
1 2 2
1 1 1
2 2 1
A
???
??
? ??
??
5.3 迭代解法的几种形式
5.3.1 Jacobi迭代法
? 方程组 Ax=b
A可写成 A=D-L-U
其中,D=diag[a11,a22,…,a nn],-L,-U分别为 A的
严格下、上三角部分（不包括对角线元素）,
由 Ax=b ＝ > x=Bx+f
由此可构造迭代法,
x(k+1)=Bx(k)+f
其中,B=D-1(L+U)=I-D-1A,f=D-1b,
function y=jacobi(a,b,x0)
D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1);
B=D\(L+U); f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=B*x0+f;
n=n+1;
end
n
? 例：用 Jacobi方法求解,
设 x(0)=0,精度为 10-6。
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> jacobi(a,b,[0;0;0])
n =
11
ans =
0.9958
0.9579
0.7916
12
1 2 3
23
10 9
10 2 7
2 10 6
xx
x x x
xx
???
?
? ? ? ??
? ? ? ?
?
5.3.2 Gauss-Seidel迭代法
由原方程构造迭代方程
x(k+1)=G x(k)+f
其中,G=(D-L)-1 U,f=(D-L)-1 b
D=diag[a11,a22,…,a nn],
-L,-U分别为 A的严格下、上三角部
分（不包括对角线元素）,
function y=seidel(a,b,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
G=(D-L)\U ;f=(D-L)\b;
y=G*x0+f; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=G*x0+f;
n=n+1;
end
n
? 例：对上例用 Gauss-Seidel迭代法求解
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> seidel(a,b,[0;0;0])
n =
7
ans =
0.9958
0.9579
0.7916
例：分别用 Jacobi和 G-S
法迭代求解,看是否收敛。
1
2
3
1 2 2 9
1 1 1 7
2 2 1 6
x
x
x
?? ? ? ? ? ?
? ? ? ? ? ?
?? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
>> a=[1 2 -2; 1 1 1; 2 2 1]; b=[9; 7; 6];
>> jacobi(a,b,[0;0;0])
n =
4
ans =
-27
26
8
>> seidel(a,b,[0;0;0])
n =
1011
ans =
1.0e+305 *
-Inf
Inf
-1.7556
5.3.3 SOR迭代法
在很多情况下,J法和 G-S法收敛较
慢，所以考虑对 G-S法进行改进。于是引
入一种新的迭代法－逐次超松弛迭代法
（ Succesise Over-Relaxation),记为 SQR法。
迭代公式为,
X(k+1)= (D-wL)-1((1-w)D+wU)x(k)
+ w(D-wL)-1 b
其中,w最佳值在 [1,2）之间，不易计
算得到，因此 w通常有经验给出。
function y=sor(a,b,w,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
M=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w;
y=M*x0+f; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=M*x0+f;
n=n+1;
end
n
例：上例中，当 w=1.103时，用 SOR法求解
原方程。
>> a=[10 -1 0; -1 10 -2; 0 -2 10];
>> b=[9; 7; 6];
>> sor(a,b,1.103,[0;0;0])
n =
8
ans =
0.9958
0.9579
0.7916
5.3.4 两步迭代法
当线性方程系数矩阵为对称正定时，
可用一种特殊的迭代法来解决，其迭代
公式为,
（ D-L） x(k+1/2) =U x(k) +b
（ D-U） x(k+1)=Lx(k+1/2) +b
=>
x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 b
x(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 b
function y=twostp(a,b,x0)
D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);
G1==(D-L)\U; f1=(D-L)\b;
G2==(D-U)\L; f1=(D-U)\b;
y=G1*x0+f1; y=G2*y+f2; n=1;
while norm(y-x0)>=1.0e-6
x0=y;
y=G1*x0+f1; y=G2*y+f2;
n=n+1;
end
n
? 例：求解方程组
>> a=[10 -1 2 0; -1 11 -1 3; 2 -1 10 3; 0 3 -1 8];
>> b=[6; 25; -11; 15];
>> twostp(a,b,[0; 0; 0; 0])
n =
7
ans =
1.0791
1.9824
-1.4044
0.9560
1
2
3
4
10 1 2 0 6
1 11 1 3 25
2 1 10 1 11
0 3 1 8 15
x
x
x
x
? ??? ? ? ?
??? ? ? ?
?? ??? ? ? ?
?
??? ? ? ?? ? ?
??? ? ? ?? ? ? ?
?? ? ? ???
5.4 线性方程组的符号解法
在 MATLAB的 Symbolic Toolbox中提供了
线性方程的符号求解函数，如
linsolve(A,b)
等同于 X = sym(A)\sym(b),
solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN ')
? 例,
>> A=sym('[10,-1,0;-1,10,-2;0,-2,10]')；
>> b=('[9; 7; 6]')；
>> linsolve(A,b)
ans =
[ 473/475]
[ 91/95]
[ 376/475]
>> vpa(ans)
ans =
[,99578947368421052631578947368421]
[,95789473684210526315789473684211]
[,79157894736842105263157894736842]
例,
>> [x,y] = solve('x^2 + x*y + y = 3','x^2
- 4*x + 3 = 0','x','y')
x =
[ 1]
[ 3]
y =
[ 1]
[ -3/2]
5.5 稀疏矩阵技术
? 稀疏矩阵的建立,
– 格式 S=sparse(i,j,s,m,n)
生成一 mxn阶的稀疏矩阵，以向量 i和 j为坐标的位置上对应
元素值为 s。
– 例,
>> n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n)
a1 =
(1,1) 4
(2,2) 4
(3,3) 4
(4,4) 4
(5,5) 4
– 例,
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n)
a2 =
(2,1) 1
(3,2) 1
(4,3) 1
(5,4) 1
>> full(a2)
– ans =
0 0 0 0 0
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
? 例,n=5,建立主对角线上元素为 4，两条次对角
线为 1的三对角阵。
>> n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n);
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n);
>> a=a1+a2+a2'
a =
(1,1) 4
(2,1) 1
(1,2) 1
(2,2) 4
(3,2) 1
(2,3) 1
(3,3) 4
(4,3) 1
(3,4) 1
(4,4) 4
(5,4) 1
(4,5) 1
(5,5) 4
>> full(a)
ans =
4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4
– 格式 A=spdiags(B,d,m,n)
生成一 mxn阶的稀疏矩阵，使得 B的列放在由 d指定的
位置。
– 例,
>> n=5
>> b=spdiags([ones(n,1),4*ones(n,1),ones(n,1)],…
[-1,0,1],n,n);
>> full(b)
ans =
4 1 0 0 0
1 4 1 0 0
0 1 4 1 0
0 0 1 4 1
0 0 0 1 4
– 格式,spconvert(dd)
对于无规律的稀疏矩阵，可使用此命令
由外部数据转化为稀疏矩阵。
– 调用形式为：先用 load函数加载以行表示对
应位置和元素值的,dat文本文件，再用此命
令转化为稀疏矩阵。
– 例：无规律稀疏矩阵的建立。
首先编制文本文件 sp.dat如下,
5 1 5.00
3 5 8.00
4 4 2.00
5 5 0
>> load sp.dat
>> spconvert(sp)
ans =
(5,1) 5
(4,4) 2
(3,5) 8
>> full(ans)
ans =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 8
0 0 0 2 0
5 0 0 0 0
? 稀疏矩阵的计算,
同满矩阵比较，稀疏矩阵在算法上有很大的不
同。具体表现在存储空间减少，计算时间减少。
例：比较求解下面方程组 n＝ 1000时两种方法的差别。
1
2
4 1 1
1 4 1
1
1 4 1
n
x
x
x
??? ? ? ?
??? ? ? ?
??? ? ? ??
??? ? ? ?
??? ? ? ?? ? ? ???
? ? ? ???
O
MO O M
>> n=1000;
>> a1=sparse(1:n,1:n,4*ones(1,n),n,n);
>> a2=sparse(2:n,1:n-1,ones(1,n-1),n,n);
>> a=a1+a2+a2';
>> b=ones(1000,1);
>> tic; x=a\b; t1=toc
t1 =
0.4800
>> a=full(a);
>> tic; x=a\b; t2=toc
t2 =
1.3220
5.6 矩阵的特征值问题
5.6.1一般矩阵的特征值与特征向量
? 格式,d=eig (A)
只求解特征值。
? 格式,[V,D]=eig (A)
求解特征值和特征向量。
? 例,
直接求解,
>> A=[16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1];
>> eig(A)
ans =
34.0000
8.9443
-8.9443
0.0000
精确解,
>> eig(sym(A))
ans =
[ 0]
[ 34]
[ 4*5^(1/2)]
[ -4*5^(1/2)]
高精度数值解,
>> vpa(ans,70)
ans =
[ 0]
[ 34.]
[ 8.9442719099991587856366946749251049417624734384461028
97083588981642084]
[ -8.94427190999915878563669467492510494176247343844610
28 97083588981642084]
? 同时求出特征值与特征向量,
直接求解,
>> [v,d] = eig(A)
v =
-0.5000 -0.8236 0.3764 -0.2236
-0.5000 0.4236 0.0236 -0.6708
-0.5000 0.0236 0.4236 0.6708
-0.5000 0.3764 -0.8236 0.2236
d =
34.0000 0 0 0
0 8.9443 0 0
0 0 -8.9443 0
0 0 0 0.0000
解析解,
>> [v,d]=eig(sym(A))
v =
[ -1,1,-8*5^(1/2)-17,8*5^(1/2)-17 ]
[ -3,1,4*5^(1/2)+9,-4*5^(1/2)+9 ]
[ 3,1,1,1 ]
[ 1,1,4*5^(1/2)+7,-4*5^(1/2)+7 ]
d =
[ 0,0,0,0]
[ 0,34,0,0]
[ 0,0,4*5^(1/2),0]
[ 0,0,0,-4*5^(1/2)]
5.6.2 矩阵的广义特征向量问题
若 B=I，则化成普通矩阵特征值问题。
（避免有重特征根时，特征向量矩阵为奇异
矩阵。）
? 格式,d=eig (A)
求解广义特征值。
? 格式,[V,D]=eig (A)
求解广义特征值和特征向量。
? 例,
直接求解,
>> A=[5,7,6,5; 7,10,8,7; 6,8,10,9; 5,7,9,10];
>> B=[2,6,-1,-2; 5,-1,2,3; -3,-4,1,10; 5,-2,-3,8];
>> [V,D] = eig(A,B)
V =
0.3697 -0.3741 + 0.6259i -0.3741 - 0.6259i 1.0000
0.9948 -0.0674 - 0.2531i -0.0674 + 0.2531i -0.6090
0.7979 0.9239 + 0.0264i 0.9239 - 0.0264i -0.2316
1.0000 -0.6599 - 0.3263i -0.6599 + 0.3263i 0.1319
D =
4.7564 0 0 0
0 0.0471 + 0.1750i 0 0
0 0 0.0471 - 0.1750i 0
0 0 0 -0.0037
检验,
>> norm(A*V-B*V*D)
ans =
1.3897e-014
符号运算工具箱中的 eig( )函数不支持广
义特征值的运算。
