第四章 线性代数问题求解 矩阵 线性方程组的直接解法 线性方程组的迭代法 线性方程组的符号解法 稀疏矩阵技术 特征值与特征向量.

Presentation on theme: "第四章 线性代数问题求解 矩阵 线性方程组的直接解法 线性方程组的迭代法 线性方程组的符号解法 稀疏矩阵技术 特征值与特征向量."— Presentation transcript:

1 第四章 线性代数问题求解 矩阵 线性方程组的直接解法 线性方程组的迭代法 线性方程组的符号解法 稀疏矩阵技术 特征值与特征向量

2 4.1 矩阵 4.1.1特殊矩阵的输入 数值矩阵的输入 零矩阵、幺矩阵及单位矩阵 生成nn方阵：
4.1 矩阵 4.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))

3 随机元素矩阵 若矩阵随机元素满足[0,1]区间上的均匀分布 生成nm阶标准均匀分布为随机数矩阵： A=rand(n,m) 生成nn阶标准均匀分布为随机数方阵： A=rand(n)

4 对角元素矩阵 已知向量生成对角矩阵： A=diag(V) 已知矩阵提取对角元素列向量： V＝diag(A) 生成主对角线上第k条对角线为V的矩阵： A=diag(V,k)

5 例：diag( )函数的不同调用格式 >> C=[1 2 3]; V=diag(C) % 生成对角矩阵 V = 1 0 0
>> V1=diag(V)' % 将列向量通过转置变换成行向量 V1 = >> C=[1 2 3]; V=diag(C,2) % 主对角线上第 k条对角线为C的矩阵

6 生成三对角矩阵: >> V=diag([ ])+diag([2 3 4],1)+diag([5 4 3],-1) V =

7 Hilbert矩阵及逆Hilbert矩阵
生成n阶的Hilbert矩阵： A=hilb(n) 求取逆Hilbert矩阵： B=invhilb(n)

8 Hankel(汉克 ) 矩阵 其中：第一列的各个元素定义为C向量，最后一行各个元素定义为R。H为对称阵。 H1=hankel(C)
由 Hankel 矩阵反对角线上元素相等得出一下三角阵均为零的Hankel 矩阵

9 Vandermonde(范德蒙)矩阵

10 伴随矩阵 其中：P(s)为首项系数为一的多向式。

11 符号矩阵的输入 数值矩阵A转换成符号矩阵： B=sym(A) 例： >> A=hilb(3) A = >> B=sym(A) B = [ 1, 1/2, 1/3] [ 1/2, 1/3, 1/4] [ 1/3, 1/4, 1/5]

12 4.1.2 矩阵基本概念与性质 行列式 格式 ：d=det(A) 例：求行列式
>> A=[ ; ; ; ]; det(A) ans =

13 例： >> tic, A=sym(hilb(20)); det(A), toc ans = 1/ elapsed_time = 2.3140 高阶的Hilbert矩阵是接近奇异的矩阵。

14 格式：r=rank(A) ％用默认的精度求数值秩 r=rank(A， ) ％给定精度下求数值秩
矩阵的迹 格式： t=trace(A) 矩阵的秩 格式：r=rank(A) ％用默认的精度求数值秩 r=rank(A， ) ％给定精度下求数值秩 矩阵的秩也表示该矩阵中行列式不等于0的子式的最大阶次。可证行秩和列秩（线性无关的）应相等。

15 例 >> 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) ％数值方法 13 >> H=sym(hilb(20)); rank(H) % 解析方法，原矩阵为非奇异矩阵 20

16 矩阵范数

17 矩阵的范数定义： 格式： N=norm(A) ％求解默认的2范数 N=norm(A，选项) ％选项可为1,2,inf等

18 例：求一向量、矩阵的范数 >> a=[16 2 3 13];
>> [norm(a), norm(a,2), norm(a,1), norm(a,Inf)] ans = e e e e+001 >> A=[ ; ; ; ]; >> [norm(A), norm(A,2), norm(A,1), norm(A,Inf)] 符号运算工具箱未提供norm( )函数，需先用double( )函数转换成双精度数值矩阵，再调用norm( )函数。

