孤立特征值情况的矩阵摄动法算例(matlab编程)
- 格式:doc
- 大小:365.20 KB
- 文档页数:8
matlab中eig函数用法
Matlab的eig函数是一个多功能的命令,用于计算矩阵的特征值和特征向量。
它可以被用来求解一些振动类型的简小问题,解决斯坦纳网络的传播问题,判断系统的稳定性以及识别线性系统的状态。
eig函数是用来计算矩阵特征值和特征向量的函数。
简单地说,特征值是一个
复数,表示变换向量空间中的某一矢量。
特征向量是一个列向量,表示变换向量空间中的某一个方向。
eig函数可以求解一个n阶实矩阵A的n个特征值和n个特征
向量。
使用函数的方法如下:
[V,D] = eig(A)
其中A为一个实矩阵,V是A的特征向量矩阵,D是包含A的特征值的对角矩阵。
它们之间的关系是:
AV=VD
其中V是A的特征向量矩阵,D是A的特征值矩阵。
eig函数是一个有效的求解矩阵特征值和特征向量的函数,可以解决一些系统
分析和振动传播类型的问题,同时也可以用它来识别线性系统的状态、判断系统的稳定性等。
matlab 约束矩阵正定
约束矩阵正定性在优化问题中是非常重要的,因为正定矩阵具有很多良好的性质,能够保证优化问题的最优解存在且唯一。
在MATLAB 中,判断一个矩阵是否正定可以使用eig函数,该函数计算矩阵的特征值和特征向量。
如果一个矩阵的所有特征值均大于零,则该矩阵为正定矩阵。
因此,我们可以通过eig函数计算出约束矩阵的特征值,判断它们是否全部大于零,从而判断约束矩阵是否正定。
下面是一个简单的MATLAB代码示例:
C = [1 2 3; 2 5 6; 3 6 9]; % 定义一个矩阵C
[eigVec, eigVal] = eig(C); % 计算矩阵C的特征值和特征向
量
if all(diag(eigVal) > 0) % 判断特征值是否全部大于零
disp('矩阵C为正定矩阵'); % 如果全部大于零,则矩阵C为正定矩阵
else
disp('矩阵C不是正定矩阵'); % 如果存在小于等于零的特征值,则矩阵C不是正定矩阵
end
通过这个简单的代码示例,我们可以看到MATLAB如何判断约束
矩阵的正定性。
在优化问题中,正确地判断约束矩阵是否正定,可以帮助我们更好地处理优化问题,得到更优的解。
matlab特征分解一、Matlab特征分解简介特征分解是矩阵理论中的一个重要概念,可以将一个方阵分解为对角矩阵和相应的特征向量矩阵。
在Matlab中,可以使用eig函数进行特征分解。
二、Matlab特征分解函数eig的使用方法1.基本语法:[V,D] = eig(A),其中A为待分解的方阵,V为特征向量矩阵,D为对角矩阵。
2.示例代码:A = [1 2 3; 4 5 6; 7 8 9]; %定义一个3*3的方阵[V,D] = eig(A); %进行特征分解disp(V); %输出特征向量矩阵disp(D); %输出对角矩阵三、Matlab特征分解的应用场景1.求解线性微分方程组;2.求解差分方程组;3.求解最小二乘问题;4.求解最大化或最小化目标函数等。
四、Matlab特征值和特征向量的性质1.若A是实数矩阵,则其复共轭转置也是其伴随矩阵,即A* = A';2.A与其伴随矩阵A'之积是一个Hermite(共轭对称)矩阵,即AA' = A'A;3.特征向量矩阵V的列向量是相互正交的,即V'V = I;4.特征向量矩阵V的逆矩阵等于其转置矩阵,即V^-1 = V'。
五、Matlab特征分解的相关函数1.eig函数:用于计算方阵的特征值和特征向量;2.svd函数:用于计算奇异值分解(SVD);3.qr函数:用于计算QR分解。
六、Matlab特征分解的注意事项1.在进行特征分解时,需要注意输入的方阵必须是方阵,否则会出现错误;2.在使用eig函数进行特征分解时,需要注意输出结果中特征向量矩阵和对角矩阵并不是按顺序排列的,需要手动排序后才能正确使用。
七、Matlab特征分解实例下面给出一个简单的实例来说明Matlab中如何使用eig函数进行特征分解:A = [4 2; 1 3]; %定义一个2*2的方阵[V,D] = eig(A); %进行特征分解disp(V); %输出特征向量矩阵disp(D); %输出对角矩阵运行以上代码后,输出结果为:V =0.8944 0.70710.4472 -0.7071D =2 00 5其中,V为特征向量矩阵,D为对角矩阵。
关于matlab中的eig函数(求特征值和特征向量)在MATLAB中,eig用途:Find eigenvalues(特征值)and eigenvectors(特征向量),常用的调用格式有5种:(1) E=eig(A):求矩阵A的全部特征值,构成向量E。
(注意,第一列为对应第一个特征值的特征向量)(2) [V,D]=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的全部列向量。
(3) [V,D]=eig(A,'nobalance'):与第2种格式类似,但第2种格式中先对A 作相似变换后求矩阵A的特征值和特征向量,而格式3直接求矩阵A的特征值和特征向量。
(4) E=eig(A,B):由eig(A,B)返回N×N阶方阵A和B的N个广义特征值,构成向量E。
(5) [V,D]=eig(A,B):由eig(A,B)返回方阵A和B的N个广义特征值,构成N×N 阶对角阵D,其对角线上的N个元素即为相应的广义特征值,同时将返回相应的特征向量构成N×N阶满秩矩阵,且满足AV=BVD。
Syntax(句法)如下:d = eig(A)d = eig(A,B)[V,D] = eig(A)[V,D] = eig(A,'nobalance')[V,D] = eig(A,B)[V,D] = eig(A,B,flag)d = eig(A)和 [V,D] = eig(A)最为常用注意,第一列为对应第一个特征值的特征向量,比如:B=rand(4)B =0.5653 0.7883 0.1365 0.97490.2034 0.5579 0.3574 0.65790.5070 0.1541 0.9648 0.08330.5373 0.7229 0.3223 0.3344>> [a,b]=eig(B) %求矩阵B的全部特征值,构成对角阵b,并求B的特征向量构成a的列向量。
关于matlab中的eig函数(求特征值和特征向量)在MATLABxx,eig用途:Find eigenvalues(特征值)and eigenvectors(特征向量),常用的调用格式有5种:(1)E=eig(A):求矩阵A的全部特征值,构成向量E。
(注意,第一列为对应第一个特征值的特征向量)(2)[V,D]=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的全部列向量。
(3) [V,D]=eig(A,'nobalance'):与第2种格式类似,但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量,而格式3直接求矩阵A的特征值和特征向量。
(4)E=eig(A,B):由eig(A,B)返回N×N阶方阵A和B的N个广义特征值,构成向量E。
(5)[V,D]=eig(A,B):由eig(A,B)返回方阵A和B的N个广义特征值,构成N×N阶对角阵D,其对角线上的N个元素即为相应的广义特征值,同时将返回相应的特征向量构成N×N阶满秩矩阵,且满足AV=BVD。
Syntax(句法)如下:d = eig(A)d = eig(A,B)[V,D] = eig(A)[V,D] = eig(A,'nobalance')[V,D] = eig(A,B)[V,D] = eig(A,B,flag)d = eig(A)和[V,D] = eig(A)最为常用注意,第一列为对应第一个特征值的特征向量,比如:B=rand(4)B =0.5653 0.7883 0.1365 0.97490.2034 0.5579 0.3574 0.65790.5070 0.1541 0.9648 0.08330.5373 0.7229 0.3223 0.3344>> [a,b]=eig(B) %求矩阵B的全部特征值,构成对角阵b,并求B的特征向量构成a的列向量。
矩阵方程的基础解系是指由矩阵方程的特解和齐次方程的通解组成的一种特殊解系,它包括了矩阵方程的所有解。
在 MATLAB 中,我们可以利用一些基本的函数和方法来求解矩阵方程的基础解系。
接下来我们将介绍如何使用 MATLAB 求解矩阵方程的基础解系的基础方法。
1. 准备工作在使用 MATLAB 求解矩阵方程的基础解系之前,首先需要确保MATLAB 环境已经安装并且处于可用状态。
需要定义矩阵系数和常数向量。
2. 使用 `\` 求解特解在 MATLAB 中,可以使用 `\` 运算符来求解线性方程组的特解。
具体而言,如果矩阵方程为 `A * X = B`,那么可以使用 `X = A \ B` 来求解特解。
3. 求解齐次方程的通解同样地,在 MATLAB 中,可以使用 `null` 函数来求解齐次方程的通解。
具体而言,如果矩阵方程为 `A * X = 0`,那么可以使用 `N =null(A,'r')` 来求解齐次方程的通解。
4. 合并特解和通解将特解和通解合并起来,就得到了矩阵方程的基础解系。
具体而言,可以使用 `X = special + nullspace` 来完成合并操作。
总结通过以上步骤,我们可以使用 MATLAB 求解矩阵方程的基础解系。
首先求解特解,然后求解齐次方程的通解,最后合并特解和通解即可得到基础解系。
这样的方法可以帮助我们更好地理解和求解矩阵方程的基础解系,同时也展现了 MATLAB 在数学计算方面的强大功能和灵活性。
在实际问题中,矩阵方程的基础解系往往是非常重要的,它可以帮助我们求解各种复杂的线性方程组和矩阵方程,从而进一步应用到实际的工程和科学问题中。
掌握如何使用 MATLAB 求解矩阵方程的基础解系是非常有必要的。
希望本文的介绍能够对读者有所帮助,也欢迎大家多多交流和共享关于 MATLAB 求解矩阵方程的经验和技巧。
在实际工作和学习中,矩阵方程的基础解系是一个非常重要的概念。
基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例杨海涛【摘要】基于Matlab的数学实验,目的是为了提高普通高校学生的数学兴趣和动手能力,为培养解决实际数学问题的能力而设计的实验方法和步骤.利用矩阵的奇异值分解svd函数的数值算法设计了一种在实验课上用svd算法仿真的实验方案.【期刊名称】《内蒙古民族大学学报(自然科学版)》【年(卷),期】2016(031)006【总页数】4页(P465-468)【关键词】数学实验;数值计算;Matlab;奇异值分解【作者】杨海涛【作者单位】内蒙古民族大学数学学院,内蒙古通辽 028043【正文语种】中文【中图分类】O241.4利用计算机计算解决实际问题时需要建立数学模型,同时也需要相应数学软件的支持,为了更好的掌握数学软件这一有利的工具,本文通过Matlab中的矩阵的奇异值分解〔1〕svd函数设计了如何学习和掌握该函数的实验方法,从而可以进一步学习其它函数.目前,Matlab在计算仿真和图形图像的压缩处理〔2~4〕有着广泛的应用,希望这个实验方法能够给正在学习数学软件的学生起到抛砖引玉的作用.1.1 理论基础任意复长方矩阵的奇异值分解定理如下(Autonne-Eckart-Young定理):〔5〕设A为一个m×n的实矩阵(m≥n),则存在正交矩阵V和U,使得:其中∑=diag(σ1,σ2,…,σr),而且σ1≥σ2≥…≥σr≥0.奇异值分解在统计分析、信号与图像处理、系统理论和控制中有着广泛的应用. 1.2 SVD的数值算法〔6〕求解一般矩阵奇异值问题最常用的算法可分为两大类:一类是QR分解的方法;另一类是Jacobi旋转的方法.奇异值分解的QR分解法分为两个阶段:一是矩阵的二重对角化,即利用Householder变换把矩阵A变换成二重对角矩阵.这个过程就是追赶(chasing)过程〔7〕.二是QR分解,即保持二重对角矩阵形式不变,利用正交变换使对角线元素逐渐减小,使得矩阵接近对角矩阵.1.2.1 svd函数.对于matlab中的svd函数不是开源的,但通过调研得知该函数的算法就是采用propack中的lansvd函数的算法.其算法过程如下:对于任意的矩阵Am×n可以利用lanczos过程将其双对角化,即存在m阶正交矩阵Q和n阶正交矩阵W满足:其中为一双对角阵,且αi, βi∈R (i=1, 2, …, n),这里假设了m>n.上面的分解是逐步的利用Householder左右变换得到的.算法中利用了已知函数完成的双对角化过程.这里如果我们已经知道Bn的svd分解为其中Vn+1为n+1阶正交矩阵,的广义对角矩阵Un为n阶正交矩阵.则有令V=Qn+1Vn+1,其前r列为对应的左奇异向量,令U=WUn,其前r列为对应的右奇异向量.因此算法的关键在于计算矩阵Bn的svd分解.因为Bn是一个双对角矩阵,它的奇异值分解计算比较简单,运算量较小,利用Givens变换,算法复杂度仅为O(n).在我们这个程序中,直接调用了库函数svd 计算Bn的奇异值.1.2.2 svds函数.对于matlab中的svds函数,算法流程如下:(1)对于m×n阶矩阵A,构造(m+n)×(m+n)阶矩阵,其中O是稀疏矩阵. (2)利用eigs函数计算出B矩阵的特征值矩阵D,以及D对应的正交阵W. (3)根据输入参数,按要求输出最大特征值、最接近某参数的特征值或对特征值进行排序等.(4)根据设定的阀值,对矩阵进行筛选得到U矩阵,对矩阵筛选得到V矩阵,对D矩阵进行筛选得到S矩阵.(5)根据S中元素的大小对U、V、S进行排序.(6)输出结果.有如下J-W定理(Jordan-Wielandt定理):若σ1≥σ2≥…≥σp-1≥σp是矩阵Am×n的奇异值(其中,p=min{m,n}),则上述B矩阵具有2p个非零特征值-σ1,…,-σp,σp,…,σ1和|m-n|个零特征值,与非零特征值(±σj)相对应的特征向量为若m≠n,另有特征向量分别对应于m>n或m<n.svds函数的基本思想就是构造出这样的矩阵B,并且求出其特征值和特征向量,然后根据函数重构的条件(所需奇异值的数量、最大最小等),求出所需的奇异值和对应的向量.2.1 矩阵生成生成10个不同的m×n的矩阵,这里取m=120 , n=50,对这10个矩阵进行svd分解.矩阵生成的方法是先做出矩阵的奇异值,然后用QR分解做出10个120阶正交阵和10个50阶正交阵,再利用奇异值分解公式生成10个矩阵A.10个矩阵的奇异值(均为25个)如表1所示.2.2 仿真结果2.2.1 计算时间svd和svds两种方法对十个矩阵进行奇异值分解消耗的时间如表2,单位是秒.可见,svd函数运行时间一般要小于svds函数运行时间.2.2.2 最小奇异值计算准确率.通过对比最小特征值的计算准确度,将每种算法特征值的最小值与矩阵最开始生成时的最小特征值进行比较.比较结果两种方法对10个矩阵的svd分解运算最小特征值的误差均几乎为0,都在10-15次方的量级左右,见表3.其中delta_min_svd和delta_min_svds是matlab程序中两个变量值.这个误差已经和截断误差非常接近了,因此可以认为两种方法都能较好的进行奇异值计算.2.2.3 左右奇异值向量计算准确率.由于svd分解后U、V矩阵的不唯一性,这里我们只看生成的U、V的正交性是不是依然比较完好.通过计算每一列和它的转置相乘减去单位阵也都在10-15这样的量级,因此svd和svds对于正交性的保持都是非常好的.综上所述,通过上述设计的方法和步骤,完全可以掌握matlab中的M-文件和.mat文件的使用方法,并且能熟练掌握奇异值分解的理论知识和具体计算方法,才能为以后svd应用打下基础,这样的实验方法才有实际的意义,才是真正的综合设计实验.本文的宗旨在于引导学生如何用数学软件完成一个实验设计,同时掌握这个实验的目的、相关内容实验的方法和手段.〔1〕尹芳黎,杨雁莹,王传栋,等.矩阵奇异值分解及其在高维数据处理中的应用〔J〕.数学的实践与认识,2011,45(15): 171-176.〔2〕胡乡峰,卫金茂.基于奇异值分解(SVD)的图像压缩〔J〕.东北师大学报:自然科学版,2006,38(3):36-39.〔3〕曾超,张陈,徐永利.基于奇异值分解的图像压缩及其matlab实现〔J〕.科技信息,2010,14:484.〔4〕向培素.一种自适应AP算法的matlab实现〔J〕.西南民族大学学报:自然科学版,2014,40(16):877-882.〔5〕张贤达.矩阵分析与应用〔M〕.北京:清华大学出版社,2004.358-360. 〔6〕姜健飞,胡良剑,唐俭.数值分析及其matlab实验〔M〕.北京:科学出版社,2004.6,219-224.〔7〕山晓东,李园.一类不可约L-矩阵的预条件AOR迭代法〔J〕.湖北民族学院学报:自然科学版,2014,32(2):128-132.。
《应用计算方法教程》m a t l a b作业二-CAL-FENGHAI.-(YICAI)-Company One1作业六6-1 试验目的 计算特征值,实现算法试验内容:随机产生一个10阶整数矩阵,各数均在-5和5之间。
(1) 用MATLAB 函数“eig ”求矩阵全部特征值。
(2) 用幂法求A 的主特征值及对应的特征向量。
(3) 用基本QR 算法求全部特征值(可用MATLAB 函数“qr ”实现矩阵的QR 分解)。
原理幂法:设矩阵A 的特征值为12n ||>||||λλλ≥⋅⋅⋅≥并设A 有完全的特征向量系12,,,n χχχ⋅⋅⋅(它们线性无关),则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim()k j k k jV V λ→∞-=,其中()k j V 表示向量的第j 个分量。
为避免逐次迭代向量k V 不为零的分量变得很大(1||1λ> 时)或很小(1||1λ< 时),将每一步的k V 按其模最大的元素进行归一化。
具体过程如下:选择初始向量0V ,令1max(),,,1kk k k k k kV m V U V AU k m +===≥,当k 充分大时1111,max()max()k k U V χλχ+≈≈。
QR 法求全部特征值:11111222111,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==⋅⎧⎪⋅==⋅⎪=⋅⋅⋅⎨⋅⋅⋅⎪⎪⋅==⋅⎩ 由于此题的矩阵是10阶的,上述算法计算时间过长,考虑采用改进算法——移位加速。
迭代格式如下:1k k k kk k k k A q I Q R A R Q q I +-=⎧⎨=+⎩ 计算k A 右下角的二阶矩阵()()1,11,()(),1,k k nn n n k k n n n n a a a a ----⎛⎫ ⎪ ⎪⎝⎭的特征值()()1,k k n n λλ-,当()()1,k k n n λλ-为实数时,选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。
matlab的arnoldi方法
Arnoldi方法是一种用于求解线性方程组的迭代算法,也可以用于计算矩阵的奇异值分解(SVD)。
在MATLAB中,Arnoldi方法可以通过内置函数
`arnoldi`实现。
以下是使用Arnoldi方法的基本步骤:
1. 定义矩阵A和向量b,其中A是要求解的线性方程组的系数矩阵,b是方程组的常数项向量。
2. 调用`arnoldi`函数,将A和b作为输入参数传递给函数。
3. `arnoldi`函数将返回一个包含Arnoldi分解结果的对象,其中包含三个部分:
`H`: 包含Arnoldi向量的矩阵,H的第i列是对应的Arnoldi向量。
`Q`: 正交矩阵,将原向量空间映射到包含Arnoldi向量的子空间。
`[V,j]`: 包含一组归一化列的矩阵V和对应的指标j,这些列是原矩阵A在子空间中的投影。
4. 可以使用返回的Arnoldi分解结果来求解线性方程组或计算矩阵的奇异值分解。
例如,可以使用返回的Q和H矩阵来求解Ax = b,计算矩阵A的奇异值分解等。
需要注意的是,Arnoldi方法是一种迭代算法,需要指定迭代的次数或收敛条件。
在MATLAB中,可以通过设置`arnoldi`函数的参数来指定迭代的次数或收敛条件。
matlab eigs函数用法Matlab是一个非常强大的数值计算软件,拥有丰富的函数库,eigs函数就是其中之一。
eigs函数是用于计算矩阵的特征值和特征向量的函数,并且它还可以计算大规模矩阵的部分特征值和特征向量。
在这篇文章中,我们将详细介绍eigs函数的用法。
我们将从基本的语法开始,逐步深入到其功能的各个方面。
希望这篇文章能够帮助那些对eigs函数感兴趣的读者更好地理解和使用它。
一、基本语法eigs函数的基本语法如下:[V,D] = eigs(A):计算矩阵A的全部特征值和特征向量,并将其分别存储在矩阵D和矩阵V中。
[V,D] = eigs(A,k):计算矩阵A的最大(或最小)k个特征值和特征向量。
如果不指定k的值,默认为6。
[V,D] = eigs(A,k,sigma):计算矩阵A的与指定值sigma最接近的k个特征值和特征向量。
[V,D] = eigs(A,B):计算广义特征值问题的特征值和特征向量,其中A和B是正方形矩阵。
[V,D,flag] = eigs(___):返回计算信息的标识符,用于指示计算是否成功。
二、计算矩阵的全部特征值和特征向量我们首先从计算矩阵的全部特征值和特征向量开始使用eigs函数。
假设我们有一个3x3的矩阵A,代码如下:A = [1 2 3; 4 5 6; 7 8 9];我们可以使用以下代码计算矩阵A的全部特征值和特征向量:[V,D] = eigs(A);这里,V和D分别是矩阵A的特征向量和特征值,可以使用disp函数来显示它们:disp(V);disp(D);结果如下:V =0.3868 0.5086 -0.44720.5820 0.2634 -0.63250.7771 -0.9820 0.5270D =16.1168 0 00 -1.1168 00 0 -0.0000从结果中可以看出,V是一个3x3的矩阵,其中每列是矩阵A的一个特征向量。
D是一个3x3的对角矩阵,其中对角线上的元素是矩阵A的特征值。
matlab特征值分解和奇异值分解特征值分解函数 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矩阵。
说明⼀般特征值问题是求解⽅程:解的问题。
⼴义特征值问题是求⽅程:解的问题。
奇异值分解函数 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。
奇异值分解压缩图像clear all;close all;clc;a=imread('C:\Users\ranji\Desktop\rgb_image.jpg');imshow(mat2gray(a))[m n]=size(a);a=double(a);%r=rank(a);[s v d]=svd(a(:,:,1)); %取⼀个分量%re=s*v*d';re=s(:,:)*v(:,1:1)*d(:,1:1)';figure;imshow(mat2gray(re));imwrite(mat2gray(re),'C:\Users\ranji\Desktop\1.jpg')re1=s(:,:)*v(:,1:20)*d(:,1:20)';figure;imshow(mat2gray(re));imwrite(mat2gray(re),'C:\Users\ranji\Desktop\2.jpg')re=s(:,:)*v(:,1:80)*d(:,1:80)';figure;imshow(mat2gray(re));imwrite(mat2gray(re),'C:\Users\ranji\Desktop\3.jpg')re=s(:,:)*v(:,1:150)*d(:,1:150)';figure;imshow(mat2gray(re));imwrite(mat2gray(re),'C:\Users\ranji\Desktop\4.jpg')不同特征值进⾏重构的效果。
1.% 取指定特征值对应的特征向量2.clc; clear all; close all;3. A = [1 3 74. 3 8 95. 5 4 6];6.k = 2;7.[V, D] = eig(A);8.% 返回的V是以列向量对应的特征向量9.% D是对角线上为特征值的矩阵10.D = diag(D);11.[D, I] = sort(D, 'descend');12.% 得到了对应的排序13.if k > length(D)14. k = length(D);15.end16.temp(1) = {D(1 : k)};17.temp(2) = {V(:, I(1 : k))};18.fprintf('\n 前k个最大特征值 : \n');19.celldisp(temp(1));20.fprintf('\n 前k个最大特征值对应的特征向量 : \n');21.celldisp(temp(2));[v,d]=eig(A);A为你的矩阵,V为特征向量矩阵,D为特征值矩阵,然后对D求最大值即可得最大特征根![V,D] = EIG(X) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that X*V = V*D.V是特征向量,D是特征值实例:矩阵:1 2/3 7/3 7/33/2 1 3/2 3/23/7 2/3 1 3/23/7 2/3 2/3 1>> format rat>> A=[1 2/3 7/3 7/33/2 1 3/2 3/23/7 2/3 1 3/23/7 2/3 2/3 1]A =1 2/3 7/3 7/33/2 1 3/2 3/23/7 2/3 1 3/23/7 2/3 2/3 1>> [V,D]=eig(A)V =1793/2855 504/3235 - 146/235i 504/3235 + 146/235i 1990/4773670/1079 -3527/5220 -3527/5220 -509/9594350/11989 1160/4499 + 287/3868i 1160/4499 - 287/3868i -350/647838/2819 181/3874 + 1179/4852i 181/3874 - 1179/4852i 1238/2467D =810/197 0 0 00 -93/4229 + 455/674i 0 00 0 -93/4229 - 455/674i 00 0 0 -149/2201******************************************************************************* **********如何归一化求权重呢?>> a=[1 3 5;1/3 1 3; 1/5 1/3 1]a =1.0000 3.0000 5.00000.3333 1.0000 3.00000.2000 0.3333 1.0000>> [V,D]=eig(a)V =0.9161 0.9161 0.91610.3715 -0.1857 + 0.3217i -0.1857 - 0.3217i0.1506 -0.0753 - 0.1304i -0.0753 + 0.1304iD =3.0385 0 00 -0.0193 + 0.3415i 00 0 -0.0193 - 0.3415i**************************************************************************>> a=[1 2 4 8 6 6 8;1/2 1 2 6 4 4 8;1/4 1/2 1 4 2 4 6;1/8 1/6 1/4 1 2 2 4;1/6 1/4 1/2 1/2 1 1 4;1/6 1/4 1/4 1/2 1 1 2;1/8 1/8 1/6 1/4 1/4 1/2 1]a =1.00002.0000 4.0000 8.0000 6.0000 6.0000 8.0000 0.5000 1.0000 2.0000 6.0000 4.0000 4.0000 8.0000 0.2500 0.5000 1.0000 4.0000 2.0000 4.0000 6.0000 0.1250 0.1667 0.2500 1.0000 2.0000 2.0000 4.0000 0.1667 0.2500 0.5000 0.5000 1.0000 1.0000 4.0000 0.1667 0.2500 0.2500 0.5000 1.0000 1.0000 2.0000 0.1250 0.1250 0.1667 0.2500 0.2500 0.5000 1.0000>> rats(a)ans =1 2 4 8 6 6 81/2 1 2 6 4 4 81/4 1/2 1 4 2 4 61/8 1/6 1/4 1 2 2 41/6 1/4 1/2 1/2 1 1 41/6 1/4 1/4 1/2 1 1 21/8 1/8 1/6 1/4 1/4 1/2 1>> [V,D]=eig(a)V =0.7884 0.8327 0.8327 0.8083 0.8083 -0.5119 + 0.3865i -0.5119 - 0.3865i0.4894 0.3216 + 0.2636i 0.3216 - 0.2636i 0.1760 +0.0792i 0.1760 - 0.0792i 0.6783 0.67830.3038 0.0883 + 0.2728i 0.0883 - 0.2728i -0.4630 +0.1038i -0.4630 - 0.1038i -0.2011 - 0.2400i -0.2011 + 0.2400i0.1404 -0.1620 + 0.1018i -0.1620 - 0.1018i 0.0620 -0.0510i 0.0620 + 0.0510i -0.0006 + 0.1021i -0.0006 - 0.1021i0.1215 -0.0627 - 0.0658i -0.0627 + 0.0658i 0.0367 +0.2360i 0.0367 - 0.2360i -0.0531 + 0.0357i -0.0531 - 0.0357i0.0975 -0.0303 - 0.0476i -0.0303 + 0.0476i 0.0488 -0.1148i 0.0488 + 0.1148i 0.0231 - 0.1221i 0.0231 + 0.1221i0.0508 0.0030 - 0.0590i 0.0030 + 0.0590i -0.0561 -0.0454i -0.0561 + 0.0454i 0.0102 + 0.0197i 0.0102 - 0.0197iD =7.3899 0 0 0 0 0 00 -0.0008 +1.5369i 0 0 0 0 00 0 -0.0008 -1.5369i 0 0 0 00 0 0 -0.1624 +0.6552i 0 0 00 0 0 0 -0.1624 - 0.6552i 0 00 0 0 0 0 -0.0317 + 0.2040i 00 0 0 0 0 0 -0.0317 - 0.2040i。
16. 数值计算—线代篇一、行列式det(A)——矩阵A的行列式;inv(A)——矩阵A的逆;rank(A)——矩阵A的秩;B(: , i)=b——将向量b赋给矩阵B的第i行;[A, eye(5)]——在矩阵A右端,拼接5阶单位矩阵;[U,s]=rref(A)——对矩阵A作行变换,U返回A的最简行阶梯形矩阵,s为行向量存储U的各行首个非0元所在列号,length(s)即为A的秩;例1 用初等行变换法求矩阵123221343A⎛⎫⎪= ⎪⎪⎝⎭的逆。
代码:format short g % 省略小数位多余的0A=[1 2 3; 2 2 1; 3 4 3];B=rref([A,eye(3)])% 对矩阵[A,I]进行初等行变换,得到最简行阶梯矩阵Bif(rank(B(:,1:3))==3)% 判断B的前3列是否为单位阵,若是取出后3列,即A逆A1=B(:,4:6)elsedisp('A不可逆');end运行结果:B = 1 0 0 1 3 -20 1 0 -1.5 -3 2.50 0 1 1 1 -1A1 = 1 3 -2-1.5 -3 2.51 1 -1例2 解方程2232113221051327132x x -=- 代码: syms x;A=[3 2 1 1;3 2 2-x^2 1;5 1 3 2;7-x^2 1 3 2];D=det(A)f=factor(D) % 对行列式D 进行因式分解X=solve(D) % 求方程“D =0”的解运行结果:D = -3*(x^2 - 1)*(x^2 - 2)f =-3*(x - 1)*(x + 1)*(x^2 - 2)X = -112^(1/2)-2^(1/2)二、向量组的线性相关性例3 向量组123452434713216,,,,313157534170ααααα⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-----⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=====-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦求它的秩和一个最大线性无关组,并用来表示其它向量。
在MATLAB中实现Truncated Power Iteration算法,首先需要了解Truncated Power Iteration是一种用于计算矩阵最大特征值的迭代算法。
这种算法在每次迭代中只保留最大的特征值,而忽略其他较小的特征值。
以下是一个简单的MATLAB实现:matlab复制代码function[eigval, eigvec] = truncated_power_iteration(A, tol,max_iter)A: 输入的矩阵tol: 精度要求max_iter: 最大迭代次数初始化v = rand(size(A, 1), 1);v = v / norm(v); 单位化向量eigval = zeros(max_iter, 1);eigvec = zeros(size(A, 1), max_iter);eigvec(:, 1) = v; 初始化特征向量矩阵for i = 1:max_iterw = A * eigvec(:, i); 计算新的向量eigval(i) = w' * v; 计算特征值if abs(eigval(i) - eigval(i-1)) < tol 检查收敛性break;endw = w / eigval(i); 单位化向量eigvec(:, i+1) = w; 更新特征向量矩阵endend这个函数接受一个矩阵A,一个精度要求tol和一个最大迭代次数max_iter作为输入,并返回计算出的特征值和特征向量。
在函数中,我们首先初始化一个随机向量v,然后开始迭代。
在每次迭代中,我们计算新的向量w,并使用它来更新特征值和特征向量。
我们检查收敛性,如果相邻的两个特征值之间的差小于给定的精度要求,我们就停止迭代。
注意,这个实现可能并不适用于所有情况,特别是当矩阵A的条件数非常大时。
在这种情况下,可能需要使用更复杂的算法或者预处理步骤来改善结果。
设原始特征值问题是 )~1(00000nixMxKiii (1)
相应的正交规范条件是 )~1,(000njixMxijjTi (2)
先考虑一个五个自由度的系统. 其中,0M、0K分别是原系统的nn阶对称质量阵、刚度阵,0i是特征值,且200ii,
0i是固有频率,0ix是相应的特征向量。
011110.5
M
02112112112111K
系统结构参数改变后,相应的质量阵、刚度阵均有相应的变化,设系统修改后的质量阵、刚度阵分别为
1010KKKMMM
(3)
其中取
111110.5
M 110000K
显然,新系统(即结构修改后的系统或称摄动系统)的特征值问题及相应的正交规范条件是
iiixMMxKK)()(1010
(4)
ijjTixMMx)(10 (5) 其中,i为新系统的特征值,ix 是相应的特征向量。 根据摄动理论,可将i、ix按小参数展开成幂级数(因此,胡海昌院士称其为小参数法),即 )(32210iiii (6)
)(32210iiiixxxx (7)
现在来确定1i,2i,1ix,2ix,将(6)、(7)式代入(4)式得 ))()(())((2210102210221010iiiiiiiiixxxMMxxxKK 展开上式,略去)(3,并比较的同次幂系数可得 000000:iiixMxK (8-1)
00101010001101:iiiiiiiixMxMxMxKxK (8-2)
00201110111020011202:iiiiiiiiiiiixMxMxMxMxMxKxK (8-3)
将(6)、(7)式代入(5)式得
ijjjjTiiixxxMMxxx))(()(2210102210 展开上式,并令ji(从后面的分析可得,ji的关系式在推导中没有被直接利用),比较的同次幂系数可得
1:0000iTixMx (9-1)
0:0100011001iTiiTiiTixMxxMxxMx (9-2)
0:0021010111102002iTiiTiiTiiTiiTixMxxMxxMxxMxxMx (9-3)
至此,我们已推得了进行摄动分析(即求解i,ix中的1i,2i,1ix,2ix)所需的全部基本方程。实际上,(8-1)、(9-1)即为(1)、(2)式,是显然满足的。于是,可在原系统的特征值问题的基础上,通过(8-2)、(9-2)求解一阶摄动1i,1ix,通过(8-3)、(9-3)求解二阶
摄动2i,2ix,代回(6)、(7)式,即求得i,ix,且具有)(3精度。
一阶摄动公式 可将1ix表示成原系统模态(向量)02010,,,nxxx的线性组合,即 nkkkixx1011 (10)
其中,1k是n个待定系数。 将(10)代入(8-2),并左乘Tkx0,得
0001010010100001010100iTkiiTkinkkkTkiiTknkkkTkxMxxMxxMxxKxxKx 利用正交规范关系(2)式,上式可简化成
000101000101001iTkiiTkiikiTkkkxMxxMxxKx (11)
ki 时,00ik,1000iTkxMx,由(11)式可得特征值的一阶摄动为
010101)(iiTiixMKx (12)
ki时,0000iTkxMx,则由(11)式可得
)()(00010101kixMKxkiiiTkk
(13)
至此,n个待定系数1k中,只有1i尚未确定,现在来求1i。 用00MxTi左乘(10)式两边,得
110001100inkkTikiTixMxxMx (14)
(14)式转置,并考虑到0M对称,且1i是一个常数,于是得 1001iiTixMx (15)
(14)、(15)式代入(9-2)可得 010121iTiixMx (16)
由(13)、(16)两式即完全确定了)~1(1nkk,也就确定了(10)式的1ix,而1i由(12)式给出,于是可得一阶摄动公式为
0010,1000010101010101)(21)()(iiTinikkkiiiTkiiiTiixxMxxxMKxxxMKx
(17)
二阶摄动公式 为了得到更精确的摄动解,需要用到二阶摄动。根据展开定理,将2ix按02010,,,nxxx展开 nkkkixx1022 (18)
将(18)式代入(8-3),并左乘Tkx0,得
0002011101110010200011010200)(iTkiiiiiiiTknkkkTkiiTknkkkTkxMxxMxMxMxxMxxKxxKx
利用(2)式,上式可简化成 000201110111000211002)(iTkiiiiiiiTkikiTkkkxMxxMxMxMxxKx
(19)
ki 时,00ik,1000iTkxMx,由(19)式可得特征值的二阶摄动为
)(0111011101102iiiiiiiTiixMxMxMxKx (20)
ki时,0000iTkxMx,则由(19)式可得
)()(000111011101102kixMxMxMxKxkiiiiiiiiTkk
(21)
现在来确定最后一个系数2i,为此,用00MxTi左乘(18)式两边,得 210200200inkkkTiiTixMxxMx (22)
(22)式转置后,注意到0M对称,2i是一个数,可得 2002iiTixMx (23)
(22)、(23)式代入(9-3)可得 )(210111011102iTiiTiiTiixMxxMxxMx (24)
于是就得到了二阶摄动公式
0011101110,100001110111011020111011101102)(21)()(iiTiiTiiTinikkkiiiiiiiiTkiiiiiiiiTiixxMxxMxxMxxxMxMxMxKxxxMxMxMxKx
编程如下 M0=diag([1,1,1,1,0.5]); K0=[2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,1]; M1=diag([1,1,1,1,0.5]); K1=diag([-1,0,0,0,0]); [x0,r0]=eig(K0,M0); %%[V,D] = eig(A,B) produces a diagonal matrix D of generalized eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that A*V = B*V*D . %%x0为特征向量矩阵,r0为特征值对角阵 %%下面来计算一阶摄动 for m=1:5 r1(m,m)=x0(:,m)'*(K1-r0(m,m)*M1)*x0(:,m); x1(:,m)=-0.5*x0(:,m)'*M1*x0(:,m)*x0(:,m); for n=1:5 if n~=m
x1(:,m)=x1(:,m)+x0(:,n)'*(K1-r0(m,m)*M1)*x0(:,m)/(r0(m,m)-r0(n,n))*M0(:,n); end end end %%再来计算二阶摄动 for m=1:5 r2(m,m)=x0(:,m)'*(K1*x1(:,m)-r0(m,m)*M1*x1(:,m)-r1(m,m)*M0*x1(:,m)-r1(m,m)*M1*x0(:,m));