孤立特征值情况的矩阵摄动法算例(matlab编程)
设原始特征值问题是
)~1(0
0000n i x M x K i i i ==λ
(1)
相应的正交规范条件是
)~1,(000n j i x M x ij
j T i ==δ
(2)
先考虑一个五个自由度的系统．
其中，0M 、0K 分别是原系统的n n ?阶对称质量阵、刚度阵，0i λ是特征值，且200i i ωλ=，
0i ω是固有频率，0i x 是相应的特征向量。
0111
10.5M ??
?
?
?= ?
?
???
021121121
12111K -?? ?-- ?
?=-- ?
-- ?
?-??
系统结构参数改变后，相应的质量阵、刚度阵均有相应的变化，设系统修改后的质量阵、刚
度阵分别为
??
?+=+=1
01
0K K K M M M εε (3)
其中取
1111
10.5M ?? ?
?
?= ?
?
???
1100
00K -?? ?
?
?= ?
?
???
显然，新系统（即结构修改后的系统或称摄动系统）的特征值问题及相应的正交规范条件是 i i i x M M x K K )()(1010ελε+=+
(4)
ij j T i x M M x δε=+)(10
(5)
其中，i λ为新系统的特征值，i x 是相应的特征向量。
根据摄动理论，可将i λ、i x 按小参数ε展开成幂级数（因此，胡海昌院士称其为小参数法），即 )(32210ελεελλλO +++=i i i i
(6)
)(32210εεεO +++=i i i i x x x x
(7)
现在来确定1i λ，2i λ，1i x ，2i x ，将（6）、（7）式代入（4）式得
))()(())((2210102210221010i i i i i i i i i x x x M M x x x K K εεελεελλεεε+++++=+++
展开上式，略去)(3
εO ，并比较ε的同次幂系数可得 000000:i i i x M x K λε=
(8-1) 00101010001101:i i i i i i i i x M x M x M x K x K λλλε++=+
(8-2)
00201110111020011202:i i i i i i i i i i i i x M x M x M x M x M x K x K λλλλλε++++=+ (8-3)
将（6）、（7）式代入（5）式得
ij j j j T i i i x x x M M x x x δεεεεε=+++++))(()(2210102210
展开上式，并令j i =（从后面的分析可得，j i ≠的关系式在推导中没有被直接利用），比较ε的同次幂系数可得 1:0000=i T i x M x ε
(9-1) 0:010*******=++i T
i i T i i T i x M x x M x x M x ε
(9-2)
0:0021010111102002=++++i T
i i T i i T i i T i i T i x M x x M x x M x x M x x M x ε
(9-3)
至此，我们已推得了进行摄动分析（即求解i λ，i x 中的1i λ，2i λ，1i x ，2i x ）所需的
全部基本方程。实际上，(8-1)、(9-1)即为（1）、（2）式，是显然满足的。于是，可在原系统的特征值问题的基础上，通过(8-2)、(9-2)求解一阶摄动1i λ，1i x ，通过(8-3)、(9-3)求解二阶摄动2i λ，2i x ，代回（6）、（7）式，即求得i λ，i x ，且具有)(3εO 精度。
一阶摄动公式
可将1i x 表示成原系统模态（向量）02010,,,n x x x 的线性组合，即
∑==n
k k k i x x 1
011α
(10)
其中，1k α是n 个待定系数。
将（10）代入(8-2)，并左乘T k x 0，得
000101001
0100
0010
1
10
i T
k i i T k i n
k k k T k i i T
k n
k k k T
k x M x x M x x M x x K x x K x λλαλα
++=+∑∑==
利用正交规范关系（2）式，上式可简化成
000101000101001i T
k i i T k i i k i T k k k x M x x M x x K x λλλαλα++=+
(11)
k i = 时，00i k λλ=，1000=i T
k x M x ，由（11）式可得特征值的一阶摄动为
010101)(i i T i i x M K x λλ-=
(12)
k i ≠时，0000=i T k x M x ，则由（11）式可得
)()(0
00101
01
k i x M K x k i i i T
k k ≠--=λλλα
(13)
至此，n 个待定系数1k α中，只有1i α尚未确定，现在来求1i α。 用00M x T i 左乘（10）式两边，得
11
0001100
i n
k k T i k i T i x M x x M x αα==∑=
(14)
（14）式转置，并考虑到0M 对称，且1i α是一个常数，于是得
1001i i T i x M x α=
(15)
（14）、（15）式代入(9-2)可得
01012
1i T i i x M x -=α
(16)
由（13）、（16）两式即完全确定了)~1(1n k k =α，也就确定了（10）式的1i x ，而1i λ由
（12）式给出，于是可得一阶摄动公式为
??
???-??????--=-=∑≠=0010,10000101010
10101)(21)()(i i T i n i k k k i i i T k i i i T i i x x M x x x M K x x x M K x λλλλλ (17)
二阶摄动公式
为了得到更精确的摄动解，需要用到二阶摄动。根据展开定理，将2i x 按02010,,,n x x x 展开
∑==n
k k k i x x 1
022α
(18)
将（18）式代入(8-3)，并左乘T
k x 0，得
00201110111001
020001101
0200
)(i T k i i i i i i i T k n
k k k T k i i T
k n
k k k T k x M x x M x M x M x x M x x K x x K x λλλλαλα++++=+∑∑==
利用（2）式，上式可简化成
000
201110111000211002)(i T k i i i i i i i T
k i k i T
k k k x M x x M x M x M x x K x λλλλλαλα++++=
+ (19)
k i = 时，00i k λλ=，1000=i T
k x M x ，由（19）式可得特征值的二阶摄动为
)(0111011101102i i i i i i i T i i x M x M x M x K x λλλλ---=
(20)
k i ≠时，0000=i T k x M x ，则由（19）式可得
)()(0
00111011101102
k i x M x M x M x K x k i i i i i i i i T
k k ≠----=λλλλλα
(21)
现在来确定最后一个系数2i α，为此，用00M x T
i 左乘（18）式两边，得
21
0200
200
i n
k k k T i i T i x M x x M x αα==∑=
(22)
（22）式转置后，注意到0M 对称，2i α是一个数，可得
2002i i T i x M x α=
(23)
（22）、（23）式代入(9-3)可得
)(2
10111011102i T i i T i i T i i x M x x M x x M x ++-=α
(24)
于是就得到了二阶摄动公式
???
?
?????
++-??????----=---=∑≠=0011101110,100001110111011020111011101102)(21)()(i i T
i i T i i T i n
i k k k i i i i i i i i T k i i i i i i i i T i i x x M x x M x x M x x x M x M x M x K x x x M x M x M x K x λλλλλλλλλ
编程如下
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));
x2(:,m)=-0.5*(x0(:,m)'*M1*x1(:,m)+x1(:,m)'*M0*x1(:,m)+x1(:,m)'*M1 *x0(:,m))*x0(:,m);
for n=1:5
if n~=m
x2(:,m)=x2(:,m)+x0(:,n)'*(K1*x1(:,m)-r0(m,m)*M1*x1(:,m)-r1(m,m)*M 0*x1(:,m)-r1(m,m)*M1*x0(:,m))/(r0(m,m)-r0(n,n))*M0(:,n);
end
end
end
%%下面计算结构参数改变5%时的固有频率的结果
e=0.05;
[x,r]=eig(K0+e*K1,M0+e*M1);
w0=sqrt(r0)%%初始解
w=sqrt(r)%%精确解
w1=sqrt(r0+0.05*r1)%%一阶解
w2=sqrt(r0+0.05*r1+0.05^2*r2)%%二阶解
在这里只算了结构参数改变5%时的固有频率的结果，其结果为
w0 =
0.3129 0 0 0 0
0 0.9080 0 0 0
0 0 1.4142 0 0
0 0 0 1.7820 0
0 0 0 0 1.9754
w =
0.3022 0 0 0 0
0 0.8789 0 0 0
0 0 1.3732 0 0
0 0 0 1.7356 0
0 0 0 0 1.9273
w1 =
0.3018 0 0 0 0
0 0.8776 0 0 0
0 0 1.3711 0 0
0 0 0 1.7331 0
0 0 0 0 1.9249 w2 =
0.3024 0 0 0 0
0 0.8792 0 0 0
0 0 1.3733 0 0
0 0 0 1.7356 0
0 0 0 0 1.9274 下面给出本算例答案供参考
幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量 —一 .幂法 1. 幕法简介： 当矩阵A 满足一定条件时，在工程中可用幕法计算其主特征值 (按模最大) 及其特征向量。矩阵A 需要满足的条件为： ⑴I 1 I I 2|n |- 0, i 为A 的特征值 (2)存在n 个线性无关的特征向量，设为 X i ,X 2，…，X n 1.1计算过程： n 对任意向量x (0)，有x (0)八：-M —不全为0,则有 i 4 X (k 岀)=Ax (k)= = A k 岀乂。) n n A k 1 aq a 扌1 5 i =1 i =1 ■k 1 2 可见，当 1 — 1 越小时，收敛越快；且当k 充分大时，有 ? "1 2算法实现 ⑶.计算x Ay,… max(x); ⑷若| ?一十：；,输出-,y,否则,转(5) (5)若N ,置k 「k 1^ -,转3,否则输出失败信息，停机. 3 matlab 程序代码 (冲1 %叫 x (k 1) [x (k) k 二 u x (k) > (k+1) 1，对应的特征向量即是 x (1).输入矩阵A ，初始向量X ，误差限 最大迭代次数N (k) 0; y (k) max(abs(x (k ))
k=1; z=0; y=x0./max(abs(x0)); x=A*y; % z相当于■ %规范化初始向量%迭代格式 b=max(x); % b相当于： if abs(z-b)
第6章 求解线性代数方程组和计算矩阵特征值的迭代法
数值计算与MATLAB 1
《数值计算与MATLAB 》 第6章求解线性代数方程组和计算矩阵特征值的迭代法§1 求解线性代数方程组的迭代法 §2 方阵特征值和特征向量的计算 §3 矩阵一些特征参数的MATLAB计算
《数值计算与MATLAB 》 6.1 求解线性代数方程组的迭代法 1、迭代法的基本原理 如果线性方程组Ax=b的系数矩阵A非奇异，则方程组有唯一解。把这种方程中的方阵A分解成两个矩阵之差：A=C-D 若方阵C是非奇异的，把A它代入方程Ax=b中，得出 (C-D)x=b，两边左乘C-1，并令 M=C-1D，g= C-1b，移项可将方程Ax=b变换成： x=Mx+g 据此便可构造出迭代公式： x k+1 =Mx k+g， M=C-1D称为迭代矩阵。
《数值计算与MATLAB 》2. 雅可比（Jacobi）迭代法 如果方程组Ax=b的系数矩阵A非奇异，a ii ≠0，若可以把A 分解成： A=D-L-U=D+(-L)+(-U)， D=diag(a11,a22,…,a nn)； -L是严格下三角阵； -U是严格上三角矩阵； x= D-1((L+U)x +b)=D-1(L+U)x+ D-1b x k+1=D-1((L+U)x k+b)= D-1(L+U)x k + D-1b M M=D-1(L+U)称为雅可比迭代矩阵
《数值计算与MATLAB 》 ? ? ? ? ? ? ? ? ? ? ? ? = 6 7- 4 1 2 1- 2 6- 3- 1 1 5- 1 2 A ? ? ? ? ? ? ? ? ? ? ? ? = 6 1- 3- 2 D ? ? ? ? ? ? ? ? ? ? ? ? = 7 4- 1- 2- 1- L ? ? ? ? ? ? ? ? ? ? ? ? = 2- 6 1- 5 1- U M=D-1(L+U)= ? ? ? ? ? ? ? ? ? ? ? ? 7/6 2/3 - 1/6 - 2 2 2- 1/3 1/2 - 5/2 1/2 -
幂法求矩阵主特征值
!程序说明：幂法求矩阵主特征值 !日期：2010年11月30日 PROGRAM Matrix_EigenValue PARAMETER(N=3) REAL ARR(N,N) CALL INPUT(ARR,N) CALL MATEV(ARR,N) END PROGRAM SUBROUTINE INPUT(ARR,N) REAL ARR(N,N) OPEN(1,FILE='MAT.TXT') READ(1,*)((ARR(I,J),J=1,N),I=1,N) END SUBROUTINE SUBROUTINE MATEV(ARR,N) PARAMETER(EPS=1E-7) REAL :: ARR(N,N),X(N),X1(N),MAX=0 INTEGER :: K=0,P=0 X=RESHAPE((/1,1,1/),(/3/)) WRITE(1,*) ' 迭代次数 U（规范化向量） & & MAX(V)（主特征值）' DO WHILE(P/=N) WRITE(1,'(I6,A, ENDDO K=K+1 ENDDO END SUBROUTINE 输出结果： 1 1 0.5 1 1 0.25 0.5 0.25 2 迭代次数 U（规范化向量） MAX(V)（主特征值） 0 ( 1.000000 1.000000 1.000000 ) 0.000000 1 ( 0.909091 0.81818 2 1.000000 ) 2.750000 2 ( 0.837607 0.743590 1.000000 ) 2.659091 3 ( 0.799016 0.703035 1.000000 ) 2.604701 4 ( 0.77741 5 0.680338 1.000000 ) 2.575267 5 ( 0.765108 0.66740 6 1.000000 ) 2.558792 6 ( 0.758025 0.659963 1.000000 ) 2.549406 7 ( 0.753925 0.655655 1.000000 ) 2.544003 8 ( 0.751544 0.653153 1.000000 ) 2.540876 9 ( 0.750158 0.651697 1.000000 ) 2.539060 10 ( 0.749351 0.650848 1.000000 ) 2.538003 11 ( 0.748880 0.650354 1.000000 ) 2.537387 12 ( 0.748606 0.650065 1.000000 ) 2.537028 13 ( 0.748445 0.649897 1.000000 ) 2.536819 14 ( 0.748352 0.649799 1.000000 ) 2.536697 15 ( 0.748298 0.649741 1.000000 ) 2.536626 16 ( 0.748266 0.649708 1.000000 ) 2.536584 17 ( 0.748247 0.649688 1.000000 ) 2.536560 18 ( 0.748236 0.649677 1.000000 ) 2.536546 19 ( 0.748230 0.649670 1.000000 ) 2.536537 20 ( 0.748226 0.649667 1.000000 ) 2.536533 21 ( 0.748224 0.649664 1.000000 ) 2.536530 22 ( 0.748223 0.649663 1.000000 ) 2.536528 23 ( 0.748222 0.649662 1.000000 ) 2.536527 24 ( 0.748222 0.649662 1.000000 ) 2.536527 25 ( 0.748222 0.649662 1.000000 ) 2.536526 26 ( 0.748221 0.649661 1.000000 ) 2.536526 摘 要 根据现代控制理论课程的特点, 提出并利用MATLAB 设计了现代控制理论课程的实验, 给出了设计的每个实验的主要内容及使用到的MATLAB 函数, 并对其中的一个实验作了详细说明。通过这些实验, 将有助于学生理解理论知识, 学习利用MATLAB 解决现代控制理论问题。 关键词：现代控制理论、MATLAB 、仿真。 1设计目的、内容及要求 1.1设计目的 本课程设计以自动控制理论、现代控制理论、MATLAB 及应用等知识为基础，求连续系统对应的离散化的系统，并用计算系数阵按模最大的特征根法判别离散系统的稳定性，目的是使学生在现有的控制理论的基础上，学会用MATLAB 语言编写控制系统设计与分析的程序，通过上机实习加深对课堂所学知识的理解，掌握一种能方便地对系统进行离散化的实现和分析系统的稳定性的设计的工具。 1.2设计内容及要求 1 在理论上对连续系统离散化推导出算法和计算公式 2 画出计算机实现算法的框图 3 编写程序并调试和运行 4 以下面的系统为例，进行计算 ??????????----=041020122A ，?? ?? ? ?????=100B ，[]111-=c 5 分析运算结果 6 幂法迭代精度为ep=0.001，离散系统展开项数为20 7 程序应具有一定的通用性，对不同参数能有兼容性。 2算法选择及推导 2.1连续系统离散化算法 书P67离散化意义 已知被控对象的状态方程为： ()()()()()()t t u t y t t u t =+=+ x Ax B Cx D 对方程求解，得： 0()()0()()()o t t t t t t e t e u d τττ --=+?A A x x B 设0t kT =，(1)t k T =+，代入上式，得： H 公式 若省略T 则为{ ? +-++Φ=+T k kT d 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 T T AT T k kT T k A AT )()()()()()()(])1([(: )()(0 ,1,,)1()()()(0 )1(])1[(+==+=+Φ=====-=-+=?==Φ=???+-+τττττ τ 数值分析 幂法求矩阵A按模最大的特征值及其 特征向量 幂法的主要思想 设 n n ij R a A ?∈=)( ,其特征值为i λ ,对应特征向量为),,,1(n i x i =即 i i i x Ax λ= ),,1(n i =，且 x 1，······，x n 线性无关。求矩阵A 的主特征值及对应的特征向量。 幂法的基本思想: 任取一个非零初始向量 v 0 ∈R n 且v 0≠0， 由矩阵A 的乘幂构造一向量序列： 称{ v k }为迭代向量， A 特征值中 λ1为强占优，即▕ λ1▕＞▏λ2 ▏＞······＞▏λn ▏, {x 1，x 2，······，x n }线性无关，即{x 1，x 2，······，x n }为R n 中的一 个基，于是对任意的初始向量v 0 ∈R n 且 v 0≠0有展开式。 (v 0 用{x i } 的线性组合表示) （且设01≠α） 则 当k =2,3,… 时，v k = A v k-1 = A k v ? ?? 1Av v =0 212v A Av v ==01 1 v A Av v k k k ++==) ,,1,0(n k =∑==n i i i x v 1 α)(221101n n x x x A v A v ααα+++==n n x A x A x A ααα+++=2211n n n x x x λαλαλα+++=222111) (111 +≡x k αλk ε 其中 由假设▕ λ1▕＞▏λ2 ▏＞······＞▏λn ▏，得 ，从而 即,0lim =∞→k k ε且收敛速度由比值||12λλ=r 确定。 所以有 说明，当k 充分大时，有1 11 x v k k αλ≈，或 k k v 1λ 越来越接近特征 向量 规范化幂法的算法 ①输入矩阵 A 、初始向量v (0)，误差 eps ，实用中一般取 v (0)=（1，1，···，1）T ； ②k ←1； ③计算 v (k) ←Au (k-1)； ④m k ←max{ v (k) }，m k-1 ←{ v (k-1) }； ⑤u (k) ←v (k)/ m k ； ⑥如果▕ m k - m k-1▕＜eps ，则显示特征值λ1←和对应的特征 向量x (1)，终止； ⑦k=k+1，转③。 n k n n k k x x )()(1 2122λλαλλαε++=),,2(1||1 n i i =<λλ ),,,2(0)(lim 1n i k i k ==∞→λλ111 lim x v k k k αλ=∞ →。 11x α 实验十四特征值和特征向量 【实验目的】 1.了解特征值和特征向量的基本概念。 2.了解奇异值分解的基本概念。 3.学习、掌握MATLAB软件有关命令。 【实验内容】 计算特征值和特征向量 【实验准备】 1.特征值和特征向量的基本概念 A是n n?矩阵，如果λ满足Ax xλ =,则称λ是矩阵A的特征值，x 是矩阵A的特征向量。如果A是实对称矩阵，则特征值为实数，否则，特制值为复数。 2.矩阵的奇异值分解 3.矩阵特征值、奇异值分解的MATLAB命令 MATLAB中主要用eig求矩阵的特征值和特征向量，用svd求矩阵的奇异值分解。 eig(A)计算矩阵A的特征值 [X，D]=eig(A) D的对角线元素是特征值，X是矩阵，它的列是相应的特征向量。 s=svd(A)假设矩阵A的行数大于列数，则s是矩阵A的n个奇异值构成的向量。 [U，S，D]=svd(A)U,S,D为矩阵A的奇异值分解三对组。 【实验重点】 1.特征值与特征向量的计算 2.矩阵的奇异值分解 【实验难点】 1.矩阵的奇异值分解 【实验方法与步骤】 练习1求矩阵 31 13 A - ?? =?? -?? 的特征值和特征向量。 相应的MATLAB代码和计算结果为 A=[3-1；-1 3] A= 3 -1 -1 3 eig(A) %A的特征值 ans= 4 2 [X,D]=eig(A) %D的对角线元素是特征值，X是矩阵X= -0.7071 -0.7071 0.7071 -0.7071 D= 4 0 0 2 练习2求矩阵 23 45 84 A ?? ?? =?? ?? ?? 的奇异值分解。 相应的MATLAB代码和计算结果为 A=[2 3；4 5；8 4] A= 2 3 4 5 8 4 s=svd(A) %s是矩阵A的2个奇异值构成的向量s= 11.2889 2.5612 [U,S,V]=svd(A) %给出简洁方式的奇异值分解结果U= 0.3011 0.4694 -0.8301 0.5491 0.6263 0.5534 0.7796 -0.6224 -0.0692 S= 11.2889 0 0 2.5612 数学建模作业 计算机学院信计1102班姜圣涛 （1）幂法求矩阵最大特征值及特征向量： 程序为： #include while(1){ max=X[0]; for(i=0;i 用matlab求矩阵的特征值和特征向量 我要计算的矩阵： 1 3 5 1/3 1 3 1/5 1/3 1 [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/3 3/2 1 3/2 3/2 3/7 2/3 1 3/2 3/7 2/3 2/3 1 >> format rat >> A=[1 2/3 7/3 7/3 3/2 1 3/2 3/2 3/7 2/3 1 3/2 3/7 2/3 2/3 1] A = 1 2/3 7/3 7/3 3/2 1 3/2 3/2 3/7 2/3 1 3/2 3/7 2/3 2/3 1 >> [V,D]=eig(A) V = 1793/2855 504/3235 - 146/235i 504/3235 + 146/235i 1990/4773 670/1079 -3527/5220 -3527/5220 -509/959 4350/11989 1160/4499 + 287/3868i 1160/4499 - 287/3868i -350/647 838/2819 181/3874 + 1179/4852i 181/3874 - 1179/4852i 1238/2467 D = 810/197 0 0 0 0 -93/4229 + 455/674i 0 0 0 0 -93/4229 - 455/674i 0 0 0 0 -149/2201 ***************************************************************************************** 如何归一化求权重呢？ >> a=[1 3 5;1/3 1 3; 1/5 1/3 1] a = 1.0000 3.0000 5.0000 0.3333 1.0000 3.0000 0.2000 0.3333 1.0000 >> [V,D]=eig(a) V = 0.9161 0.9161 0.9161 0.3715 -0.1857 + 0.3217i -0.1857 - 0.3217i 0.1506 -0.0753 - 0.1304i -0.0753 + 0.1304i D = 幂法反幂法求解矩阵大小特征值及其对应的特征向量 ————————————————————————————————作者：————————————————————————————————日期： 数值计算解矩阵的按模最大最小特征值及对应的特征向量 一.幂法 1. 幂法简介： 当矩阵A 满足一定条件时，在工程中可用幂法计算其主特征值(按模最大)及其特征向量。矩阵A 需要满足的条件为： (1) 的特征值为A i n λλλλ,0||...||||21 ≥≥≥> (2) 存在n 个线性无关的特征向量，设为n x x x ,...,,21 1.1计算过程： i n i i i u x x αα,1 ) 0()0(∑==，有对任意向量不全为0，则有 1 11111221 12111 1 1 11 1 011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k n i i k i i n i i i k )(k (k))(k αλλλλλα++++=+=+++≈? ? ????+++======∑∑ 可见，当||1 2 λλ越小时，收敛越快；且当k 充分大时，有1)11 11)11111λαλαλ=??????==+++(k )(k k (k k )(k x x u x u x ，对应的特征向量即是)(k x 1+。 2 算法实现 . ,, 3,,1 , ).5() 5(,,,,||).4();max(,).3() (max(;0,1).2(,).1()() () (停机否则输出失败信息转置若转否则输出若计算最大迭代次数，误差限，初始向量输入矩阵βλβεβλβλε←+←<<-←←= ←←k k N k y x Ay x x abs x y k N x A k k k 3 matlab 程序代码 matlab矩阵运算和数组运算 作者：佚名教程来源：网络点击数：1368 更新时间：2010-5-3 矩阵运算和数组运算是Matlab的数值运算中的两大类运算。矩阵运算是按矩阵运算法 则进行的运算；数组运算无论是何种运算操作都是对元素逐个进行。 矩阵运算和数组运算指令对照汇总 矩阵运算指令指令含义数组运算指令指令含义 A' 矩阵转置 A.+B 对应元素相加 A+B 矩阵相加 A.-B 对应元素相减 A-B 矩阵相减 A.*B 同维数组对应元素相乘 s+B 标量加矩阵 s.*A A的每个元素乘s s-B,B-s 标量矩阵相减 A./B A的元素被B的对应元素除 A*B 矩阵相乘 B.\A 同上 A/B A右除B s./B, B.\s s 分别被B的元素除 B\A A左除B A.^n A的每个元素自乘n 次 inv(A) 矩阵求逆 log(A) 对A的每个元素求对数 A^n 矩阵的n次幂 sqrt(A) 对A的每个元素求平方根 f(A) 求A的各个元素的函数值 例： a=[1 2 3; 4 5 6; 7 8 9];b=[1 2 3; 3 2 1;1 4 5]; c=[1 1 1;2 3 1;1 0 2]; d=a*c^2+b d = 32 31 36 82 79 82 128 129 134 3.4 矩阵函数和数组函数 3.4.1 基本数组函数 数组函数是对各个元素的函数设计的。 f(.)基本函数表 函数名称功能函数名称功能 sin 正弦 acosh 反双曲余弦 cos 余弦 atanh 反双曲正切 tan 正切 acoth 反双曲余切 cot 余切 asech 反双曲正割 sec 正割 acsch 反双曲余割 csc 余割 fix 朝零方向取整 asin 反正弦 ceil 朝正无穷大方向取整 acos 反余弦 floor 朝负无穷大方向取整 atan 反正切 round 四舍五入到整数 atan2 四象反正切 rem 除后取余数 acot 反余切 sign 符号函数 asec 反正割 abs 绝对值 acsc 反余割 angle 复数相角 sinh 双曲正弦 imag 复数虚部 cosh 双曲余弦 real 复数实部 tanh 双曲正切 conj 复数共轭 coth 双曲余切 log10 常用对数 sech 双曲正割 log 自然对数 csch 双曲余割 exp 指数 asinh 反双曲正弦 aqrt 平方根 f(.)特殊函数表 函数名称功能函数名称功能 bessel 第一、第二类Bessel函数 erf 误差函数 一、矩阵的表示 在MATLAB中创建矩阵有以下规则： a、矩阵元素必须在”[ ]”内； b、矩阵的同行元素之间用空格（或”,”）隔开； c、矩阵的行与行之间用”;”（或回车符）隔开； d、矩阵的元素可以是数值、变量、表达式或函数； e、矩阵的尺寸不必预先定义。 二，矩阵的创建： 1、直接输入法 最简单的建立矩阵的方法是从键盘直接输入矩阵的元素，输入的方法按照上面的规则。建立向量的时候可以利用冒号表达式，冒号表达式可以产生一个行向量，一般格式是：e1:e2:e3，其中e1为初始值，e2为步长，e3为终止值。还可以用linspace函数产生行向量，其调用格式为：linspace(a,b,n) ，其中a和b是生成向量的第一个和最后一个元素，n是元素总数。 2、利用MATLAB函数创建矩阵 基本矩阵函数如下： (1) ones()函数：产生全为1的矩阵，ones(n)：产生n*n维的全1矩阵，ones(m,n)：产生m*n 维的全1矩阵； (2) zeros()函数：产生全为0的矩阵； (3) rand()函数：产生在（0，1）区间均匀分布的随机阵； (4) eye()函数：产生单位阵； (5) randn()函数：产生均值为0，方差为1的标准正态分布随机矩阵。 3、利用文件建立矩阵 当矩阵尺寸较大或为经常使用的数据矩阵，则可以将此矩阵保存为文件，在需要时直接将文件利用load命令调入工作环境中使用即可。同时可以利用命令reshape对调入的矩阵进行重排。reshape(A,m,n)，它在矩阵总元素保持不变的前提下，将矩阵A重新排成m*n的二维矩阵。 二、矩阵的简单操作 1．获取矩阵元素 可以通过下标（行列索引）引用矩阵的元素，如Matrix(m,n)。 也可以采用矩阵元素的序号来引用矩阵元素。 矩阵元素的序号就是相应元素在内存中的排列顺序。 在MATLAB中，矩阵元素按列存储。 序号(Index)与下标(Subscript )是一一对应的，以m*n矩阵A为例，矩阵元素A(i,j)的序号为(j-1)*m+i。 其相互转换关系也可利用sub2ind和ind2sub函数求得。 2．矩阵拆分 幂法求矩阵最大特征值 摘要 在物理、力学和工程技术中的很多问题在数学上都归结为求矩阵特征值的问题，而在某些工程、物理问题中，通常只需要求出矩阵的最大的特征值（即主特征值）和相应的特征向量，对于解这种特征值问题，运用幂法则可以有效的解决这个问题。 幂法是一种计算实矩阵A的最大特征值的一种迭代法，它最大的优点是方法简单。对于稀疏矩阵较合适，但有时收敛速度很慢。 用java来编写算法。这个程序主要分成了三个大部分：第一部分为将矩阵转化为线性方程组；第二部分为求特征向量的极大值；第三部分为求幂法函数块。其基本流程为幂法函数块通过调用将矩阵转化为线性方程组的方法，再经过一系列的验证和迭代得到结果。 关键词：幂法；矩阵最大特征值；j ava；迭代 POWER METHOD TO CALCULATE THE MAXIMUM EIGENV ALUE MATRIX ABSTRACT In physics, mechanics and engineering technology of a lot of problems in math boil down to matrix eigenvalue problem, and in some engineering, physical problems, usually only the largest eigenvalue of the matrix (i.e., the main characteristics of the value) and the corresponding eigenvectors, the eigenvalue problem for solution, using the power law can effectively solve the problem. Power method is A kind of computing the largest eigenvalue of real matrix A of an iterative method, its biggest advantage is simple.For sparse matrix is right, but sometimes very slow convergence speed. Using Java to write algorithms.This program is mainly divided into three most: the first part for matrix can be converted to linear equations;The second part is the eigenvector of the maximum;The third part is the exponentiation method of function block.Its basic process as a power law function block by calling the method of matrix can be converted to linear equations, then after a series of validation and iteration to get the results. Key words: Power method; Matrix eigenvalue; Java; The iteration Basic Matrix Operations 一、实验目的 1、掌握向量和矩阵的创建方法； 2、掌握向量和矩阵元素的索引方法； 3、掌握向量和矩阵的基本操作； 4、利用MATLAB编写程序进行矩阵运算。 二、基础知识 1、常见数学函数 函数名数学计算功能函数名数学计算功能 Abs(x) 实数的绝对值或复数的幅值floor(x) 对x朝-∞方向取整 Acos(x) 反余弦arcsin x gcd(m，n）求正整数m和n的最大公约数 acosh(x) 反双曲余弦arccosh x imag(x) 求复数x的虚部 angle(x) 在四象限内求复数 x 的相角lcm(m，n) 求正整数m和n的最小公倍数 asin(x) 反正弦arcsin x log(x) 自然对数（以e为底数） asinh(x) 反双曲正弦arcsinh x log10(x) 常用对数（以10为底数） atan(x) 反正切arctan x real(x) 求复数x的实部 atan2(x,y) 在四象限内求反正切Rem(m，n) 求正整数m和n的m/n之余数 atanh(x) 反双曲正切arctanh x round(x) 对x四舍五入到最接近的整数 ceil(x) 对x朝+∞方向取整sign(x) 符号函数：求出x的符号 conj(x) 求复数x的共轭复数sin(x) 正弦sin x cos(x) 余弦cos x sinh(x) 反双曲正弦sinh x cosh(x) 双曲余弦cosh x sqrt(x) 求实数x的平方根：x exp(x) 指数函数xe tan(x) 正切tan x fix(x) 对x朝原点方向取整tanh(x) 双曲正切tanh x 2、常量与变量 系统的变量命名规则：变量名区分字母大小写；变量名必须以字母打头，其后可以是任意字母，数字，或下划线的组合。此外，系统内部预先定义了几个有特殊意义和用途的变量，见下表： 特殊的变量、常量取值 C语言课程设计报告 课程名称：计算机综合课程设计 学院：土木工程学院 设计题目：矩阵特征值分解 级别： B 学生姓名： 学号： 同组学生：无 学号：无 指导教师： 2012年 9 月 5 日 C语言课程设计任务书 （以下要求需写入设计报告书） 学生选题说明: 以所发课程设计要求为准，请同学们仔细阅读； 本任务书提供的设计案例仅供选题参考；也可自选，但难易程度需难度相当； 鼓励结合本专业（土木工程、力学）知识进行选题，编制程序解决专业实际问题。 限2人选的题目可由1-2人完成（A级）；限1人选的题目只能由1人单独完成（B级）；设计总体要求: 采用模块化程序设计； 鼓励可视化编程； 源程序中应有足够的注释； 学生可自行增加新功能模块（视情况可另外加分）； 必须上机调试通过； 注重算法运用，优化存储效率与运算效率； 需提交源程序(含有注释)及相关文件(数据或数据库文件)； (cpp文件、txt或dat文件等) 提交设计报告书，具体要求见以下说明。 设计报告格式： 目录 1．课程设计任务书（功能简介、课程设计要求）； 2．系统设计（包括总体结构、模块、功能等，辅以程序设计组成框图、流程图解释）； 3．模块设计（主要模块功能、源代码、注释（如函数功能、入口及出口参数说明，函数调用关系描述等）； 4．调试及测试：（调试方法，测试结果的分析与讨论，截屏、正确性分析）； 5．设计总结：（编程中遇到的问题及解决方法）； 6．心得体会及致谢； 参考文献 1.课程设计任务书 功能简介： a）输入一个对称正方矩阵A，从文本文件读入； b）对矩阵A进行特征值分解，将分解结果：即U矩阵、S矩阵输出至文本文件； c）将最小特征值及对应的特征向量输出至文本文件； d）验证其分解结果是否正确。 提示：A=USU T,具体算法可参考相关文献。 功能说明： 矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。 注：以三阶对称矩阵为例 2.系统设计 3.模块设计 #include 西京学院数学软件实验任务书 课程名称数学软件实验班级数0901 学号0912020107 姓名李亚强 实验课题矩阵特征值及相应特征向量的Jacobi法，QR法 实验目的熟悉矩阵特征值及相应特征向量的Jacobi法，QR法运用Matlab/C/C++/Java/Maple/Mathematica等其中实验要求 一种语言完成 实验内容矩阵特征值及相应特征向量的Jacobi法，QR法成绩教师 实验十三实验报告 一、实验名称：矩阵特征值及相应特征向量的Jacobi法，QR法。 二、实验目的：熟悉矩阵特征值及相应特征向量的Jacobi法，QR 法。 三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成程序设计。 四、实验内容： %矩阵特征值及相应特征向量的Jacobi法 function [D,R]=Jacobi(A,eps) if nargin==2 eps=1.0e-5; end n=length(A); R=eye(n); while 1 Amax=0; for l=1:n-1 for k=l+1:n if abs(A(l,k))>Amax Amax=abs(A(l,k)); i=l;j=k; end end end if Amax t=1; else t=sign(d)/(abs(d)+sqrt(d^2+1)); end c=1/sqrt(t^2+1); s=c*t; for l=1:n if l==i Aii=A(i,i)*c^2+A(j,j)*s^2+2*A(i,j)*s*c; Ajj=A(i,i)*s^2+A(j,j)*c^2-2*A(i,j)*s*c; A(i,j)=(A(j,j)-A(i,i))*s*c+A(i,j)*(c^2-s^2); A(j,i)=A(i,j); A(i,i)=Aii; A(j,j)=Ajj; elseif l~=j Ail=A(i,l)*c+A(j,l)*s; Ajl=-A(i,l)*s+A(j,l)*c; A(i,l)=Ail;A(l,i)=Ail; A(j,l)=Ajl;A(l,j)=Ajl; end Rli=R(l,i)*c+R(l,j)*s; Rlj=-R(l,i)*s+R(l,j)*c; R(l,i)=Rli; R(l,j)=Rlj; end end D=diag(diag(A)); %矩阵特征值及相应特征向量的QR法 function l=qrtz(A,M) for(i=1:M) [q,r]=qr(A); A=r*q; l=diag(A); end 数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序 矩阵的特征值与特征向量的计算 摘要 物理，力学，工程技术中的很多问题在数学上都归结于求矩阵特征值的问题，例如振动问题（桥梁的振动，机械的振动，电磁振动等）、物理学中某些临界值的确定问题以及理论物理中的一些问题。矩阵特征值的计算在矩阵计算中是一个很重要的部分，本文使用幂法和反幂法分别求矩阵的按模最大，按模最小特征向量及对应的特征值。 幂法是一种计算矩阵主特征值的一种迭代法，它最大的优点是方法简单，对于稀疏矩阵比较合适，但有时收敛速度很慢。其基本思想是任取一个非零的初始向量。由所求矩阵构造一向量序列。再经过所构造的向量序列求出特征值和特征向量。 反幂法用来计算矩阵按模最小特征向量及其特征值，及计算对应于一个给定近似特征值的特征向量。本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。然后经过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。 关键词：矩阵；特征值；特征向量；冥法；反冥法 THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIX ABSTRACT Physics, 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 求权重程序 A=[1,1/3,3,1;3,1,7,3;1/3,1/7,1,1/5;1,1/3,5,1]; %求得x为特征向量矩阵，y为特征值矩阵[x,y]=eig(A); %找到y中对应最大的特征值所在列m [maxy,m]=find(y==max(max(y))); %w即为矩阵A的权重！ w=x(:,m)/sum(x(:,m)) 结果： w= 0.1881 0.5347 0.0597 0.2175 求最大特征值 b=[1,1/3,3,1;3,1,7,3;1/3,1/7,1,1/5;1,1/3,5,1] [v,d]=eig(b) r=abs(sum(d)); n=find(r==max(r)); max_d_b=d(n,n)%最大特征根 max_v_b=v(:,n)%最大特征根所对应的特征向量 结果： b= 1.00000.3333 3.0000 1.0000 3.0000 1.00007.0000 3.0000 0.33330.1429 1.00000.2000 1.00000.3333 5.0000 1.0000 v= 0.3083-0.4859-0.1383-0.0804i-0.1383+0.0804i 0.87650.8503-0.8693-0.8693 0.09790.00000.0514+0.1126i0.0514-0.1126i 0.35650.20240.3162-0.3216i0.3162+0.3216i d= 4.0571000 0-0.000000 00-0.0285+0.4804i0 000-0.0285-0.4804i max_d_b=4.0571最大特征根 max_v_b=最大特征根所对应的特征向量 0.3083 0.8765 0.0979 0.3565matlab求矩阵特征值特征向量 乘幂法
幂法求矩阵A按模最大的特征值及其特征向量
matlab实验十四__特征值和特征向量
数学建模 用幂法 和法 根法求特征值特征向量
用MATLAB求矩阵特征值
幂法反幂法求解矩阵大小特征值及其对应的特征向量
matlab矩阵运算和数组运算
矩阵在matlab中的基本命令
幂法求矩阵最大特征值
MATLAB基本矩阵运算
雅克比法求矩阵特征值特征向量
数学实验“矩阵特征值及相应特征向量的Jacobi法,QR法”实验报告(内含matlab程序)
数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序
求判断矩阵权重以及最大特征值 MATLAB程序