19 特征多项式 格式： C=poly(A) 例：>> A=[ ; ; ; ]; >> poly(A) ％直接求取 ans = e e+001 e e+003 e-012 >> A=sym(A); poly(A) ％运用符号工具箱 x^4-34*x^3-80*x^2+2720*x

20 矩阵多项式的求解

21 符号多项式与数值多项式的转换 格式： f=poly2sym(P) 或 f=poly2sym(P,x) 格式： P=sym2poly(f)

22 例： >> 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 =

23 矩阵的逆矩阵 格式： C=inv(A) 例： >> format long; H=hilb(4); H1=inv(H) H1 =
1.0e+003 *

24 >> norm(H*inv(H)-eye(size(H)))
检验： >> H*H1 ans = 计算误差范数： >> norm(H*inv(H)-eye(size(H))) e-013 >> H2=invhilb(4); norm(H*H2-eye(size(H))) e-014

25 >> H=hilb(10); H1=inv(H); norm(H*H1-eye(size(H)))
ans = >> H2=invhilb(10); norm(H*H2-eye(size(H))) e-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 = e-018. >> H2=invhilb(13); norm(H*H2-eye(size(H))) 对接近于奇异矩阵，高阶一般不建议用inv( )，可用符号工具箱。

26 >> H=sym(hilb(7)); inv(H)
ans = [ 49, , , , , , ] [-1176, , , , , , ] [8820, , , , , , ] [-29400, , , , , , ] [48510, , , , , , ] [-38808, , , , , , ] [12012, , , , , , ] >> H=sym(hilb(30)); norm(double(H*inv(H)-eye(size(H))))

27 例：奇异阵求逆 >> 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 = e-017. B = 1.0e+014 * >> norm(A*B-eye(size(A))) ％检验 ans = >> A=sym(A); inv(A) ％奇异矩阵不存在一个相应的逆矩阵，用符号工具箱的函数也不行 ??? Error using ==> sym/inv Error, (in inverse) singular matrix

28 同样适用于含有变量的矩阵求逆。 例： >> syms a1 a2 a3 a4; >> C=[a1 a2;a3 a4]； >> inv(C) ans = [ -a4/(-a1*a4+a2*a3), a2/(-a1*a4+a2*a3)] [ a3/(-a1*a4+a2*a3), -a1/(-a1*a4+a2*a3)]

29 矩阵的相似变换与正交矩阵 其中：A为一方阵，B矩阵非奇异。 相似变换后，X矩阵的秩、迹、行列式与特征值等均不发生变化，其值与A矩阵完全一致。
对于一类特殊的相似变换满足如下条件，称为正交基矩阵。

30 例： >> A=[5,9,8,3; 0,3,2,4; 2,3,5,9; 3,4,5,8]; >> Q=orth(A)
>> norm(Q'*Q-eye(4)) ans = 4.6395e-016 >> norm(Q*Q'-eye(4)) 4.9270e-016

31 例： >> A=[16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1]; >> Q=orth(A) ％A为奇异矩阵，故得出的Q为长方形矩阵 Q = >> norm(Q'*Q-eye(3)) ans = 1.0140e-015

32 4.2 线性方程组直接解法 4.2.1线性方程组直接求解－矩阵除法
关于线性方程组的直接解法，如Gauss消去法、选主元消去法、平方根法、追赶法等等，在MATLAB中，只需用“／”或“\”就解决问题。它内部实际包含着许许多多的自适应算法，如对超定方程用最小二乘法，对欠定方程时它将给出范数最小的一个解，解三对角阵方程组时用追赶法等等。 格式： x=A\b

33 例：解方程组 >> A=[.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129; .3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927]; >> b=[ ]'; >> x=A\b; x' ans =

34 4.2.2线性方程组直接求解－判定求解

35

36 例： >> A=[ ; ; ; ]; B=[5 1; 4 2; 3 3; 2 4]; >> C=[A B]; rank(A), rank(C) ans = 4 >> x=inv(A)*B x =

37 检验 精确解 >> 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))

