振动理论第次大作业.docVIP

中国科学技术大学 本科实验报告 课程名称 振动理论 实验名称 第二次大作业 姓 名 学 院 工程科学学院 系 别 近代力学系 专 业 理论与应用力学 年 级 学 号 任课教师 陈海波 2014年 12 月 24日 《振动理论》第二次大作业报告 流程图如下： matlab程序代码如下： %*****************龙格库塔法求解多自由度楼层振动问题**************** %K：刚度矩阵 M：质量矩阵 eig_vector:inv(M)*K的特征向量也就是模态矩阵 %k2: eig_vector*K*eig_vector为正交化刚度矩阵 %m2: eig_vector*M*eig_vector为正交化质量矩阵 %x存储q,y存储q的导数 %weiyi:保存楼层的位移随时间变化的值 %****************************************************************** M=0,0,0,0;00,0,0;0,00,0;0,0,00;0,0,0,0,7880000]; K=[35000000000.0000,0000,0,0,0; 0000,31500000000.0000,-14000000000.0000,0,0; 0,-14000000000.0000,28000000000.0000,-14000000000.0000,0; 0,0,-14000000000.0000,24500000000.0000,-10500000000.0000; 0,0,0,-10500000000.0000,10500000000.0000]; A=M\K; [eig_vector,namda]=eig(A); omg=zeros(5,1); for i=1:5 omg(i)=sqrt(namda(i,i)); eig_vector(:,i)=eig_vector(:,i)/eig_vector(5,i); end m2= eig_vector*M*eig_vector; k2= eig_vector*K*eig_vector; for i=1:5 mm(i)=m2(i,i); kk(i)=k2(i,i); end %初始化，输入求解区间最大值t_end，和求解步长h t_start=0; t_end=input(请输入要求解区间的最大值（最小值已默认置为0）：); h=input(请输入差分计算是的时间步长h:); n=(t_end-t_start)/h; x1=zeros(n+1,5); y=zeros(n+1,5);%存储x的导数的值 t=t_start:h:t_end;%生成t t=t; %下面是微分方程的初值 x1(1,:)=0; y(1,:)=0; %下面进行循环求解 for j=1:5 for i=1:n T1=t(i); T2=t(i)+h/2; T3=t(i)+h/2; T4=t(i)+h; X1=x1(i,j); Y1=y(i,j);F1=f(T1,X1,Y1,j,omg(j),mm(j),eig_vector); X2=x1(i,j)+Y1*h/2; Y2=y(i,j)+F1*h/2; F2=f(T2,X2,Y2,j,omg(j),mm(j),eig_vector); X3=x1(i,j)+Y2*h/2;Y3=y(i,j)+F2*h/2;F3=f(T3,X3,Y3,j,omg(j),mm(j),eig_vector); X4=x1(i,j)+Y3*h;Y4=y(i,j)+F3*h;F4=f(T4,X4,Y4,j,omg(j),mm(j),eig_vector); x1(i+1,j)=x1(i,j)+(h/6)*(Y1+2*Y2+2*Y3+Y4); y(i+1,j)=y(i,j)+(h/6)*(F1+2*F2+2*F3+F4); end end %*****************龙格库塔法求解**************** %******************下面进行绘图****************************