MATLP运用.docVIP

卢彬彬 工程力学一阶摄动公式： 二阶摄动公式： 新系统的值： 取质量阵，刚度阵 计算原系统特征值和特征向量 m0=[1 0 0;0 1 0;0 0 2]; k0=[2 -1 0;-1 3 -2;0 -2 2]; [x0,a0]=eig(k0,m0) x0 = 0.2691 -0.8782 -0.3954 0.5008 -0.2231 0.8363 0.5817 0.2992 -0.2685 a0 = 0.1392 0 0 0 1.7459 0 0 0 4.1149 原系统的特征值为孤立特征值 在这种情况下，当、很小时，新系统的特征值及相应的均只有小变化。根据摄动理论，可将、按小参数展开成幂级数 现在来确定，，， 取 根据一阶摄动公式求， m1=[1 -1 0;-1 2 1;0 1 2]; k1=[1 0 2;0 1 -1;2 -1 1]; for i=1:3 a1(i,i)=x0(:,i)*(k1-a0(i,i)*m1)*x0(:,i); end a1 a1 = 0.4875 0 0 0 -0.9223 0 0 0 -6.0652 for i=1:3 for b=1:3 if i~=b s1=0; y1=(x0(:,b)*(k1-a0(i,i)*m1)*x0(:,i))/(a0(i,i)-a0(b,b))*x0(:,b); s1=s1+y1; end end end for i=1:3 x1(:,i)=s1-0.5*(x0(:,i)*m1*x0(:,i))*x0(:,i); end x1 = 0.9617 1.4024 1.5502 -0.0937 0.3564 -0.5016 -0.8542 -0.4777 -0.1427 for i=1:3 a2(i,i)=x0(:,i)*(k1*x1(:,i)-a0(i,i)*m1*x1(:,i)-a1(i,i)*m0*x1(:,i)-a1(i,i)*m1*x0(:,i)); end a2 = 0.7213 0 0 0 0.4836 0 0 0 15.8159 根据二阶摄动公式求， for i=1:3 a2(i,i)=x0(:,i)*(k1*x1(:,i)-a0(i,i)*m1*x1(:,i)-a1(i,i)*m0*x1(:,i)-a1(i,i)*m1*x0(:,i)); end a2 = 0.7213 0 0 0 0.4836 0 0 0 15.8159 for i=1:3 for b=1:3 if i~=b s2=0; y2=(x0(:,b)*(k1*x1(:,i)-a0(i,i)*m1*x1(:,i)-a1(i,i)*m0*x1(:,i)-a1(i,i)*m1*x0(:,i))/(a0(i,i)-a0(b,b)))*x0(:,b); s2=s2+y2; end end end for i=1:3 x2(:,i)=s2-0.5*(x0(:,i)*m1*x1(:,i)+x1(:,i)*m0*x1(:,i)+x1(:,i)*m1*x0(:,i))*x0(:,i); end x2 = -0.7673 -0.5366 -1.5170 0.0519 -0.1363 1.0265 0.6459 0.1828 -0.0909 1.当=0.01时， 一阶摄动得到的特征值和特征向量 e=0.01 y1=a0+e*a1 y1 = 0.1441 0 0 0 1.7367 0 0 0 4.0543 t1=x0+e*x1 t 1= 0.2787 -0.8642 -0.3799 0.4998 -0.2196 0.8313 0.5732 0.2944 -0.2699 二阶摄动得到的特征值和特征向量 2=a0+e*a1+e^2*a2 y2 = 0.1441 0 0