38 原方程组对应的齐次方程组的解 求取A矩阵的化零矩阵： 格式： Z=null(A) 求取A矩阵的化零矩阵的规范形式： 格式： Z=null(A, ‘ r ’)

39 例： 判断可解性 >> 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 = >> Z=null(A,'r') % 解出规范化的化零空间 Z = ％

40 >> x0=pinv(A)*B % 得出一个特解
0.9542 %全部解 验证得出的解 >> a1=randn(1); a2=rand(1); % 取不同分布的随机数 >> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(A*x-B) ans = 4.4409e-015

41 解析解 >> 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]

42 验证得出的解 >> a1=randn(1); a2=rand(1); % 取不同分布的随机数 >> x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(double(A*x-B)) ans = 通解 >> 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]

43 摩尔－彭罗斯广义逆求解出的方程最小二乘解不满足原始代数方程。

44 4.2.3 线性方程组的直接求解分析 LU分解

45

46

47 L是一个单位下三角矩阵，u是一个上三角矩阵， p是代表选主元的置换矩阵。 故：Ax=y => PAx=Py
格式 [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

48 例：对A进行LU分解 >> A=[1 2 3; 2 4 1; 4 6 7]; >> [l,u,p]=lu(A)
u = p =

49 >> [l,u]=lu(A) ％ l＝P-1 L
u =

50 QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。 求得正交矩阵Q和上三角阵R，Q和R满足A=QR。 格式： [Q,R] = qr(A)

51 例： >> A =[ ; ; ; ]; >> [Q,R] = qr(A) Q = R =

52 Cholesky(乔里斯基 )分解 若矩阵A为 n阶对称正定阵,则存在唯一的对角元素为正的三角阵D，使得

53 格式： D=chol(A)

54 例：进行Cholesky分解。 >> A=[16 4 8; ; ]; >> D=chol(A) D =

55 利用矩阵的LU、QR和cholesky分解求方程组的解
A*X=b 变成 L*U*X=b 所以 X=U\(L\b) 这样可以大大提高运算速度。 例：求方程组 的一个特解。 解： >> A=[4 2 -1;3 -1 2;11 3 0]; >> B=[2 10 8]'; >> D=det(A) D =

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

57 >> X=U\(L\B) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = e-017. X = 1.0e+016 * % 结果中的警告是由于系数行列式为零产生的。 % 可以通过A*X验证其正确性。 1.4862 1.3511 >> A*X ans = 8

58 （2）Cholesky分解 若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积， 方程 A*X=b 变成 R’*R*X=b 所以 X=R\(R’\b) （3）QR分解 对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR 方程 A*X=b 变形成 QRX=b 所以 X=R\(Q\b) 这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。

59 在线性方程组的迭代求解中，要用到系数矩阵A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。
上三角变换： 格式 triu(A,1) 对角变换： 格式 diag(A) 下三角变换： 格式 tril(A,-1) 例：对此矩阵做三种变换。

60 >> A=[1 2 -2;1 1 1;2 2 1]; ％ >> triu(A,1) ans = >> tril(A,-1) >> b=diag(A); b'

61 4.3 迭代解法的几种形式 5.3.1 Jacobi迭代法 方程组 Ax=b A可写成 A=D-L-U
其中:D=diag[a11,a22,…,ann], -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.

62 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; n=n+1; end n

63 例：用Jacobi方法求解， 设x(0)=0,精度为10-6。 >> a=[ ; ; ]; >> b=[9; 7; 6]; >> jacobi(a,b,[0;0;0]) n = 11 ans = 0.9958 0.9579 0.7916

64 4.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,…,ann], -L、-U分别为A的严格下、上三角部分（不包括对角线元素）.

65 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

