高手们就帮帮我吧,确实不会啊,纠结一下午了
-
- function vtb5(tf,delt)
- close all;clc
- fid1=fopen('disp','wt');
- m=2*[1 0 0;0 1 0;0 0 1];
- c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
- k=50*[2 -1 0;-1 2 -1;0 -1 2];
- x0=[1 1 1]';
- v0=[1 1 1]';
- bita=1/6;
- md=m+delt/2*c+bita*delt^2*k;
- [E,F]=eig(m\k);
- diag(sqrt(F));
- for t=0:delt:tf;
- f=[2.00*sin(3.754*t) -2.00*cos(2.2*t) 1.00*sin(2.8*t)]';
- if t==0;xdd0=m\(f-k*x0-c*v0);
- else
- xdd=md\(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0));
- x=md\(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);
- xd=v0+delt/2*(xdd0+xdd);
- xdd0=xdd;v0=xd;x0=x;
- fprintf(fid1,'%10.4f',x0);
- end
- end
- fid2=fopen('disp','rt');
- n=tf/delt;
- x=fscanf(fid2,'%f',[3,n]);
- t=1:n;
- figure('numbertitle','off','name','自由度1的位移','pos',[450 180 400 420]);
- plot(t,x(1,t)),grid,xlabel('时间*0.1秒'),title('自由度1的位移与时间的关系')
- figure('numbertitle','off','name','自由度2的位移','pos',[350 160 400 420]);
- plot(t,x(2,t)),grid,xlabel('时间*0.1秒'),title('自由度2的位移与时间的关系')
- figure('numbertitle','off','name','自由度3的位移','pos',[250 140 400 420]);
- plot(t,x(3,t)),grid,xlabel('时间*0.1秒'),title('自由度3的位移与时间的关系')
复制代码
[ 本帖最后由 ChaChing 于 2010-3-24 23:29 编辑 ] |