实验4(行列式、矩阵与线性变换).ppt

1、数学实验,实验四 行列式、矩阵与线性变换,实验目的,学会用matlab软件对矩阵进行一些数值计算; 学会用matlab软件解线性方程组； 掌握逆矩阵的一种应用：整数逆矩阵加密、解密方法； 熟悉三维空间中的线性变换，加深对正交变换保持距离不变性的理解。,det(A） 行列式计算 inv(A) 矩阵的逆 rank(A) 矩阵的秩 trace(A) 矩阵的迹 范数 n=norm(V) norm(V,P) = sum(abs(V).P)(1/P) norm(V) = norm(V,2) norm(V,inf) = max(abs(V) norm(V,-inf) = min(abs(V),Matlab软

2、件对矩阵操作命令,Matlab软件对矩阵操作命令,v=triu(x) 输入矩阵x，输出v是x的上三角阵 v=tril(x) 输入矩阵x，输出v是x的下三角阵 v=triu(x,k) 输入矩阵x，输出v是x的上三角阵, k0,对角线朝上平移k行； k0,对角线朝上平移k行； k0,对角线朝下平移k行；,d=eig(A) 计算矩阵A的特征值 v,d=eig(A) 特征向量和特征值， 满足AV=VD. rref(A) 通过初等行变换将A化为 最简阶梯形矩阵 rrefmovie 可观察初等行变换的过程 reshape(A,m,n) 将A改写成m行n 列的矩阵,Matlab软件对矩阵操作命令,基本操作,

3、示例1、产生一个4阶随机矩阵，执行下面的操作： （1）求其行列式，检验其是否可逆；若可逆，求其逆矩阵； （2）计算该矩阵的特征值、特征向量； （3）将该矩阵化为行最简单的阶梯形； （4）求出其各阶顺序主子式； （5）验证矩阵的特征值之和等于矩阵主对角元之和，特征值之积等于矩阵的行列式。,示例1：（程序）,clear;clc; a=rand(4) %产生一个4阶随机矩阵； hls=det(a) %(1)计算行列式； njz=inv(a) %(1)计算矩阵逆； tzxl,tzz=eig(a) %(2)计算特征向量和特征值； jtx=rref(a) %(3)化为行最简单的阶梯形；,示例1：（程序续）

4、,zzs=; for i=1:4 zzs=zzs,det(a(1:i,1:i); end zzs %(4)求顺序主子式； tzzh=sum(eig(a) %(5)求特征值之和； zdjh=trace(a) %(5)求主对角元之和； tzzj=prod(eig(a) %(5)求特征值之积；,示例2,解：（程序） syms a b c d A=1 a a2 a3;1 b b2 b3;1 c c2 c3;1 d d2 d3; det(A) inv(A),LU分解 QR分解 Cholesky分解 Schur分解,矩阵的分解,LU分解法,矩阵的LU分解又称三角分解，它的目的是将一个矩阵分解成一个下三角矩

5、阵L和一个上三角矩阵U的乘积，即A=LU。只要矩阵非奇异，这种分解是可以进行的。 MATLAB提供lu函数来求矩阵的LU分解，其格式如下： L,U = lu(A) %U为上三角阵，L为下三角阵或其变换形式，满足LU=A。 L,U,P = lu(A) %U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，满足LU=PA。,QR分解法,矩阵的QR分解，它将矩阵分解为一个正交矩阵与一个 上三角矩阵的乘积。,MATLAB提供qr函数来求矩阵的LU分解，其格式如下： Q,R = qr(A) %求得正交矩阵Q和上三角阵R，Q和R满足A=QR。 Q,R,E = qr(A) %求得正交矩阵Q和上三角阵R，E

6、为单位矩阵的变换形式，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时减少了内在数字信息的损耗。,Cholesky分解法,当系数矩阵A正定且对称时,矩阵可分解为A=RTR.。其中R为上三角矩阵。,MATLAB提供chol函数来求矩阵的Cholesky分解，其格式如下： R = chol(X) %如果X为n阶对称正定矩阵，则存在一个实的非

7、奇异上三角阵R，满足R*R = X；若X非正定，则产生错误信息。 R,p = chol(X) %不产生任何错误信息，若X为正定阵，则p=0，R与上相同；若X非正定，则p为正整数，R是有序的上三角阵。,Schur分解法,Schur分解后：AX=b的解为：X=U (T (UTb),MATLAB提供schur函数来求矩阵的Schur分解，其格式如下： T = schur(A) %产生schur矩阵T，即T的主对角线元素为特征值的三角阵。 T = schur(A,flag) %若A有复特征根，则flag=complex，否则flag=real。 U,T = schur(A,) %返回正交矩阵U和sch

8、ur矩阵T，满足A = U*T*U。,稀疏矩阵(sparse matrix)是其元素大部分为零的矩阵。在科学与工程领域中求解线性模型时经常出现大型的稀疏矩阵。在使用计算机存储和操作稀疏矩阵时，经常需要修改标准算法以利用矩阵的稀疏结构。由于其自身的稀疏特性，通过压缩可以大大节省稀疏矩阵的内存代价。更为重要的是，由于过大的尺寸，标准的算法经常无法操作这些稀疏矩阵。 稀疏矩阵的计算速度更快,因为MATLAB只对非零元素进行操作,这是稀疏矩阵的一个突出的优点. 对于矩阵Amn的每个元素aij,知道其行号i和列号j就可以确定其位置.因此对于稀疏矩阵可以用一个结点来存储一个非0元素.该结点可以定义如下:i

9、,j,aij 该结点由3个域组成,i:行号,j:列号;aij元素值.,稀疏矩阵技术,S = sparse(A) %将矩阵A转化为稀疏矩阵形式，即由A的非零元素和下标构成稀疏矩阵S。若A本身为稀疏矩阵，则返回A本身。 S = sparse(m,n) %生成一个mn的所有元素都是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) %生成一个mn的稀疏矩阵，(i,j)对应位置元素为

10、si，m = max(i)且n =max(j)。 S = sparse(i,j,s,m,n,nzmax) %生成一个mn的含有nzmax个非零元素的稀疏矩阵S，nzmax的值必须大于或者等于向量i和j的长度。,稀疏矩阵的创建,函数 full 格式 A=full(S) %S为稀疏矩阵，A为满矩阵。,将稀疏矩阵转化为满矩阵, 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, S=sparse(1:5,1:5,4:8) S = (1,1) 4 (2,2) 5 (3,3) 6 (4,4) 7 (5,5) 8,结果,示例3,

11、示例3：程序,n=10;b=1:n; a1=sparse(1:n,1:n,4,n,n)； a2=sparse(2:n,1:n-1,1,n,n)； a=a1+a2+a2； tic;x=ab;t1=toc，aa=full(a) tic;xx=aab;t2=toc x，xx,注意技巧,计算所用时间,示例4、给定向量组：,求该向量组的秩和一个最大无关组， 并将其余向量用最大无关组线性表出。,向量组的线性相关性分析,a=1,1,2,-1 b=0,2,1,-4 c=1,1,0,-1 d=2,0,3,2 A=a,b,c,d format rat Rref_A=rref(A),Rref_A = 1 0 0 2

12、 0 1 0 -1 0 0 1 0 0 0 0 0,我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断： 若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解； 若系数矩阵的秩rn，则可能有无穷解； 线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。,线性方程组的求解,这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 直接法；另一类是解大型稀疏矩阵 迭代法。 1利用矩阵除法求线性方程组的特解（或一个解） 方程：AX=b 解法：X=Ab 2利用矩

13、阵的LU、QR和cholesky分解求方程组的解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。 （1）LU分解： A*X=b 变成L*U*X=b 所以X=U(Lb) 这样可以大大提高运算速度。 （2）Cholesky分解：若A为对称正定矩阵，方程 A*X=b 变成 R*R*X=b所以X=R(Rb) （3）QR分解 ：方程 A*X=b 变形成 QRX=b 所以 X=R(Qb),求线性方程组的唯一解,示例5,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) %求

14、秩 X=AB %求解 运行后结果如下 R_A = 5 X = 2.2662 -1.7218 1.0571 -0.5940 0.3188,这就是方程组的解,方法一：, 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,方法二：,function kleim(A,b) m,n=size(A); if m=n end X else d

15、isp(Do not work use rule of kleim!) end,克兰姆法则程序命令：( kleim.m),在Matlab中，函数null用来求解零空间，即满足AX=0的解空间，实际上是求出解空间的一组基（基础解系），其格式如下： z = null % z的列向量为方程组的正交规范基，满足Z*Z=I。 z = null(A, r) % z的列向量是方程AX=0的有理基。,求线性齐次方程组的通解,示例6,A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3; format rat %指定有理式格式输出 B=null(A,r) %求解空间的有理基 syms k1 k2 X=k

16、1*B(:,1)+k2*B(:,2) %写出方程组的通解 pretty(X) %让通解表达式更加精美,B= 2 5/3 -2 -4/3 1 0 0 1,X= 2*k1+5/3*k2 -2*k1-4/3*k2 k1 k2,方法一：,A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3; B=rref(A) 运行结果是： B = 1 0 -2 -5/3 0 1 2 4/3 0 0 0 0,方法二：,非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。 因此，步骤为： 第一步：判断AX=b是否有解，若有解则进行第二步 第二步：求AX=b的一个特解 第三步：求AX=0的通解 第四步：

17、AX=b的通解= AX=0的通解+AX=b的一个特解。,求线性非齐次方程组的通解,判断下面的线性方程组是否有解，若有解求其通解。,示例7,jtx_matrix = 1.0000 0 -1.5000 0.7500 1.2500 0 1.0000 -1.5000 -1.7500 -0.2500 0 0 0 0 0,变量的最后一列是方程组的一个特解，由第3、4列可得相应齐次方程组的一个基础解系。因此方程组的通解为,解：A=1,1,-3,-1;3,-1,-3,4;1,5,-9,-8; b=1;4;0 B=A b; jtx_matrix=rref(B),方法一：,A=1,1,-3,-1;3,-1,-3,

18、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-3 1 -1;-3 -1 1 ; Lamda=eig(A); U,T=schur(A),称函数f(x)=ax+b为仿射映射；特别地，当b=0时，称其为线性映射。仿射变换实际上就是经线性变换后再做一次同维的“平移”变换。 仿射变换可以通过一系列的原子变换的复合来实现，包括：平移（Translation）、缩放（Scale）、翻转（Flip）、旋转（Rotation）和剪切（Shear）。,空间中的仿射变换,二维空间中的线性变换,在二维空间中

19、，线性变换可以用 22 的变换矩阵表示。,（一）旋转变换 绕原点逆时针旋转度角的变换公式是 x = xcos ysin 与 y = xsin+ ycos，用矩阵表示为：,（二）反射变换 为了沿经过原点的直线反射向量，假设 (ux, uy) 为直线方向的单位向量。变换矩阵为： 按照不经过原点的直线的反射是仿射变换，而不是线性变换。 特例：关于x轴对称、关于y轴对称、关于直线y=x对称。,二维空间中的线性变换,（三）缩放变换 缩放公式为x=sxx与y=syy，用矩阵表示为： 当sx=sy时，就是以原点为位似中心的位似变换。 当sxsy时，这就是我们平时所说的拉伸变换。,二维空间中的线性变换,（四）

20、剪切变换 切变有两种可能的形式，平行于 x 轴的切变为 x = x + ky 与 y = y，矩阵表示为： 平行于 y 轴的切变为 x = x 与 y = y + kx，矩阵表示为： “剪切变换”又称“错切变换”，指的是类似于四边形不稳定性那种性质，街边小商店那种铁拉门都见过吧？想象一下上面铁条构成的菱形拉动的过程，那就是“错切”的过程。,二维空间中的线性变换,（五）投影变换 为了将向量正投影到一条经过原点的直线，假设 (ux, uy) 是直线方向的单位向量，变换矩阵为： 同反射一样，正投影到一条不经过原点的直线的变换是仿射变换，而不是线性变换。但是透视投影不是线性变换 平行投影也是线性变换，

21、也可以用矩阵表示。 特例：投影到x轴、投影到y轴,二维空间中的线性变换,示例9,已知矩阵,按要求完成下面的实验任务： （1）选取下面不同的变换矩阵T，画出数组w经变换后的轨迹。,clear;clc;clf; w=1,1.5,2,2.5,3,3.5,4,4.5,5;2,1,0,1,2,1,0,1,2; T1=0 1;1 0;T2=1 2;0 1;T3=1 0;2 1; T4=1 0;0 0;T5=0 0;0 1; T6=1/2 sqrt(3)/2;-sqrt(3)/2 1/2; w1=T1*w;w2=T2*w;w3=T3*w;w4=T4*w; w5=T5*w;w6=T6*w; WW=w;w1;w

22、2;w3;w4;w5;w6;hold on for i=1:7 subplot(3,3,i),plot(WW(2*i-1,:),WW(2*i,:) end hold off,示例9：程序,示例9,空间中的正交变换,示例10,按要求完成下面的实验任务： （1）随机生成一个三阶矩阵，判断其是否可逆，若可逆将其列向量组正交化，作成一个正交矩阵； （2）任意选取一个三维非零向量x作为初始迭代向量，用（1）产生的正交矩阵做变换矩阵T，做迭代xk+1=Txk，k=0,1,2, 将迭代生成的向量序列 表示在同一个坐标系中，你能知道他们的分布状态吗？请自己编程检验你的想法。 （3）编程计算向量序列 相应的长度

23、序列 ，检验正交变换的保距性。,clear;clc;clf; a=rand(3) %产生3阶随机阵; a_hls=det(a) %判断a是否可逆; if det(a)=0 T=orth(a) %将a正交化; else disp 矩阵a不可正交化 end,示例10：程序（1）,x=rand(3,1) %产生一个列向量; tx=x; xlcd=norm(x); n=100; for k=1:n x=T*x; tx=tx,x; xlcd=xlcd,norm(x); end plot3(tx(1,:),tx(2,:),tx(3,:),r*) hold on plot3(tx(1,:),tx(2,:),

24、tx(3,:),b-) xlcd,示例10：程序（2、3）,示例10-2,迭代法,迭代法是将求一组解转换为求一个近似解序列的过程，并用最终的近似解来逼近真实解。迭代法要考虑以下3个重要的问题： （1）迭代的初始值 （2）迭代算法 （3）迭代的收敛性,逐次逼近法,对于n阶线性方程组Ax=b的系数矩阵A（假设A是非奇异的）进行如下分解： A=Q-C 其中Q为非奇异的。 则方程组可变换为： x=Bx+r 其中B=Q-1C，r=Q-1b。 这种迭代方法称为逐次逼近法。取不同的Q、C得到不同的迭代算法。通常该方法收敛的充要条件是：迭代矩阵的谱半径小于1。,里查森(Richason)迭代法,里查森迭代法可

25、以说是最简单的迭代法。,迭代法公式为：,对于线性方程组Ax=b，如果A为非奇异方阵，即aii0(i=1,2,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为： x=D-1(L+U)x+D-1b 与之对应的迭代公式为： xk+1=D-1(L+U)xk+D-1b 这就是Jacobi迭代公式。如果序列xk+1收敛于x，则x必是方程Ax=b的解。,Jacobi迭代法,function y,n=jacobi(A,b,x0,eps) if nargin=3 eps=1.0e-6; elseif nargin=eps x0=y;y=B*x0+f;n=n+1; end,Jacobi迭代法的MATLAB函数文件,在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Serdel迭代公式。和Jacobi迭代