matlab 摄动,孤立特征值情况的矩阵摄动法算例(matlab编程).doc

孤立特征值情况的矩阵摄动法算例(matlab编程).doc

设原始特征值问题是(1))~1(00nixMKiii相应的正交规范条件是(2)),(0jixijjTi先考虑一个五个自由度的系统．其中，、0K分别是原系统的阶对称质量阵、刚度阵，是特征值，且Mn0i，是固有频率，是相应的特征向量。20iii0ix010.5M02121K系统结构参数改变后，相应的质量阵、刚度阵均有相应的变化，设系统修改后的质量阵、刚度阵分别为(3)10KM其中取110.510K显然，新系统(即结构修改后的系统或称摄动系统)的特征值问题及相应的正交规范条件是(4)iiixMxK)()(1010(5)ijjTi10其中，为新系统的特征值，是相应的特征向量。iix根据摄动理论，可将、按小参数展开成幂级数(因此，胡海昌院士称其为小参数法)ii，即(6))(3210iiii(7)iiiixx现在来确定，，，，将(6)、(7)式代入(4)式得1i2i1i2i))()(())((2101021000iiiiiiiiixxMxxK展开上式，略去，并比较的同次幂系数可得3(8-1)000:iiixK(8-2)0111101:iiiiiiiixx(8-3)0200201202:iiiiiiiiiiiiMMxMKx将(6)、(7)式代入(5)式得ijjjjTiiixx))(()(21010210展开上式，并令(从后面的分析可得，的关系式在推导中没有被直接利用)，比ji较的同次幂系数可得(9-1)1:00iTixM(9-2)0:101iTiiiiTix(9-3):210202iTiiiiTiiiiTixxM至此，我们已推得了进行摄动分析(即求解，中的，，，)所需的iix1i2i1ix2i全部基本方程。实际上，(8-1)、(9-1)即为(1)、(2)式，是显然满足的。于是，可在原系统的特征值问题的基础上，通过(8-2)、(9-2)求解一阶摄动，，通过(8-3)、(9-3)求解1ii二阶摄动，，代回(6)、(7)式，即求得，，且具有精度。2iixiix)(3一阶摄动公式可将表示成原系统模态(向量)的线性组合，即1ix0210,,nx(10)kki10其中，是n个待定系数。1k将(10)代入(8-2)，并左乘，得Tkx0010101001100iTkiiTkinkkTiiknkTxMxxMKKx利用正交规范关系(2)式，上式可简化成(11)0101010101iTkiiTkiikiTkkx时，，，由(11)式可得特征值的一阶摄动为iii(12)01101)(iiTiixMKx时，，则由(11)式可得ki0iTkxM(13))()(01101kixkiiTk至此，n个待定系数中，只有尚未确定，现在来求。1k1i1i用左乘(10)式两边，得0MxTi(14)1100inkkTiiTixMx(14)式转置，并考虑到对称，且是一个常数，于是得0i(15)101iiTix(14)、(15)式代入(9-2)可得(16)0112iTiixM由(13)、(16)两式即完全确定了，也就确定了(10)式的，而)~(nk1ix由(12)式给出，于是可得一阶摄动公式为1i(17)010,1001100)(2)(iTinikkkiiTkiiiiixMxKx二阶摄动公式为了得到更精确的摄动解，需要用到二阶摄动。根据展开定理，将按展2ix0210,,nx开(18)nkkixx1022将(18)式代入(8-3)，并左乘，得Tk00201101010201)(iTkiiiiiiiTknkTiikkxMxMxxMxKK利用(2)式，上式可简化成(19)0201101002)(iTkiiiiiiiTkikxxx时，，，由(19)式可得特征值的二阶摄动为iii(20))(011010102iiiiiiiTiixMxKx时，，则由(19)式可得ki0iTkM(21))()(00111012kixxxkiiiiiiik现在来确定最后一个系数，为此，用左乘(18)式两边，得2iMTi(22)2102020inkkiiTixx(22)式转置后，注意到对称，是一个数，可得0M2i(23)20iiTix(22)、(23)式代入(9-3)可得(24))(21011010iTiiTiiTiixM于是就得到了二阶摄动公式011010,1011010102)()(iTiiTiiTinikkkiiiiiiiTkiiiiiiiiiixMxxKxx编程如下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)producesadiagonalmatrixDofgeneralizedeigenvaluesandafullmatrixVwhosecolumnsarethecorrespondingeigenvectorssothatA*V=B*V*D.%%x0为特征向量矩阵，r0为特征值对角阵%%下面来计算一阶摄动=1:5r1(m,m)=x0(:,m) *(K1-r0(m,m)*M1)*x0(:,m);x1(:,m)=-0.5*x0(:,m) *M1*x0(:,m)*x0(:,m);forn=1:5ifn~=mx1(:,m)=x1(:,m)+x0(:,n) *(K1-r0(m,m)*M1)*x0(:,m)/(r0(m,m)-r0(n,n))*M0(:,n);endendend%%再来计算二阶摄动=1:5r2(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);forn=1:5ifn~=mx2(:,m)=x2(:,m)+x0(:,n) *(K1*x1(:,m)-r0(m,m)*M1*x1(:,m)-r1(m,m)*M0*x1(:,m)-r1(m,m)*M1*x0(:,m))/(r0(m,m)-r0(n,n))*M0(:,n);endendend%%下面计算结构参数改变5%时的固有频率的结果e=0.05;[x,r]=eig(K0+e*K1,M0+e*M1);w0=sqrt(r0)%%初始解w=sqrt(r)%%精确解