矩阵计算-MATLAB-幂法程序
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量-附Matlab程序
矩阵的特征值与特征向量的计算摘要物理,力学,工程技术中的很多问题在数学上都归结于求矩阵特征值的问题,例如振动问题(桥梁的振动,机械的振动,电磁振动等)、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分,本文使用幂法和反幂法分别求矩阵的按模最大,按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法,它最大的优点是方法简单,对于稀疏矩阵比较合适,但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值,及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词:矩阵;特征值;特征向量;冥法;反冥法THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACTPhysics, mechanics, engineering technology in a lot of problems in mathematics are attributed to matrix eigenvalue problem, such as vibration (vibration of the bridge, mechanical vibration, electromagnetic vibration, etc.) in physics, some critical values determine problems and theoretical physics in some of the problems. Matrix eigenvalue calculation is a very important part in matrix computation. In this paper, we use the power method and inverse power method to calculate the maximum of the matrix, according to the minimum characteristic vector and the corresponding characteristic value.Power method is an iterative method to calculate the eigenvalues of a matrix. It has the advantage that the method is simple and suitable for sparse matrices, but sometimes the convergence rate is very slow. The basic idea is to take a non - zero initial vector. Construct a vector sequence from the matrix of the matrix. Then the eigenvalues and eigenvectors are obtained by using the constructed vector sequence.The inverse power method is used to calculate the minimum feature vectorsand their eigenvalues of the matrix, and to calculate the eigenvalues of the matrix. In this paper, we use the inverse power method to calculate the minimum eigenvalue of a matrix and its corresponding eigenvalues. The basic idea of calculating the minimum characteristic vector of a matrix is to transform it to the maximum characteristic vector of the modulus of the inverse matrix. Then, according to the model, the minimum feature vector of the original matrix is introduced.Key words: Matrix;Eigenvalue;Eigenvector;Iteration methods;目录1 引言 (1)2 相关定理。
matlab幂法求特征值与特征向量 -回复
matlab幂法求特征值与特征向量-回复Matlab幂法求特征值与特征向量Matlab是一种常用的数学软件,它提供了一系列强大的数值计算工具和函数,旨在简化数学建模和计算的过程。
其中,求解特征值与特征向量是矩阵分析与线性代数中的重要问题之一。
在此,我们将介绍如何使用Matlab中的幂法来求解矩阵的特征值与特征向量。
特征值与特征向量是矩阵分析的基本概念。
给定一个矩阵A,如果存在一个非零向量x,使得Ax=λx,其中λ是一个实数,则称λ为A的特征值,x 为相应于特征值λ的特征向量。
在Matlab中,计算矩阵的特征值与特征向量可以使用`eig`函数。
这个函数能够计算矩阵所有特征值的值,其中特征值按照降序排列。
对于复杂特征值,这个函数会返回具有相应特征向量的V矩阵。
然而,幂法是一种迭代方法,可用于估计矩阵A的最大特征值λ和相应的特征向量x。
幂法的基本思想是利用矩阵的特征值分解性质中最大特征值的绝对值大于其他特征值的绝对值,从而将问题简化为求解最大特征值及其特征向量。
下面,我们将以以下步骤详细介绍如何使用Matlab中的幂法求解矩阵的特征值与特征向量:步骤1:定义初始向量x0首先定义一个非零的初始向量x0。
该向量可以是随机生成的,或者是具有合理初始值的向量。
步骤2:计算矩阵的迭代利用初始向量x0和矩阵A,计算下一个迭代向量x1。
具体而言,使用x0得到x1通过以下公式计算:x1 = A * x0步骤3:归一化迭代向量计算归一化的迭代向量x1。
这可以通过除以向量中的最大元素来完成。
归一化向量可以确保以后的计算产生可靠结果。
x1 = x1 / max(x1)步骤4:计算特征值估计计算特征值的估计值λ。
这可以通过计算x1的无穷范数与x0的无穷范数之比来实现:λ= norm(x1,Inf) / norm(x0,Inf)步骤5:收敛判断判断计算得到的特征值估计是否收敛。
这可以通过设定一个容差值来实现,在误差满足一定条件时停止迭代计算。
矩阵幂和矩阵指数函数的计算方法
矩阵幂和矩阵指数函数的计算方法矩阵幂和矩阵指数函数是矩阵运算中比较重要的两个概念。
在矩阵幂和矩阵指数函数的计算过程中,我们需要用到一些特殊的算法和方法。
本文将介绍矩阵幂和矩阵指数函数的概念、计算方法和应用等方面的内容,帮助读者更好地了解和掌握这两个概念。
一、矩阵幂的概念对于一个$n$阶矩阵$A$,设$k$为一个自然数,则$A^k$表示$k$次幂。
即:$A^k=\underbrace{A \times A \times \cdots \times A}_{k\text{个} A}$其中,当$k=0$时,$A^k$等于$n$阶单位矩阵$I_n$。
矩阵幂的计算过程中,我们需要用到矩阵乘法的定义。
对于两个$n$阶矩阵$A$和$B$,它们的乘积$AB$定义为:$AB=[c_{ij}]=\sum_{k=1}^na_{ik}b_{kj}$其中,$c_{ij}$表示矩阵的第$i$行第$j$列的元素,$a_{ik}$和$b_{kj}$分别表示第$i$行第$k$列的元素和第$k$行第$j$列的元素。
二、矩阵幂的计算方法矩阵幂的计算方法有两种:直接幂法和快速幂法。
1. 直接幂法直接幂法是一种比较简单的计算矩阵幂的方法。
对于一个$n$阶矩阵$A$和一个自然数$k$,我们可以通过$k-1$次连乘的方式计算出$A^k$的值。
即:$A^k=\underbrace{A \times A \times \cdots \times A}_{k-1\text{个} A} \times A$由此可见,计算矩阵幂的直接幂法需要进行$k-1$次矩阵乘法运算,时间复杂度为$O(kn^3)$。
2. 快速幂法快速幂法是计算矩阵幂的高效方法,它能够有效地减少运算次数,提高计算效率。
该方法基于指数的二进制表示,通过不断地平方和乘以相应的权值,最终计算出矩阵幂的值。
具体步骤如下:(1)将指数$k$转换成二进制数,例如,$k=13$转换成二进制数为$1101$。
幂法和反幂法的matlab实现
幂法和反幂法的matlab实现幂法求矩阵主特征值及对应特征向量摘要矩阵特征值的数值算法,在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题,所以说研究利用数学软件解决求特征值的问题是非常必要的。
实际问题中,有时需要的并不是所有的特征根,而是最大最小的实特征根。
称模最大的特征根为主特征值。
幂法是一种计算矩阵主特征值(矩阵按模最大的特征值)及对应特征向量的迭代方法,它最大的优点是方法简单,特别适用于大型稀疏矩阵,但有时收敛速度很慢。
用java来编写算法。
这个程序主要分成了四个大部分:第一部分为将矩阵转化为线性方程组;第二部分为求特征向量的极大值;第三部分为求幂法函数块;第四部分为页面设计及事件处理。
其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法,再经过一系列的验证和迭代得到结果。
关键字:主特征值;特征向量;线性方程组;幂法函数块POWER METHOD FOR FINDING THE EIGENVALUES AND CORRESPONDING EIGENVECTORS OF THEMATRIXABSTRACTNumerical algorithm for the eigenvalue of matrix, in science and engineering technology, alot of problems in mathematics are attributed matrix characteristic value problem, so that studies using mathematical software to solve the eigenvalue problem is very necessary. In practical problems, sometimes need not all eigenvalues, but the maximum and minimum eigenvalue of real. The characteristic value of the largest eigenvalue of the modulus maximum.Power method is a calculation of main features of the matrix values (matrix according to the characteristics of the largest value) and the corresponding eigenvector of iterative method. It is the biggest advantage is simple method, especially for large sparse matrix, but sometimes the convergence speed is very slow.Using java to write algorithms. This program is divided into three parts: the first part is the matrix is transformed into linear equations; the second part for the sake of feature vector of the maximum; the third part isthe exponentiation function block. The fourth part is the page design and eventprocessing .The basic process is a power law function block by calling the matrix is transformed into linear equations method, after a series of validation and iteration results.Power method for finding the eigenvalues and corresponding eigenvectors of the matrixKey words: Main eigenvalue; characteristic vector; linear equations; power function block、目录1幂法......................................................... . (1)1.1幂法的基本理论和推导 (1)1.2幂法算法的迭代向量规范化 (2)2概要设计........................................................ (3)2.1设计背景 (3)2.2运行流程........................................... . (3)2.3运行环境........................................... (3)3程序详细设计 (4)3.1矩阵转化为线性方程组……..………………………………………. .43.2特征向量的极大值 (4)3.3求幂法函数块............….....…………...…......…………………………3.4界面设计与事件处理..........................................................................4运行过程及结果................................................ (6)4.1 运行过程....................................... ..................………………………………………. .64.2 运行结果................................................ .. (6)4.3 结果分析.......................................... (6)5结论 (7)参考文献 (8)附录 (56)1 幂法设实矩阵nn ijaA ⨯=)(有一个完备的特征向量组,其特征值为n λλλ ,,21,相应的特征向量为nx x x ,,21。
matlab矩阵的n次方的计算程序
文章标题:深度解析:如何编写高效的Matlab矩阵的n次方计算程序在现代科学和工程领域,矩阵的n次方计算是一个频繁出现的问题。
无论是在数学建模、信号处理、图像处理还是优化问题中,都离不开对矩阵的高效运算。
在Matlab中,作为最常用的科学计算软件之一,我们经常需要编写高效的矩阵的n次方计算程序来提高计算效率。
1. 背景介绍矩阵的n次方计算是指将一个矩阵自乘n次,即A^n。
而在Matlab 中,有多种计算矩阵的n次方的方法,包括直接计算、对角化和特征值分解等。
然而,不同的方法在不同的情况下都有其适用性和性能差异。
2. 直接计算法直接计算法是指通过循环将矩阵连乘n次来得到矩阵的n次方。
这种方法简单直接,但在n较大时计算量会很大,效率不高。
3. 对角化法对角化法是将矩阵对角化,然后计算对角矩阵的n次方,最后再将结果反变换回原矩阵。
这种方法适用于对角化后易于计算的情况,但对于一般矩阵来说,可能会增加计算复杂度。
4. 特征值分解特征值分解是将矩阵分解为特征向量和特征值的形式,然后通过特征值的幂计算矩阵的n次方。
这种方法适用于稀疏矩阵和大规模矩阵,但在实际应用中可能会受到精度问题的影响。
5. 个人观点和理解在实际编写Matlab矩阵的n次方计算程序时,我倾向于综合利用直接计算、对角化和特征值分解等方法,根据矩阵的特点和应用的情况选择最合适的计算方式。
并且,我会结合Matlab提供的优化工具和并行计算技术,提高程序的效率和性能。
总结回顾编写高效的Matlab矩阵的n次方计算程序需要综合考虑计算方法的适用性和性能,灵活选择合适的计算方式,并利用Matlab的优化和并行计算技术来提高程序的效率。
通过深入理解矩阵的特性和应用背景,我们可以更好地编写高质量的n次方计算程序,提高程序的性能和可维护性。
以上就是针对Matlab矩阵的n次方计算程序的全面解析,希望对您有所帮助。
在实际科学和工程应用中,矩阵的n次方计算是一项非常重要的任务。
matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量
竭诚为您提供优质文档/双击可除matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量篇一:幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1.幂法简介:当矩阵a满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵a需要满足的条件为:(1)|1||2|...|n|0,i为a的特征值xn(2)存在n个线性无关的特征向量,设为x1,x2,...,1.1计算过程:n对任意向量x,有x(0)(0)iui,i不全为0,则有i1x(k1)ax(k)...ak1x(0)aαiuiαiλik1uik1i1i1nnnk12k1λ1u1()a2u2()anun11k111u1k112|越小时,收敛越快;且当k充分大时,有可见,当|1 (k1)k111u1x(k1)x(k1)(k)x1(k),对应的特征向量即是。
kxx11u12算法实现(1).输入矩阵a,初始向量x,误差限,最大迭代次数n(2).k1,0;y(k)x(k)max(abs(x(k))(3).计算xay,max(x);(4).若||,输出,y,否则,转(5)(5).若kn,置kk1,,转3,否则输出失败信息,停机.3matlab程序代码function[t,y]=lpowera,x0,eps,n)%t为所求特征值,y 是对应特征向量k=1;z=0;%z相当于y=x0./max(abs(x0));%规范化初始向量x=a*y;%迭代格式b=max(x);%b相当于ifabs(z-b) t=max(x);return;endwhileabs(z-b)>epsz=b;y=x./max(abs(x));x=a*y;b=max(x);end[m,index]=max(a(matlab用规范化乘幂法求以下矩阵的按模最大特征值及其特征向量)bs(x));%这两步保证取出来的按模最大特征值t=x(index);%是原值,而非其绝对值。
MATLAB矩阵计算大全
MATLAB与矩阵运算1.矩阵运算(1)矩阵元素的初始化:A=[1 2 3;4,5,6]A=[1 2 34 5 6](2)矩阵运算:A^2,A*A,A/B,A\B,A+B,A-B,a*Aa) 矩阵乘法:A)两个矩阵相乘A*B要求:A的列数和B的行数相等B)矩阵的数乘x*A %x与A的各个元素分别相乘C)点乘 A.*B要求:维数相同的向量或矩阵,对应元素对应相乘D)内积dot(A,B);dot(A,B,dim)% A×B=ATB要求:向量长度或矩阵维数相同(同为m x n维阵)。
b) 矩阵除法:在MATLAB中,有两种矩阵除法运算:\和/,分别表示左除和右除。
如果A矩阵是非奇异方阵,则A\B和B/A运算可以实现。
A\B等效于A矩阵的逆左乘B矩阵,也就是inv(A)*B,相当于A*x = B的解;B/A等效于A矩阵的逆右乘B矩阵,也就是B*inv(A),相当于x*A = B的解。
注意:对于含有标量的运算,两种除法运算的结果相同,如3/4和4\3有相同的值,都等于0.75。
如,设a=[10.5,25],则a/5=5\a=[2.1000 5.0000]。
对于矩阵来说,左除和右除表示两种不同的除数矩阵和被除数矩阵的关系。
对于矩阵运算,一般A\B≠B/A。
c) 矩阵的乘方一个矩阵的乘方运算可以表示成A^x,要求A为方阵,x为标量。
点运算:在MATLAB中,有一种特殊的运算,因为其运算符是在有关算术运算符前面加点,所以叫点运算。
点运算符有.*、./、.\和.^。
两矩阵进行点运算是指它们的对应元素进行相关运算,要求两矩阵的维参数相同。
(3)常见的运算rank(A): 矩阵秩的函数trace(A): 求矩阵的迹的函数det(A):求矩阵的行列式的值inv(A):求矩阵的逆A’:矩阵的转置内置矩阵函数:zeros(3,4);ones(3,4);2.矩阵的特征值与特征向量在MATLAB中,计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有3种:(1) E=eig(A):求矩阵A的全部特征值,构成向量E。
幂矩阵的计算方法
幂矩阵的计算方法
幂矩阵是一种特殊的矩阵乘法,在计算机科学和数学领域被广泛应用。它可以用于解决一些重要的问题,比如图论中的路径计算、网络分析中的节点权重计算等。
幂矩阵的计算方法是通过矩阵的乘法运算来实现的。假设我们有一个n阶矩阵A,我们想要计算A的m次幂矩阵,也就是A的m-1次幂矩阵与A的乘积。
我们需要定义矩阵的乘法运算。矩阵的乘法运算是将两个矩阵的对应们的乘积C可以表示为C = AB,其中C的元素c[i][j]的计算方式为c[i][j] = ∑a[i][k] * b[k][j],其中k的范围是从1到n。
幂矩阵的计算方法在实际应用中有广泛的用途。例如,在图论中,我们可以使用幂矩阵的计算方法来计算图中两个节点之间的路径数量。具体来说,我们可以定义一个邻接矩阵,其中矩阵的元素a[i][j]表示从节点i到节点j的边的数量。然后,我们可以计算邻接矩阵的m次幂矩阵,其中m表示两个节点之间的最短路径的长度。通过这种方式,我们可以快速有效地计算出图中任意两个节点之间的最短路径。
在计算幂矩阵的过程中,我们需要进行多次矩阵乘法运算。假设我们要计算矩阵A的m次幂矩阵,我们可以使用迭代的方法来实现。具体来说,我们首先将A赋值给一个临时矩阵B,然后进行以下操作m-1次:将B与A相乘,将结果赋值给B。最后,B就是A的m次幂矩阵。
需要注意的是,幂矩阵的计算方法要求矩阵A是一个方阵,也就是行数和列数相等。否则,矩阵的乘法运算无法进行。
数值研究分析试验幂法与反幂法matlab
数值分析试验幂法与反幂法matlab————————————————————————————————作者:————————————————————————————————日期:一、问题的描述及算法设计(一)问题的描述我所要做的课题是:对称矩阵的条件数的求解设计1、求矩阵A 的二条件数问题 A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----210121012 2、设计内容:1)采用幂法求出A 的. 2)采用反幂法求出A 的.3)计算A 的条件数 ⅡA Ⅱ2* ⅡA -1Ⅱ2=cond2(A )=/.(精度要求为10-6) 3、设计要求1)求出ⅡA Ⅱ2。
2)并进行一定的理论分析。
(二)算法设计1、幂法算法(1)取初始向量u )0((例如取u )0(=(1,1,…1)T ),置精度要求ε,置k=1.(2)计算v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k(3)若| m k = m 1-k |<ε,则停止计算(m k 作为绝对值最大特征值1λ,u )(k 作为相应的特征向量)否则置k=k+1,转(2)2、反幂法算法(1)取初始向量u )0((例如取u )0(=(1,1,…1)T ),置精度要求ε,置k=1.(2)对A 作LU 分解,即A=LU(3)解线性方程组 Ly )(k =u )1(-k ,Uv )(k =y )(k(4)计算m k =max(v )(k ), u )(k = v )(k / m k(5)若|m k =m 1-k |<ε,则停止计算(1/m k 作为绝对值最小特征值n λ,u )(k 作为相应的特征向量);否则置k=k+1,转(3).二、算法的流程图(一)幂法算法的流程图noyes开始k=0;m1=0 v=A*u [vmax,i]=max(abs(v m=v(i);u=v/mabs(m-m1)< 1e-6 index=1;break; 输出:m ,u ,index 结束m1=m;k=k+1(二)反幂法算法的流程图no yes开始输入A ;[m ,u,index] =pow_inv(A,1e-6) k=0;m1=0 v=invA*u [vmax,i]=max(abs(v)) m=v(i);u=v/mabs(m-m1)< 1e-6 index=1;break; 输出:m ,u ,index结束 m1=m;k=k+1 输入A ;三、算法的理论依据及其推导(一)幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法,也即,绝对值最大的特征值。
MATLAB中的矩阵运算函数
MATLAB中的矩阵运算函数1,round函数函数简介调用格式:Y = round(X)在matlab中round也是一个四舍五入函数。
对数组A中每个元素朝最近的方向取整数部分,并返回与A同维的整数数组B,对于一个复数参量A,则分别对其实部和虚数朝最近的方向取整数部分,并返回一复数数据B。
(1)fix(x) : 截尾取整.>>fix( [3.12 -3.12])ans =3 -3(2)floor(x):不超过x 的最大整数.(高斯取整)>>floor( [3.12 -3.12])ans =3 -4(3)ceil(x) : 大于x 的最小整数>>ceil( [3.12 -3.12])ans =4 -3(4)四舍五入取整>> round(3.12 -3.12)ans =0>> round([3.12 -3.12])ans =3 -32,reshape函数:重新调整矩阵的行数、列数、维数先给上一段代码:>> a=[1 2 3;4 5 6;7 8 9;10 11 12];>> b=reshape(a,2,6);这段代码的结果是这样的:>>a1 2 34 5 67 8 910 11 12>>b1 72 83 94 105 116 12对于 b=reshape(a,m,n);其中的规律是这样的,先把矩阵a按列拆分,然后拼接成一个大小为m*n的向量。
然后对这个向量每隔m间隔取一个元素组成一个向量b_i,之后的向量b_i+1也是这样生成,只不过第一个元素往下移一位。
这样做完之后得到m个大小为n的行向量,将这些行向量拼接即可得到矩阵b。
3,取模(mod)与取余(rem)通常取模运算也叫取余运算,它们返回结果都是余数.rem和mod 唯一的区别在于:当x和y的正负号一样的时候,两个函数结果是等同的;当x和y的符号不同时,rem 函数结果的符号和x的一样,而mod和y一样。
matlab幂法求特征值和特征向量方法实现和函数表示
matlab幂法求特征值和特征向量方法实现和函数表示1. 引言在数值分析中,求解特征值和特征向量是一项重要而且经常出现的任务。
特征值和特征向量在矩阵和线性代数中有着广泛的应用,涉及到许多领域,如机器学习、信号处理、结构动力学等。
在matlab中,幂法是一种常用的求解特征值和特征向量的方法,同时也有对应的函数可以实现这一过程。
2. 幂法的原理幂法是一种迭代方法,它利用矩阵的特征值和特征向量的性质,通过不断地迭代计算,逼近矩阵的主特征值和对应的特征向量。
具体来说,假设A是一个n阶矩阵,它的特征值λ1>λ2≥...≥λn,并且对应着线性无关的特征向量v1,v2,...,vn。
如果选择一个任意的非零初始向量x0,并进行以下迭代计算:```x(k+1) = Ax(k) / ||Ax(k)||```其中,||.||表示向量的模长。
不断迭代计算后,x(k)将收敛到矩阵A的主特征向量v1上,并且相应的特征值即为A的主特征值λ1。
3. matlab实现幂法求解特征值和特征向量在matlab中,幂法的实现也非常简单。
可以使用自带的eig函数,该函数可以直接求解矩阵的特征值和特征向量。
使用方法如下:```[V,D] = eig(A)```其中,A为待求解的矩阵,V为特征向量矩阵,D为特征值矩阵。
利用eig函数,即可一步到位地求解矩阵的特征值和特征向量,非常简单方便。
4. 函数表示幂法求解特征值和特征向量的过程可以表示为一个matlab函数。
通过封装相关的迭代算法和收敛判据,可以方便地实现幂法的函数表示。
可以定义一个名为powerMethod的函数:```matlabfunction [lambda, v] = powerMethod(A, x0, maxIter, tol)% 初始化k = 1;x = x0;% 迭代计算while k <= maxItery = A * x;lambda = norm(y, inf);x = y / lambda;% 检查收敛性if norm(A * x - lambda * x) < tolbreak;endk = k + 1;endv = x;end```利用这个函数,就可以自己实现幂法求解特征值和特征向量的过程。
matlab 矩阵运算程序
matlab 矩阵运算程序MATLAB是一种强大的数学软件,主要用于数值计算、算法开发、数据可视化和数据分析等。
在MATLAB中,矩阵运算是非常常见的操作。
以下是一个简单的MATLAB矩阵运算程序示例:```matlab创建两个矩阵A和BA = [1, 2, 3;4, 5, 6;7, 8, 9];B = [9, 8, 7;6, 5, 4;3, 2, 1];矩阵加法C = A + B;disp('矩阵A和矩阵B的和:');disp(C);矩阵减法D = A - B;disp('矩阵A和矩阵B的差:'); disp(D);矩阵乘法E = A * B;disp('矩阵A和矩阵B的乘积:'); disp(E);矩阵转置T = transpose(A);disp('矩阵A的转置:');disp(T);求矩阵的行列式det_A = det(A);disp('矩阵A的行列式:');disp(det_A);求矩阵的逆矩阵inv_A = inv(A);disp('矩阵A的逆矩阵:');disp(inv_A);求矩阵的秩rank_A = rank(A);disp('矩阵A的秩:');disp(rank_A);求矩阵的特征值eig_A = eig(A);disp('矩阵A的特征值:');disp(eig_A);```以上程序演示了MATLAB中的一些基本矩阵运算,如加法、减法、乘法、转置、求行列式、求逆矩阵、求秩和求特征值等。
您可以根据实际需求修改矩阵A和B的值,然后运行该程序以观察结果。
需要注意的是,这里的矩阵运算都是在MATLAB环境下进行的。
如果要编写比MATLAB更快的矩阵运算程序,可以尝试使用如C、C++等编程语言,并链接到高性能的数学库,如Intel的Math Kernel Library(MKL)。
matlab求矩阵特征值特征向量 乘幂法
摘 要根据现代控制理论课程的特点, 提出并利用MATLAB 设计了现代控制理论课程的实验, 给出了设计的每个实验的主要内容及使用到的MATLAB 函数, 并对其中的一个实验作了详细说明。
通过这些实验, 将有助于学生理解理论知识, 学习利用MATLAB 解决现代控制理论问题。
关键词:现代控制理论、MATLAB 、仿真。
1设计目的、内容及要求1.1设计目的本课程设计以自动控制理论、现代控制理论、MATLAB 及应用等知识为基础,求连续系统对应的离散化的系统,并用计算系数阵按模最大的特征根法判别离散系统的稳定性,目的是使学生在现有的控制理论的基础上,学会用MATLAB 语言编写控制系统设计与分析的程序,通过上机实习加深对课堂所学知识的理解,掌握一种能方便地对系统进行离散化的实现和分析系统的稳定性的设计的工具。
1.2设计内容及要求1 在理论上对连续系统离散化推导出算法和计算公式2 画出计算机实现算法的框图3 编写程序并调试和运行4 以下面的系统为例,进行计算⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----=041020122A ,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=100B ,[]111-=c 5 分析运算结果6 幂法迭代精度为ep=0.001,离散系统展开项数为207 程序应具有一定的通用性,对不同参数能有兼容性。
2算法选择及推导2.1连续系统离散化算法书P67离散化意义已知被控对象的状态方程为:()()()()()()t t u t y t t u t =+=+ xAx B Cx D对方程求解,得:0()()0()()()ott t t t t e t e u d τττ--=+⎰A A x x B设0t kT =,(1)t k T =+,代入上式,得:H 公式若省略T 则为{⎰+-++Φ=+Tk kTd kT Bu T k kt x T T k x )1()(])1[()()(])1([(ττφ不改变与离散后时刻,即得连续离散化方程则:相当于)+=(上限相当于下限设令D C kT Du kT Cx kT y kT t kT u T H kT x T G T k x Bdt t Bdt e T H t T k T t kT d dt T k t Bd e T H e T T G TTAT T k kT T k A AT )()()()()()()(])1([(:)()(0,1,,)1()()()(0)1(])1[(+==+=+Φ=====-=-+=⋅==Φ=⎰⎰⎰+-+ττττττ2.2判别离散系统的稳定性2.2.1方法选择这里选用乘幂法,即求矩阵A 按模最大的特征值和相应的特征向量的方法判别离散系统的稳定性。
如何利用MATLAB进行矩阵运算
如何利用MATLAB进行矩阵运算概述在科学和工程领域,矩阵运算是一项非常重要的技能。
MATLAB作为一种高级数值计算和数据可视化软件,提供了丰富的功能和工具来处理矩阵运算。
本文将介绍如何使用MATLAB进行矩阵运算,包括矩阵的创建、矩阵的运算、矩阵的转置和逆矩阵等。
1. 矩阵的创建在MATLAB中,矩阵可以通过不同的方式进行创建。
最常见的方法是使用"["和"]"符号。
例如,以下命令将创建一个3x3的零矩阵:A = [0 0 0; 0 0 0; 0 0 0]除了手动创建矩阵外,MATLAB还提供了一些内置的函数来创建特殊类型的矩阵。
例如,下面的代码将创建一个单位矩阵:I = eye(3)2. 矩阵的运算使用MATLAB进行矩阵运算非常简单。
可以使用标准的数学运算符来执行加法、减法、乘法和除法等操作。
以下是一些示例代码:A = [1 2 3; 4 5 6; 7 8 9]B = [9 8 7; 6 5 4; 3 2 1]C = A + B % 矩阵加法D = A - B % 矩阵减法E = A * B % 矩阵乘法除了标准的数学运算符,MATLAB还提供了一些特殊的函数来执行矩阵运算。
例如,使用"inv"函数可以计算矩阵的逆矩阵:A = [1 2; 3 4]B = inv(A) % 计算A的逆矩阵3. 矩阵的转置矩阵的转置是指将矩阵的行和列互换。
在MATLAB中,可以使用"'"符号来实现矩阵的转置。
以下是一个示例:A = [1 2 3; 4 5 6; 7 8 9]B = A' % 矩阵A的转置4. 矩阵的逆矩阵逆矩阵是指对于一个方阵A,存在一个方阵B,使得AB=BA=I,其中I是单位矩阵。
在MATLAB中,可以使用"inv"函数来计算矩阵的逆矩阵。
以下是一个示例:A = [1 2; 3 4]B = inv(A) % 计算A的逆矩阵然而需要注意的是,并非所有的矩阵都有逆矩阵。
matlab幂次方
matlab幂次方MATLAB幂次方MATLAB是一种非常强大的编程语言,广泛应用于各种领域,包括科学和工程。
在MATLAB中,幂次方是一种非常常见的数学操作,用于求解各种数学问题。
本文将介绍MATLAB中的幂次方运算,并且提供一些实用的幂次方运算实例。
一、MATLAB中幂次方的基本操作在MATLAB中,幂次方运算用“^”符号表示。
例如,如果想计算2的3次方,可以执行以下操作:>> 2^3ans =8这将返回8,即2的3次方。
同样的,如果想计算一个向量、矩阵或多维数组的幂次方,也可以使用同样的运算符。
例如,如果想计算一个2x2矩阵的2次方,可以执行以下操作:>> A = [1 2; 3 4];>> A^2ans =7 1015 22这将返回一个2x2的矩阵,它是原始矩阵A的2次方。
同样的,如果想计算一个3维数组的幂次方,也可以使用同样的运算符。
二、使用幂次方运算来解决实际问题除了基本的幂次方操作之外,MATLAB中的幂次方运算还可以用于解决许多实际问题。
以下是一些常用的实例:1. 创建一个钟形曲线钟形曲线,也称为高斯曲线,是一种非常常见的函数形状。
可以使用幂次方运算来创建这个曲线。
例如,以下代码将创建一个钟形曲线:>> x = -3:0.1:3;>> y = exp(-x.^2);>> plot(x,y);这将创建一个钟形曲线的图像,其中x的范围是从-3到3,以0.1为步长,y的值是x的平方的相反数的指数函数。
plot(x,y)函数将绘制这个曲线。
2. 使用幂次方运算来求解线性方程组线性方程组是一组n个线性方程,其中每个方程具有n个未知数。
可以使用幂次方运算来求解这个方程组。
例如,以下代码将解决以下线性方程组:3x + y - z = 7x - 2y + 4z = -11x + y - 2z = 1>> A = [3 1 -1; 1 -2 4; 1 1 -2];>> b = [7; -11; 1];>> x = A^(-1)*bx =1.00002.00003.0000这将返回一个列向量x,其中包含了线性方程组的解。
matlab 数组幂运算
matlab 数组幂运算我可以帮助你写一篇关于MATLAB 数组幂运算的1500-2000字的文章。
现在,让我们一起来逐步回答这个问题。
第一步:介绍MATLAB 数组幂运算在使用MATLAB 进行数值计算时,数组幂运算是一项非常重要的操作。
它可以对数组中的每个元素进行幂运算,并且能够快速且高效地执行。
在本文中,我们将学习如何使用MATLAB 对数组进行幂运算,并且将展示一些实际例子来进一步理解数学和编程的奇妙结合。
第二步:数组幂运算的语法和基本概念在MATLAB 中,数组幂运算使用符号“^”进行表示。
该符号放置在两个数组或矩阵之间,并且被视为一种“元素级操作”。
这意味着它将对两个数组(或矩阵)中的相应元素进行幂运算,生成一个新的结果矩阵或数组。
如果输入是两个实数数组,那么幂运算的结果将是每对相应元素的幂运算结果。
例如,如果输入为A 和B 两个数组,则结果数组C 的元素C(i,j) 将是A(i,j) 的B(i,j) 次幂。
对于复数数组,幂运算也是逐元素进行的。
这意味着对于每对相应的复数元素,幂运算将得到一对复数结果。
例如,如果有一个复数数组Z,那么Z(i,j) 的幂运算结果将是Z(i,j) 的第二个输入参数的Z(i,j) 次幂。
第三步:一些示例和实际应用理论还是有些抽象,让我们看一些具体的例子来更好地理解MATLAB 数组幂运算的用法和应用。
例子1:假设有两个输入数组A 和B,分别是2×2 的矩阵。
其中,矩阵A 的元素为[2 4; 1 3],矩阵B 的元素为[3 2; 5 2]。
现在,我们要将这两个矩阵进行幂运算。
在MATLAB 中,我们可以使用如下代码实现:matlabA = [2 4; 1 3];B = [3 2; 5 2];C = A.^B;disp(C);运行上述代码后,输出将是一个2×2 的矩阵,其中的元素为:8 161 9这是因为矩阵A 和B 的对应元素进行了幂运算。
matlab 幂次方
matlab 幂次方
Matlab是一种数学软件,可以用于进行各种数学计算和数据分析。
其中,幂次方是Matlab中的一个基本运算,可以使用“^”符号进行计算。
幂次方运算是指将一个数的某个次方作为结果输出。
例如,2的3次方等于8,即2^3=8。
在Matlab中,可以使用以下方式进行幂次方计算:
1. 直接使用“^”符号进行计算。
例如:2^3=8。
2. 使用power函数进行计算。
power函数的语法格式为:
power(x,y),其中x表示底数,y表示指数。
例如:power(2,3)=8。
除了这两种方法外,在Matlab中还有一些其他的函数可以用于幂次方运算。
例如:
1. exp函数:exp(x)表示e的x次方。
2. log函数:log(x)表示以e为底数的对数。
3. sqrt函数:sqrt(x)表示x的平方根。
需要注意的是,在Matlab中进行幂次方运算时,需要注意数据类型和精度问题。
如果数据类型不正确或者精度不够高,可能会导致结果出现错误。
总之,在Matlab中进行幂次方运算非常简单,只需要使用“^”符号或者相应的函数即可。
同时也需要注意数据类型和精度问题。
幂法及其MATLAB程序
5.2 幂法及其MATLAB 程序5.2.2 幂法的MATLAB 程序用幂法计算矩阵A 的主特征值和对应的特征向量的MATLAB 主程序function [k,lambda,Vk,Wc]=mifa(A,V0,jd,max1)lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0;while ((k<=max1)&(state==1))Vk=A*V; [m j]=max(abs(Vk)); mk=m;tzw=abs(lambda-mk); Vk=(1/mk)*Vk;Txw=norm(V-Vk); Wc=max(Txw,tzw); V=Vk;lambda=mk;state=0;if (Wc>jd)state=1;endk=k+1;Wc=Wc;endif (Wc<=jd)disp('请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('请注意:迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:') endVk=V;k=k-1;Wc;例 5.2.2 用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量,精度510-=ε.并把(1)和(2)输出的结果与例5.1.1中的结果进行比较.(1)⎪⎪⎭⎫ ⎝⎛-=4211A ; (2)⎪⎪⎪⎭⎫ ⎝⎛=633312321B ;(3)⎪⎪⎪⎭⎫ ⎝⎛--=1124111221C ;(4)⎪⎪⎪⎭⎫ ⎝⎛---=20101350144D . 解 (1)输入MATLAB 程序>>A=[1 -1;2 4]; V0=[1,1]';[k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100),[V,D] = eig (A), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =33 3.00000173836804 8.691862856124999e-007Vk = V = wuV =-0.49999942054432 -0.70710678118655 0.44721359549996 -0.894428227562941.00000000000000 0.70710678118655 -0.89442719099992 -0.89442719099992Dzd = wuD =3 1.738368038406435e-006由输出结果可看出,迭代33次,相邻两次迭代的误差W c ≈8.69 19e-007,矩阵A 的主特征值的近似值lambda ≈3.000 00和对应的特征向量的近似向量V k ≈(-0.500 00,1.00000T ), lambda 与例5.1.1中A 的最大特征值32=λ近似相等,绝对误差约为1.738 37e-006,V k 与特征向量X =T22k T )1,21(- )0(2≠k 的第1个分量的绝对误差约等于0,第2个分量的绝对值相同.由wuV 可以看出,2λ的特征向量V (:,2) 与V k 的对应分量的比值近似相等.因此,用程序mifa.m 计算的结果达到预先给定的精度510-=ε.(2) 输入MATLAB 程序>>B=[1 2 3;2 1 3;3 3 6]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(B,V0,0.00001,100), [V,D] = eig (B), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = Dzd = wuD =3 9 0 9 0Vk = wuV =0.50000000000000 0.816496580927730.50000000000000 0.816496580927731.00000000000000 0.81649658092773V =0.70710678118655 0.57735026918963 0.40824829046386-0.70710678118655 0.57735026918963 0.408248290463860 -0.57735026918963 0.81649658092773(3) 输入MATLAB 程序>> C=[1 2 2;1 -1 1;4 -12 1];V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(C,V0,0.00001,100), [V,D] = eig (C), Dzd=max(diag(D)), wuD=abs(Dzd-lambda),Vzd=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示请注意:迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =100 0.09090909090910 2.37758124193119Dzd = wuD =1.00000000000001 0.90909090909091Vk= Vzd = wuV =0.99999999999993 0.90453403373329 0.904534033733350.99999999999995 0.30151134457776 0.301511344577781.00000000000000 -0.30151134457776 -0.30151134457776由输出结果可见,迭代次数k 已经达到最大迭代次数max 1=100,并且lambda 的相邻两次迭代的误差Wc ≈2.377 58>2,由wuV 可以看出,lambda 的特征向量V k 与真值Dzd 的特征向量V zd 对应分量的比值相差较大,所以迭代序列发散.实际上,实数矩阵C 的特征值的近似值为i ,i ,010*********.000321=-==λλλ ,并且对应的特征向量的近似向量分别为X T1=1k (0.90453403373329,0.30151134457776,-0.30151134457776)T ,X =T 22k (-0.72547625011001,-0.21764287503300-0.07254762501100i, 0.58038100008801-0.29019050004400i )T ,X =T33k ( -0.72547625011001, -0.21764287503300 + 0.07254762501100i,0.58038100008801 + 0.29019050004400i)T0,0(21≠≠k k , 03≠k 是常数).(4)输入MATLAB 程序>> D=[-4 14 0;-5 13 0;-1 0 2]; V0=[1,1,1]';[k,lambda,Vk,Wc]=mifa(D,V0,0.00001,100), [V,Dt] =eig (D), Dtzd=max(diag(Dt)), wuDt=abs(Dtzd-lambda),Vzd=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下:k = lambda = Wc =19 6.00000653949528 6.539523793591684e-006Dtzd = wuDt =6.00000000000000 6.539495284840768e-006Vk = Vzd = wuV =0.79740048053564 0.79740048053564 0.797400480535640.71428594783886 0.56957177181117 0.79740021980618-0.24999918247180 -0.19935012013391 0.797403088133705.3 反幂法和位移反幂法及其MATLAB程序5.3.3 原点位移反幂法的MATLAB程序(一)原点位移反幂法的MATLAB主程序1用原点位移反幂法计算矩阵A的特征值和对应的特征向量的MATLAB主程序1 function [k,lambdan,Vk,Wc]=ydwyfmf(A,V0,jlamb,jd,max1)[n,n]=size(A); A1=A-jlamb*eye(n); jd= jd*0.1;RA1=det(A1);if RA1==0disp('请注意:因为A-aE的n阶行列式hl等于零,所以A-aE不能进行LU分解.')returnendlambda=0;if RA1~=0for p=1:nh(p)=det(A1(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意:因为A-aE的r阶主子式等于零,所以A-aE不能进行LU分解.')returnendendif h(1,i)~=0disp('请注意:因为A-aE的各阶主子式都不等于零,所以A-aE 能进行LU分解.')k=1;Wc =1;state=1; Vk=V0;while((k<=max1)&(state==1))[L U]=lu(A1); Yk=L\Vk;Vk=U\Yk; [mj]=max(abs(Vk));mk=m;Vk1=Vk/mk; Yk1=L\Vk1;Vk1=U\Yk1;[m j]=max(abs(Vk1));mk1=m;Vk2=(1/mk1)*Vk1;tzw1=abs((mk-mk1)/mk1);tzw2=abs(mk1-mk);Txw1=norm(Vk)-norm(Vk1);Txw2=(norm(Vk)-norm(Vk1))/norm(Vk1);Txw=min(Txw1,Txw2); tzw=min(tzw1,tzw2);Vk=Vk2;mk=mk1; Wc=max(Txw,tzw);Vk=Vk2;mk=mk1;state=0;if(Wc>jd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif (Wc<=jd)disp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k 已经达到最大迭代次数max1,按模最小特征值的迭代值lambda,特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下:')endhl,RA1endend[V,D]=eig(A,'nobalance'),Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;例5.3.2 用原点位移反幂法的迭代公式(5.28),根据给定的下列矩阵的特征值n λ的初始值n λ~,计算与n λ对应的特征向量n X 的近似向量,精确到0.000 1. (1)⎪⎪⎪⎭⎫ ⎝⎛----210242011,2.0~2=λ;(2)⎪⎪⎭⎫ ⎝⎛-4211,001.2~2=λ;(3)⎪⎪⎪⎭⎫ ⎝⎛--3315358215211,8.26~3=λ.解 (1)输入MATLAB 程序>> A=[1 -1 0;-2 4 -2;0 -1 2];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.2,0.0001,10000)运行后屏幕显示结果 请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = hl =3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D =1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 00.7616 1.0000 -0.7616 0.3633 0 0.2384 00.4323 -0.3200 -0.4323 1.0000 0 0 1.6367(2)输入MATLAB 程序>> A=[1 -1;2 4];V0=[20,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001,0.0001,100)运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc = hl =2 2.0020 5.1528e-007 -1.0010 -0.0010Vk = V = D =1.0000 -1.0000 0.5000 2 0-1.0000 1.0000 -1.0000 0 3(3)输入MATLAB 程序>> A=[-11 2 15;2 58 3;15 3 -3];V0=[1,1,-1]';[k,lambdan,Vk,Wc]=ydwyfmf(A,V0,8.26, 0.0001,100)运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambdan= Wc = hl =2 8.2640 6.9304e-008 -19.2600 -961.9924 -6.1256Vk = V = D =-0.7692 0.7928 0.6081 0.0416 -22.5249 0 00.0912 0.0030 -0.0721 0.9974 0 8.2640 0-1.0000 -0.6095 0.7906 0.0590 0 0 58.2609例 5.3.3 用原点位移反幂法的迭代公式(5.28),计算⎪⎪⎪⎭⎫ ⎝⎛-----=1026471725110A 的分别对应于特征值 1.001~11=≈λλ,.001 2~22=≈λλ, 001.4~33=≈λλ的特征向量1X ,2X ,3X 的近似向量,相邻迭代误差为0.001.将计算结果与精确特征向量比较. 解 (1)计算特征值 1.001~11=≈λλ对应的特征向量1X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]= ydwyfmf(A,V0,1.001, 0.001,100),[V,D]=eig(A);Dzd=min(diag(D)), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-1.00100000000000 5.98500100000000 -0.00299600100000k = lambda = RA1 =5 1.00200000000000 -0.00299600100000Vk = VD = wuV =-0.50000000000000 -0.40824829046386 0.81649658092773-0.50000000000000 -0.40824829046386 0.81649658092773-1.00000000000000 -0.81649658092773 0.81649658092773Wc = Dzd = wuD =1.378794763695562e-009 1.00000000000000 0.00200000000000 从输出的结果可见,迭代5次,特征向量1X 的近似向量1~X 的相邻两次迭代的误差Wc ≈1.379 e-009,由wuV 可以看出,1~X = Vk 与VD 的对应分量的比值相等.特征值1λ的近似值lambda ≈1.002与初始值=1~λ 1.001的绝对误差为0.001,而与 1λ的绝对误差为0.002,其中 =1X T )000000000001.000 , 000000000000.500- , 000000000000.500( -, =1~X T )000000000001.000 , 000000000000.500- , 000000000000.500(-. (2)计算特征值.001 2~22=≈λλ对应特征向量2X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,2.001, 0.001,100) ,[V,D]=eig(A); WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-2.00100000000000 -8.01299900000000 0.00200099900000k = Wc = lambda = WD =2 3.131363162302120e-007 2.00200000000016 0.00200000000016Vk = VD = wuV =-0.24999999999999 0.21821789023599 -0.87287156094401 -0.49999999999999 0.43643578047198 -0.87287156094398 -1.00000000000000 0.87287156094397 -0.87287156094397 从输出的结果可见,迭代2次,特征向量2X 的近似向量2~X 的相邻两次迭代的误差Wc ≈3.131e-007,2~X 与2X 的对应分量的比值近似相等.特征值2λ的近似值lambda ≈2.002与初始值=2~λ 2.001的绝对误差约为0.001,而lambda 与2λ的绝对误差约为0.002,其中 =2~X T )00000000000000.1,99999999999499.0,99999999999249.0(---, =2X T ) 000000000001.000- ,000000000000.500- ,99999999999-0.249( . (3)计算特征值 001.4~33=≈λλ对应特征向量3X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,4.001, 0.001,100)[V,D]=eig(A);WD=lambda-max(diag(D)),VD=V(:,3),wuV=V(:,3)./Vk,运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行L U 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-4.00100000000000 -30.00899900000000 -0.00600500099999 k = lambda = Wc = WD =2 4.00199999999990 1.996084182914842e-007 0.00199999999990Vk = VD = wuV =0.40000000000001 -0.32444284226153 -0.81110710565380 0.60000000000001 -0.48666426339229 -0.81110710565381 1.00000000000000 -0.81110710565381 -0.81110710565381 从输出的结果可见,迭代2次,特征向量3X 的近似向量3~X 的相邻两次迭代的误差Wc ≈1.996e-007,3~X 与3X 的对应分量的比值近似相等.特征值3λ的近似值 4.001~4.0022=≈λ与初始值lambda 的绝对误差近似为001.0,而lambda 与3λ的绝对误差约为0.002,其中 =3X (-0.400 000 000 000 00,-0.600 000 000 000 00,-1.000 000 000 000 00T ), =3~X T )000000000001.000 ,100000000000.600 ,10000000000.400(.(二)原点位移反幂法的MATLAB 主程序2用原点位移反幂法计算矩阵A 的特征值和对应的特征向量的MATLAB 主程序2function [k,lambdan,Vk,Wc]=wfmifa1(A,V0,jlamb,jd,max1)[n,n]=size(A); jd= jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1); lambda1=0;k=1;Wc =1;state=1; U=V0;while ((k<=max1)&(state==1))Vk=A1\U; [m j]=max(abs(Vk)); mk=m; Vk=(1/mk)*Vk;Vk1=A1\Vk;[m1 j]=max(abs(Vk1)); mk1=m1,Vk1=(1/mk1)*Vk1;U=Vk1,Txw=(norm(Vk1)-norm(Vk))/norm(Vk1);tzw=abs((lambda1-mk1)/mk1);Wc=max(Txw,tzw); lambda1=mk1;state=0;if (Wc>jd)state=1;endk=k+1;endif (Wc<=jd)disp('请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:')elsedisp('请注意迭代次数k 已经达到最大迭代次数max1, 特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:') end[V,D] =eig(A,'nobalance'),Vk=U;k=k-1;Wc;lambdan=jlamb+1/mk;例5.3.4 用原点位移反幂法的迭代公式(5.27),计算例题5.3.3,并且将这两个例题的计算结果进行比较.再用两种原点位移反幂法的MATLAB 主程序,求979999999990.999~1=λ对应的特征向量. 解 (1)计算特征值 1.001~11=≈λλ对应特征向量1X 的近似向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0,1.001,0.001,100)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = lambda = Wc =5 1.00200000000138 1.376344154436924e-006Vk’ = -0.50000000000000 -0.50000000000000 -1.00000000000000同理可得,另外与两个特征值对应的特征向量.(2)再用两种原点位移反幂法的MATLAB 主程序,求979999999990.999~1=λ对应的特征向量.输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=ydwyfmf(A,V0,0.99999999999997,0.001,100) 运行后屏幕显示结果请注意:因为A-aE 的各阶主子式都不等于零,所以A-aE 能进行LU 分解.A-aE 的秩R(A-aE)和各阶顺序主子式值hl 、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:hl =-0.99999999999997 6.00000000000045 0.00000000000010RA1 = 1.039168751049192e-013 k = 2 lambda = 1.00000000000000输入MATLAB 程序>> A=[0 11 -5;-2 17 -7;-4 26 -10];V0=[1,1,1]';[k,lambda,Vk,Wc]=wfmifa1(A,V0, 0.99999999999997,0.001,100) 运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下:k = 3 lambda = 1.00000000000000 Wc =5.412337245047640e-016Vk = 0.50000000000000 0.50000000000000 1.00000000000000 Wc = 4.317692037236759e-013 Vk =0.500000000000000.500000000000001.000000000000005.4 雅可比(Jacobi)方法及其MATLAB 程序5.4.3 雅可比方法的MATLAB 程序用雅可比方法计算对称矩阵A 的特征值和对应的特征向量的MATLAB 主程序function [k,Bk,V,D,Wc]=jacobite(A,jd,max1)[n,n]=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n); Aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(Aij);[m2 j]=max(m1);i=i(j);while ((k<=max1)&(state==1))k=k+1,aij=abs(Bk-diag(diag(Bk)));[m1 i]=max(abs(aij));[m2 j]=max(m1);i=i(j),j,Aij=(Bk-diag(diag(Bk)));mk=m2*sign(Aij(i,j)),Wc=m2,Dk=diag(diag(Bk));Pk=P0;c=(Bk(j,j)-Bk(i,i))/(2*Bk(i,j)),t=sign(c)/(abs(c)+sqrt(1+c^2)),pii=1/( sqrt(1+t^2)), pij=t/( sqrt(1+t^2)),Pk(i,i)=pii;Pk(i,j)=pij;Pk(j,j)=pii; Pk(j,i)=-pij;Pk,B1=Pk'*Bk;B2=B1*Pk; Vk=Vk*Pk,Bk=B2,if (Wc>jd)state=1;elsereturnendPk;Vk;Bk=B2;Wc;endif (k>max1)disp('请注意迭代次数k 已经达到最大迭代次数max1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:')elsedisp('请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:')endWc;k=k; V=Vk;Bk=B2;D=diag(diag(Bk));[V1,D1]=eig(A,'nobalance')例5.4.2 用雅可比方法的MATLAB 程序计算矩阵⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=12101152302756135612A 的特征值i λ和对应的特征向量i X (4,3,2,1=i ).解 (1)保存名为jacobite.m 为M 文件;(2)输入MATLAB 程序>> A=[12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12];[k,B,V,D,Wc]=jacobite(A,0.001,100)(3)运行后屏幕显示如下:k = i = j = mk = Wc =1 2 1 -56 56c = t =-0.04464285714286 -0.95635313919972pii = pij =0.72270271801843 -0.69115901308510Pk =0.72270271801843 0.69115901308510 0 0 -0.69115901308510 0.72270271801843 0 0 0 0 1.00000000000000 00 0 0 1.00000000000000Vk =0.72270271801843 0.69115901308510 0 0 -0.69115901308510 0.72270271801843 0 00 0 1.00000000000000 00 0 0 1.00000000000000Bk =65.55577579518456 0 0.78579012788509 -0.72270271801843 -0.00000000000001 -46.55577579518456 3.51888247529217 -0.691159013085100.78579012788509 3.51888247529217 5.00000000000000 1.00000000000000 -0.72270271801843 -0.69115901308510 1.00000000000000 12.00000000000000 k =i = j = mk = Wc =2 3 2 3.51888247529217 3.51888247529217c = t =-7.32558932518824 -0.06793885568129pii = pij =0.99770011455446 -0.06778260409592Pk =1.00000000000000 0 0 00 0.99770011455446 0.06778260409592 00 -0.06778260409592 0.99770011455446 00 0 0 1.00000000000000Vk =0.72270271801843 0.68956942653035 0.04684855775127 0 -0.69115901308510 0.72104058455581 0.04898667221449 00 -0.06778260409592 0.99770011455446 00 0 0 1.00000000000000Bk =65.55577579518456 -0.05326290114092 0.78398290060672 -0.72270271801843 -0.05326290114093 -46.79484464383285 0 -0.757352030626270.78398290060672 0.00000000000000 5.23906884864829 0.95085155680318 -0.72270271801843 -0.75735203062627 0.95085155680318 12.00000000000000 k = i = j = mk = Wc =3 4 3 0.95085155680318 0.95085155680318c = t =-3.55519802380213 -0.13796227443116pii = pij =0.99061693994324 -0.13666776612460Pk =1.00000000000000 0 0 00 1.00000000000000 0 00 0 0.99061693994324 0.136667766124600 0 -0.13666776612460 0.99061693994324 Vk =0.72270271801843 0.68956942653035 0.04640897492032 0.00640268773403 -0.69115901308510 0.72104058455581 0.04852702732712 0.006694899061430 -0.06778260409592 0.98833863446096 0.136353445918420 0 -0.13666776612460 0.99061693994324 Bk =65.55577579518456 -0.05326290114092 0.87539690801061 -0.60877636330628 -0.05326290114093 -46.79484464383285 0.10350561019562 -0.750245751038800.87539690801061 0.10350561019562 5.10788720522532 -0.00000000000000 -0.60877636330628 -0.75024575103880 -0.00000000000000 12.13118164342297 k =i = j = mk = Wc =4 1 3 0.87539690801061 0.87539690801061c = t =-34.52598931799430 -0.01447880833914pii = pij =0.99989519853186 -0.01447729093877Pk =0.99989519853186 0 -0.01447729093877 00 1.00000000000000 0 00.01447729093877 0 0.99989519853186 00 0 0 1.00000000000000Vk =0.72329885394465 0.68956942653035 0.03594133368062 0.00640268773403 -0.69038403871280 0.72104058455581 0.05852805174080 0.006694899061430.01430846595712 -0.06778260409592 0.98823505512105 0.13635344591842-0.00197857901214 0 -0.13665344314206 0.99061693994324Bk =65.56845049923633 -0.05175883827808 -0.00000000000000 -0.60871256264964-0.05175883827809 -46.79484464383285 0.10426586517177 -0.75024575103880-0.00000000000000 0.10426586517177 5.09521250117356 0.00881343252823-0.60871256264964 -0.75024575103880 0.00881343252823 12.13118164342297 k = i = j = mk = Wc =5 4 2 -0.75024575103880 0.75024575103880c = t =39.27114962375084 0.01272992971264pii = pij =0.99991898429114 0.01272889838836Pk =1.00000000000000 0 0 00 0.99991898429114 0 -0.012728898388360 0 1.00000000000000 00 0.01272889838836 0 0.99991898429114Vk =0.72329885394465 0.68959505973603 0.03594133368062 -0.00237529014628-0.69038403871280 0.72106738763160 0.05852805174080 -0.002483695665250.01430846595712 -0.06604148348220 0.98823505512105 0.13720519702737-0.00197857901214 0.01260946237032 -0.13665344314206 0.99053668440964Bk =65.56845049923633 -0.05950288535679 -0.00000000000000 -0.60800441437674-0.05950288535680 -46.80439521951078 0.10436960328590 0.00000000000000-0.00000000000000 0.10436960328590 5.09521250117356 0.00748552889860-0.60800441437674 0.00000000000000 0.00748552889860 12.14073221910090 k =i = j = mk = Wc =6 4 1 -0.60800441437674 0.60800441437674c = t =-43.93694931878409 -0.01137847012503pii = pij =0.99993527149402 -0.01137773361366Pk =0.99993527149402 0 0 0.011377733613660 1.00000000000000 0 00 0 1.00000000000000 0-0.01137773361366 0 0 0.99993527149402Vk =0.72327906130899 0.68959505973603 0.03594133368062 0.00585436528595-0.69031109235777 0.72106738763160 0.05852805174080 -0.010338540582940.01274645560931 -0.06604148348220 0.98823505512105 0.13735911385404-0.01324851347145 0.01260946237032 -0.13665344314206 0.99045005670500Bk =65.57536865930122 -0.05949903382392 -0.00008516835377 -0.00000000000000-0.05949903382393 -46.80439521951078 0.10436960328590 -0.00067700797883-0.00008516835377 0.10436960328590 5.09521250117356 0.00748504437150-0.00000000000000 -0.00067700797883 0.00748504437150 12.13381405903603 k =i = j = mk = Wc =7 3 2 0.10436960328590 0.10436960328590c = t =-2.486337309269764e+002 -0.00201098208240pii = pij =0.99999797798167 -0.00201097801616Pk =1.00000000000000 0 0 00 0.99999797798167 0.00201097801616 00 -0.00201097801616 0.99999797798167 00 0 0 1.00000000000000…………………………………………………………………………请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D 如下:V1 =0.68990429476497 -0.03732423222484 0.00588594854431 -0.722913771734500.72058252860300 -0.05998661236737 -0.01028322161977 0.69069289931337-0.06802029759277 -0.98795368410472 0.13841044442471 -0.012779125692250.01288885768193 0.13768088498200 0.99030407443219 0.01325486405899D1 =-46.80463661419736 0 0 00 5.09541442877727 0 00 0 12.13382202426702 00 0 0 65.57540016115307k =10B =65.57540016045945 0.00000000000175 -0.00020481967566 0.000000148628360.00000000000175 -46.80463661419739 0.00000062739984 0.00000000000000-0.00020481967566 0.00000062739984 5.09541442947090 -0.000000000007370.00000014862836 -0.00000000000000 -0.00000000000737 12.13382202426704V =0.72291389811507 0.68990429521617 0.03732177568689 0.00588595055487-0.69069269613201 0.72058252932816 0.05998894273570 -0.010283223540620.01278247108107 -0.06802028564977 0.98795364164379 0.13841044446122-0.01325533307898 0.01288885601755 -0.13768084024946 0.99030407439520D =65.57540016045945 0 0 00 -46.80463661419739 0 00 0 5.09541442947090 00 0 0 12.13382202426704Wc =6.920584967017158e-0045.5 豪斯霍尔德(Householder)方法及其MATLAB程序5.5.1 豪斯霍尔德方法及其MATLAB程序求初等反射矩阵P,使得PX的第一个分量以外的其余的分量都为零的MATLAB主程序function [xigema,rou,miou,P,PX]=Householder(X)n=size(X);nX=norm(X,2);xigema=nX*sign(X(1));rou=xigema*(xigema+X(1));miou=[xigema,zeros(1,n-1)]'+X,E=eye(n,n); C=2*miou*(miou)';P=E-C/(norm(miou,2)^2); PX=P*X;例5.5.1设向量=X()T1,2,2,确定一个初等反射矩阵P,使得PX的后两个分量为零.解输入MATLAB程序>> X=[2 2 1]'; [xigema,rou,miou,P,PX]=Householder(X)运行后屏幕显示结果P = PX =-0.6667 -0.6667 -0.3333 -3.0000-0.6667 0.7333 -0.1333 0.0000-0.3333 -0.1333 0.9333 0.00005.5.2 矩阵约化为上豪斯霍尔德矩阵及其MATLAB程序用豪斯霍尔德变换将n阶矩阵A规约成上豪斯霍尔德矩阵的MATLAB主程序function [k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)n=size(A); Ak=A;for k=1:n-2k,Sk=norm(Ak(k+1:n,k))*sign(Ak(k+1,k)),uk= Ak(k+1:n,k)+ Sk*eye(n-k,1),ck=(norm(uk,2)^2)/2,Pk= eye(n-k,n-k)-uk*uk'/ck,Uk=[eye(k,k),zeros(k,n-k);zeros(n-k, k),Pk],A1=Uk*Ak;Ak=A1,end例5.5.3 用初等反射矩阵正交相似约化实矩阵A 为上豪斯霍尔德矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=34 19- 37 78- 41- 31 11 72- 98 10.2- 78- 32-94- 21 12 1 0 1- 63- 72 1 5 2 3 17- 32 02 7 56- 51- 17 12- 34 52- 12A . 解 输入MATLAB 程序>> A=[12 -52 34 -12 17 -51;-56 7 2 0 32 -17;3 2 5 1 72 -63;-1 0 1 12 21 -94;-32 -78 -10.2 98 -72 11;31 -41 -78 37 -19 34];[k,Sk,uk,ck,Pk,Uk,Ak]=Householdrer1(A)运行后屏幕显示结果k = Sk = ck =1 -71.6310 9.1423e+003uk = Pk =-127.6310 -0.7818 0.0419 -0.0140 -0.4467 0.43283.0000 0.0419 0.9990 0.0003 0.0105 -0.0102-1.0000 -0.0140 0.0003 0.9999 -0.0035 0.0034-32.0000 -0.4467 0.0105 -0.0035 0.8880 0.108531.0000 0.4328 -0.0102 0.0034 0.1085 0.8949Uk =1.0000 0 0 0 0 00 -0.7818 0.0419 -0.0140 -0.4467 0.43280 0.0419 0.9990 0.0003 0.0105 -0.01020 -0.0140 0.0003 0.9999 -0.0035 0.00340 -0.4467 0.0105 -0.0035 0.8880 0.10850 0.4328 -0.0102 0.0034 0.1085 0.8949Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.76430.0000 1.8892 5.7655 1.6556 72.7134 -63.9112-0.0000 0.0369 0.7448 11.7815 20.7622 -93.6963-0.0000 -76.8184 -18.3655 91.0066 -79.6101 20.71910.0000 -42.1447 -70.0897 43.7749 -11.6277 24.5846k = Sk = ck =2 87.6402 7.8464e+003uk = Pk =89.5295 -0.0216 -0.0004 0.8765 0.48090.0369 -0.0004 1.0000 0.0004 0.0002-76.8184 0.8765 0.0004 0.2479 -0.4126-42.1447 0.4809 0.0002 -0.4126 0.7736Uk =1.0000 0 0 0 0 00 1.0000 0 0 0 00 0 -0.0216 -0.0004 0.8765 0.48090 0 -0.0004 1.0000 0.0004 0.00020 0 0.8765 0.0004 0.2479 -0.41260 0 0.4809 0.0002 -0.4126 0.7736Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.7643-0.0000 -87.6402 -49.9272 100.7790 -76.9476 31.4002-0.0000 -0.0000 0.7219 11.8223 20.7005 -93.6570-0.0000 0.0000 29.4202 5.9564 48.8026 -61.06030.0000 0.0000 -43.8731 -2.8860 58.8230 -20.2818…………………………………………………………………………k = Sk = ck =4 -12.2088 195.0398uk = Pk =-15.9753 -0.3085 0.951211.6133 0.9512 0.3085Uk =1.0000 0 0 0 0 00 1.0000 0 0 0 00 0 1.0000 0 0 00 0 0 1.0000 0 00 0 0 0 -0.3085 0.95120 0 0 0 0.9512 0.3085Ak =12.0000 -52.0000 34.0000 -12.0000 17.0000 -51.000071.6310 11.7128 -30.5678 -27.8930 1.6473 21.7643-0.0000 -87.6402 -49.9272 100.7790 -76.9476 31.40020.0000 -0.0000 -52.8292 -5.8754 21.3902 18.44030.0000 0.0000 0.0000 12.2088 40.2435 -106.81340.0000 0.0000 -0.0000 0.0000 64.7555 -34.09095.5.3 实对称矩阵的三对角化及其MATLAB程序将n阶实对称矩阵A规约成三对角形式的MATLAB主程序function T=house(A)[n,n]=size(A);for k=1:n-2s=norm(A(k+1:n,k),2);if (A(k+1,k)<0)s=-s;endr=sqrt(2*s*(A(k+1,k)+s));U(1:k)=zeros(1,k);U(k+1)=(A(k+1,k)+s)/r;U(k+2:n)=A(k+2:n,k)'/r;V(1:k)=zeros(1,k);V(k+1:n)=A(k+1:n,k+1:n)*U(k+1:n)';C=U(k+1:n)*V(k+1:n)';P(1:k)=zeros(1,k);P(k+1:n)=V(k+1:n)-C*U(k+1:n);A(k+2:n,k)=zeros(n-k-1,1);A(k,k+2:n)=zeros(1,n-k-1);A(k+1,k)=-s; A(k,k+1)=-s;A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-2*U(k+1:n)'*P(k+1:n)-2*P( k+1:n)'*U(k+1:n);endT=A;例5.5.4 用初等反射矩阵正交相似约化实对称矩阵A为三对角矩阵.其中⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------------------=5261215121416134237299021237312611451721253233219612371564901435612A 解 输入MATLAB 程序>> A=[12 -56 3 -14 -90 -4;-56 71 23 61 -9 -21;3 23 53 12 -72 51;-14 61 12 73 23 21;-90 -9 -72 23 -34 -61;-41 -21 51 21 -61 -52];T=house(A)运行后屏幕显示结果T =12.0000 114.5513 0 0 0 0114.5513 -43.2395 -108.2763 0 0 00 -108.2763 49.7411 -22.7766 0 00 0 -22.7766 40.2476 -89.1355 00 0 0 -89.1355 44.9606 39.30900 0 0 0 39.3090 19.29025.6 QR 方法及其MATLAB 程序5.6.5 最末元位移QR 法计算实对称矩阵特征值及其MATLAB 程序用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序function tzg=qr4(A,t,max1)[n,n]=size(A); k=0;Ak=A;tzg=zeros(n); state=1;for i=1:n;while ((k<=max1)&(state==1)&(n>1))b1=abs(Ak(n,n-1)); b2=abs(Ak(n,n));b3=abs(Ak(n-1,n-1));b4=min(b2, b3); jd=10^(-t); jd1=jd*b4;if (b1>=jd1)sk=Ak(n,n); Bk=Ak-sk*eye(n); [Qk,Rk]=qr(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp('请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,')disp(' Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解,At=Rk*Qk+sk*eye(n)迭代矩阵:')i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;elsedisp('请注意:i 表示求第i 个特征值,tzgk 是矩阵A 的特征值的近似值,k 是迭代次数,')disp(' 下面的矩阵B 是m 阶矩阵At 的(m-1)阶主子矩阵,即At 降一阶.')i,tzgk=Ak(n,n),tzg(n,1)=tzgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state==1; i=i+1;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1));tzgk=Akdisp('请注意:n 阶实对称矩阵A 的全部真特征值lamoda 和至少含t个有效数字的近似特征值tzg 如下:')lamoda=sort(eig(A))例5.6.5 用最末元位移QR 方法求下列实对称矩阵的全部近似特征值,并将计算结果与A 全部真特征值比较.其中,2 1 1 1 1 3 1 21 1 4- 21 2 2 5)1(⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=A 精度为=ε510-; ,52612151214161342372990212373126114517212532332196123715641901435612)2(⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛------------------=A 精度为=ε410-.解 (1)首先保存用最末元位移QR 方法求实对称矩阵A 全部特征值的MATLAB 主程序为M 文件,取名为qr4.m.在MATLAB 工作窗口输入程序>> A=[5 2 2 1;2 -4 1 1;2 1 3 1;1 1 1 2]; tzg=qr4(A,5,100) 运行后屏幕显示结果请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,Bk=Ak-sk*eye(n),Qk 和Rk 是Bk 的QR 分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =1k =1sk =2Bk =3 2 2 12 -6 1 12 1 1 11 1 1 0Qk =-0.70710678118655 0.38807526285317 0.12674485010490 -0.57735026918963-0.47140452079103 -0.87963726246718 0.06337242505245 0-0.47140452079103 0.20697347352169 -0.63372425052448 0.57735026918963-0.23570226039552 0.18110178933148 0.76046910062937 0.57735026918963 Rk =-4.24264068711929 0.70710678118655 -2.59272486435067 -1.649915822768610 6.44204936336256 0.28458852609232 -0.284588526092320 0 0.44360697536713 -0.443606975367130 0 0 0.00000000000000At =6.27777777777778 -3.10388935193069 -0.10455916682125 0.00000000000000-3.10388935193069 -3.65930388219545 0.01147685957127 0.00000000000000-0.10455916682125 0.01147685957127 1.38152610441767 0.00000000000000 -0.00000000000000 0.00000000000000 0.00000000000000 2.00000000000000 请注意:i 表示求第i 个特征值,tzgk 是矩阵A 的特征值的近似值,k 是迭代次数,下面的矩阵B 是m 阶矩阵At 的(m-1)阶主子矩阵,即At 降一阶.i =1tzgk =2.00000000000000k =1sk =2B =6.27777777777778 -3.10388935193069 -0.10455916682125-3.10388935193069 -3.65930388219545 0.01147685957127-0.10455916682125 0.01147685957127 1.38152610441767请注意:下面的i 表示求第i 个特征值,k 是迭代次数,sk 是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =2k =2sk =1.38152610441767Bk =4.89625167336011 -3.10388935193069 -0.10455916682125-3.10388935193069 -5.04082998661312 0.01147685957127-0.10455916682125 0.01147685957127 0Qk =-0.84445320114929 -0.53537837009187 0.016394874396770.53532568873289 -0.84460953959679 -0.007818734217300.01803324849744 0.00217404228940 0.99983502413586Rk =-5.79813264571247 -0.07718952005739 0.094439180886190 5.91931326753920 0.046285251232420 0 -0.00180396892170At =6.23815929000691 3.16959512520840 -0.000032531419853.16959512520840 -3.61788172311421 -0.00000392190472-0.00003253141985 -0.00000392190472 1.37972243310730请注意:i表示求第i个特征值,tzgk是矩阵A的特征值的近似值,k是迭代次数,下面的矩阵B是m阶矩阵At的(m-1)阶主子矩阵,即At降一阶.i =2tzgk =1.37972243310730k =2sk =1.38152610441767B =6.23815929000691 3.169595125208403.16959512520840 -3.61788172311421请注意:下面的i表示求第i个特征值,k是迭代次数,sk是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =3k =3sk =-3.61788172311421Bk =9.85604101312112 3.169595125208403.16959512520840 0Qk =-0.95198403663348 -0.30614766697629-0.30614766697629 0.95198403663348Rk =-10.35315786173815 -3.017403961789690 -0.97036415284199At =7.16193047323385 0.297074721510000.29707472151000 -4.54165290634115请注意:下面的i表示求第i个特征值,k是迭代次数,sk是原点位移量,Bk=Ak-sk*eye(n),Qk和Rk是Bk的QR分解,At=Rk*Qk+sk*eye(n)迭代矩阵:i =3k =4sk =-4.54165290634115Bk =11.70358337957500 0.297074721510000.29707472151000 0。
matlab反幂法
matlab反幂法
Matlab反幂法(inverse power method)是一种用于求解矩阵的特征值和特征向量的迭代算法。
它可以用于求解实对称矩阵或者复共轭对称矩阵的最小特征值和对应的特征向量。
算法步骤如下:
1.选取一个与所需特征值最接近的初始特征向量x0,通常可以选取一个随机的向量。
2.通过将矩阵A减去所需特征值的一个估计量来构造一个新的矩阵B,即B=A-μI,其中μ是我们的估计量,I为单位矩阵。
3.通过求解线性方程组Bx1=y来计算新的特征向量x1,其中y是单位向量。
4.将x1进行归一化处理,得到新的特征向量x1/||x1||。
5.计算特征值的新估计量μ1=x1^T Ax1/ x1^T x1。
如果μ1与原估计量μ的差距小于要求精度,则停止迭代,否则返回第2步。
反幂法的核心思想是将原问题转化为求解矩阵的倒数,从而求得矩阵的最小特征
值和对应的特征向量。
该方法在求解大型稀疏矩阵的特征值问题方面具有很好的效果。