66 例：对上例用Gauss-Seidel迭代法求解
>> b=[9; 7; 6]; >> seidel(a,b,[0;0;0]) n = 7 ans = 0.9958 0.9579 0.7916 例：分别用Jacobi和G-S 法迭代求解,看是否收敛。

67 >> 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]) 1011 1.0e+305 * -Inf Inf

68 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通常有经验给出。

69 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

70 例：上例中，当w=1.103时，用SOR法求解原方程。
>> a=[ ; ; ]; >> b=[9; 7; 6]; >> sor(a,b,1.103,[0;0;0]) n = 8 ans = 0.9958 0.9579 0.7916

71 4.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

72 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

73 例：求解方程组 >> a=[ ; ; ; ]; >> b=[6; 25; -11; 15]; >> twostp(a, b, [0; 0; 0; 0]) n = 7 ans = 1.0791 1.9824 0.9560

74 4.4 线性方程组的符号解法 在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，如 linsolve(A,b)
等同于 X = sym(A)\sym(b). solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN ')

75 例： >> 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) [ ] [ ] [ ]

76 例： >> [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0','x','y') x = [ 1] [ 3] y = [ 1] [ -3/2]

77 4.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) (2,2) (3,3) (4,4) (5,5)

78 例： >> a2=sparse(2:n, 1:n-1,ones(1,n-1),n,n) a2 = (2,1) (3,2) (4,3) (5,4) >> full(a2) ans =

79 例：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) (2,1) (1,2) (2,2) (3,2) (2,3) (3,3) (4,3)

80 (3,4) (4,4) (5,4) (4,5) (5,5) >> full(a) ans =

81 格式 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 =

82 格式： spconvert(dd) 对于无规律的稀疏矩阵，可使用此命令由外部数据转化为稀疏矩阵。 调用形式为：先用load函数加载以行表示对应位置和元素值的.dat文本文件，再用此命令转化为稀疏矩阵。 例：无规律稀疏矩阵的建立。 首先编制文本文件sp.dat如下：

83 >> load sp.dat >> spconvert(sp) ans = (5,1) (4,4) (3,5) >> full(ans)

84 稀疏矩阵的计算： 同满矩阵比较，稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少，计算时间减少。
例：比较求解下面方程组n＝1000时两种方法的差别。

85 >> 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

86 4.6 矩阵的特征值问题 4.6.1一般矩阵的特征值与特征向量
4.6 矩阵的特征值问题 一般矩阵的特征值与特征向量 格式： d=eig (A) 只求解特征值。 格式： [V, D]=eig (A) 求解特征值和特征向量。

87 例： 直接求解： >> A=[ ; ; ; ]; >> eig(A) ans = 8.9443 0.0000

88 精确解： 高精度数值解： >> eig(sym(A)) ans = [ 0] [ 34] [ 4*5^(1/2)]
[ ] [ ] [ 4*5^(1/2)] [ -4*5^(1/2)] 高精度数值解： >> vpa(ans,70) [ ] [ ] [ ] [ ]

89 同时求出特征值与特征向量： 直接求解： >> [v, d] = eig(A) v = d =

90 解析解： >> [v,d]=eig(sym(A)) v = [ -1, , *5^(1/2)-17, *5^(1/2)-17 ] [ -3, , *5^(1/2)+9, *5^(1/2)+9 ] [ 3, , , ] [ 1, , *5^(1/2)+7, *5^(1/2)+7 ] d = [ , , , ] [ , , , ] [ , , 4*5^(1/2), ] [ , , , -4*5^(1/2)]

91 4.6.2 矩阵的广义特征向量问题 若B=I，则化成普通矩阵特征值问题。 格式： d=eig (A，B) 求解广义特征值。
格式： [V, D]=eig (A，B) 求解广义特征值和特征向量。

92 例： 直接求解： >> 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 = i i i i i i i i

93 >> norm(A*V-B*V*D) ans = 1.3897e-014
i i 检验： >> norm(A*V-B*V*D) ans = 1.3897e-014 符号运算工具箱中的eig( )函数不支持广义特征值的运算。

94