结构动力学大作业对于如下结构，是研究质量块的质量变化和在简支梁上位置的变化对整个系统模态的影响。
1以上为一个简支梁结构。
集中质量块放于梁上，质量块距简支梁的左端点距离为L.将该简支梁简化为欧拉伯努利梁，并离散为N 个单元。
每个单元有两个节点，四个自由度。
单元的节点位移可表示为：]1122,,,e v v δθθ⎡=⎣则单元内一点的挠度可计作：带入边界条件：1332210)(x a x a x a a x v +++=01)0(a v x v ===3322102)(L a L a L a a v L x v +++===110d d a x vx ===θ2321232d d L a L a a xv Lx ++===θ10v a =[]1234N N N N N =建立了单元位移模式后，其动能势能均可用节点位移表示。
单元的动能为：00111()222l l T T Tke e e e e y E dx qN Ndxq q mq t ρρ∂===∂⎰⎰ 其中m 为单元质量阵，并有：lT m N Ndx ρ=⎰带入公式后积分可得：2222156225413224133541315622420133224l l l l l l l m l l l l l l ρ-⎡⎤⎢⎥-⎢⎥=⎢⎥-⎢⎥---⎣⎦单元势能可表示为2220011()()222Tl lTT e pe e e e q y E EI dx EI N N dxq q Kq x ∂''''===∂⎰⎰ 其中K 为单元刚度矩阵，并有()lT K EI N N dx ''''=⎰2232212612664621261266264l l l l l l EI k l l l l l l l -⎡⎤⎢⎥-⎢⎥=⎢⎥---⎢⎥-⎣⎦以上为单元类型矩阵，通过定义全局位移矩阵，可以得到系统刚度矩阵和系统质量矩11θ=a )2(1)(3211222θθ+--=Lv v L a )(1)(22122133θθ++-=Lv v L a 1232133222231)(θ⎪⎪⎭⎫⎝⎛+-+⎪⎪⎭⎫ ⎝⎛+-=L x L x x v L x L x x v 22232332223θ⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎭⎫ ⎝⎛-+L x L x v L x L x 24231211)()()()()(θθx N v x N x N v x N x v +++=阵。
当集中质量块加到简支梁上时，可以认为系统的刚度矩阵不变，而质量矩阵在 有基础上有所增加。
当集中质量块位于第i 个单元内，距单元左节点的位移为0x ，则其动能可表示为：000()()2T Te e m V N x N x δδ=则质量块附加到总体刚度矩阵的子矩阵尾为：000()()T m N x N x例如当质量为0m 的质量块位于单元的右节点上时，其质量阵可表示为00000000000100000m m ⎡⎤⎢⎥⎢⎥'=⎢⎥⎢⎥⎣⎦根据以上理论，编写matlab 程序。
（相关参数：密度7860KG/M3， 长度L=1m ，截面尺寸0.02m*0.02m ）以下为matlab 源代码： 主程序clear; clc;Beam_InputData541; %输入相关参数，划分了40个单 %元k=zeros(No_nel*No_dof,No_nel*No_dof); %初始化单元刚度阵和质量阵 m=zeros(No_nel*No_dof,No_nel*No_dof); kk=zeros(Sys_dof,Sys_dof);mm=zeros(Sys_dof,Sys_dof); %初始化总体质量阵和刚度阵index=zeros(No_nel*No_dof,1); foriel=1:No_el nd(1)=iel; nd(2)=iel+1; leng=0.025;[k,m]=BeamElement11(prop,leng); %调用子函数形成单元刚度矩阵 %质量矩阵index=femEldof(nd,No_nel,No_dof);kk=femAssemble1(kk,k,index);mm=femAssemble1(mm,m,index); %调用子函数形成总体刚度矩阵 end %和质量矩阵Mzhiliangkuai=1/5*prop(3)*prop(6)*1; %将集中质量块加进总体质量矩Ms=zeros(82,82); %之中Ms(41,41)=Mzhiliangkuai; %此为质量块为于梁中央时的处mm=Ms+mm; 理方法kk(1,:)=[]; %引入边界条件，划去相应的kk(80,:)=[]; %行和列kk(:,1)=[];kk(:,80)=[];mm(1,:)=[];mm(80,:)=[];mm(:,1)=[];mm(:,80)=[];[V,D]=eig(kk,mm); %求解特征值和特征向量%或使用雅克比迭代求解特征值与特征向量，以下为代码%function[aa,v]=jac(a)%此程序用jacobi方法求实对称矩阵的全部特征值和特征向量%输入x：n x n矩阵%----------------------------------------------------%% n=length(a);% aa=zeros(1,n);% u=zeros(n,n);% l=0;% v是n阶单位矩阵% v=eye(n,n);% while (1)% l=l+1;% fm=0.0;% for i=1:n% for j=1:n% if i==j% continue;% else% d=abs(a(i,j));% if d>fm% fm=d;% p=i;% q=j;% end% end% end% end% if fm<0.0001% break;% end% x=2*a(p,q)*(sign(a(p,p)-a(q,q)));% y=abs(a(p,p)-a(q,q));% if y~=0% cn=sqrt((1+y/sqrt(x^2+y^2))/2.0);% else% cn=1/sqrt(2);% end% if y~=0% sn=x/sqrt(x^2+y^2)/(2*cn);% elseif y==0 & a(p,q)>0% sn=1/sqrt(2);% elseif y==0 & a(p,q)<0% sn=-1/sqrt(2);% end% %更新a% a_t=a;% for i=1:n% if i~=p & i~=q% a(i,p)=a_t(i,p)*cn+a_t(i,q)*sn;% a(i,q)=-a_t(i,p)*sn+a_t(i,q)*cn;% a(p,i)=a(i,p);% a(q,i)=a(i,q);% end% end% a(p,p)=a_t(p,p)*cn*cn+2*a_t(p,q)*cn*sn+a_t(q,q)*sn*sn; % a(q,q)=a_t(p,p)*sn*sn-2*a_t(p,q)*sn*cn+a_t(q,q)*cn*cn; % a(p,q)=0.0;% a(q,p)=a(p,q);% %更新v% v_t=v;% for i=1:n% v(i,p)=v_t(i,p)*cn+v_t(i,q)*sn;% v(i,q)=-v_t(i,p)*sn+v_t(i,q)*cn;% end% end% for i=1:n% aa(1,i)=a(i,i);% end%end%for i=1:10;% w(i,1)=sqrt(D(i,i))/(2*pi);%end%E=prop(1); Iz=prop(8); rho=prop(3)*prop(6); L=1;%i=(1:10)';%omega2=i.*i*pi^2*sqrt(E*Iz/(rho*L^4));%omega3=omega2/(2*pi);%for i=1:39;% x1(i)=i*0.025;% y1(i)=V(2*i,3);%end%plot(x1,y1);各子函数：Beam_InputData541No_el=40; number of elementsNo_nel=2; No_dof=2;No_node=(No_nel-1)*No_el+1;Sys_dof=No_node*No_dof;prop(1)=2.1e11; % elastic modulusprop(3)=7860; prop(4)=0.02;prop(5)=0.02; prop(6)=0.02*0.02;prop(7)=0.02*0.02^3/12; prop(8)=0.02*0.02^3/12;BeamElement11function [k,m]=BeamElement11(prop,leng);%,Opt_mass)E=prop(1);u=prop(2);rho=prop(3);A=prop(6); Iz=prop(8);G=E/(2*(1+u));c=E*Iz/(leng^3);k0=[12 6*leng -12 6*leng;6*leng 4*leng^2 -6*leng 2*leng^2;-12 -6*leng 12 -6*leng;6*leng 2*leng^2 -6*leng 4*leng^2];k=c*k0;mass=rho*A*leng;m0=[156 22*leng 54 -13*leng;22*leng 4*leng^2 13*leng -3*leng^2;54 13*leng 156 -22*leng;-13*leng -3*leng^2 -22*leng 4*leng^2];m=mass/420*m0;femAssemble1function [kk]=femAssemble1(kk,k,index)eldof = length(index);for i=1:eldofii=index(i);for j=1:eldofjj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);endendfunction [index]=femEldof(nd,No_nel,No_dof)k=0;for i=1:No_nelstart = (nd(i)-1)*No_dof;for j=1:No_dofk=k+1;index(k)=start+j;endend首先我分析了无集中质量块时梁的模态特性，提取了前十阶固有频率和前四阶振型。
然后分析了当集中质量块分别为简支梁质量1/40,1/20.1/10,1/5,1/2,1 倍和2倍时这八种工况。
