【求助，matlab动力学编程出错】
%--------------------------------------------------------
%Ex.2.2 plane 9-bar truss
%--------------------------------------------------------
%To compute the dynamic response of node 4 on the vertical direction
%
%
E=2.0e11;
A=2.5e-3;
density=7860;
node_number=6;
elment_number=9;
nc=; %node_coordinate
en=; % element_node
ed(1:node_number,1:2)=1; % elment_displacement
dof=0;
for loopi=1:node_number
for loopj=1:2
dof=dof+1;
ed(loopi,loopj)=dof;
end
end
ek=E*A*; % elment_stiffness_matrix
em=(density*A)/6*[2 0 1 0; ... %elment_mass_matrix
0 2 0 1; ...
1 0 2 0; ...
0 1 0 2];
k(1:dof,1:dof)=0; %tructural_stiffness_matrix
m=k; % structural_mass_matrix, same sizewithk
theta(1:9)=0;
el(1:9)=0;
e2s(1:4)=0; %index of transform the elment displament number to structural
for loopi=1:elment_number
for zi=1:2
e2s((zi-1)*2+1)=ed(en(loopi,zi),1);
e2s((zi-1)*2+2)=ed(en(loopi,zi),2);
end
el(loopi)=sqrt((nc(en(loopi,1),1)-nc(en(loopi,2),1))^2+(nc(en(loopi,1),2)-nc(en(loopi,2),2))^2);
theta(loopi)=asin((nc(en(loopi,1),2)-nc(en(loopi,2),2))/el(loopi));
lmd=;
t=;
dk=t'*ek*t/el(loopi);
dm=t'*em*t*el(loopi);
for jx=1:4
for jy=1:4
k(e2s(jx),e2s(jy))=k(e2s(jx),e2s(jy))+dk(jx,jy);
m(e2s(jx),e2s(jy))=m(e2s(jx),e2s(jy))+dm(jx,jy);
end
end
end
c=0*k+0*m;
nt=1500;dt=0.0001;
time=0:dt:nt*dt;
q0=zeros(dof,1);
dq0=zeros(dof,1);
bcdof=zeros(dof,1);
fd=zeros(dof,nt);
for i=1:nt
fd(6,i)=200;
end
=TransResp4(k,c,m,fd,bcdof,nt,dt,q0,dq0);
plot(time, dsp(8,:))
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
——————————————————————————————————————————————
用到的函数 wilson
——————————————————————————————————————————————
function =TransResp4(kk,cc,mm,fd,bcdof,nt,dt,q0,dq0)
%--------------------------------------------------------------------------
%Purpose:
% The function subroutine TransResp4.m calculates transient response of
% a structural system using Wilson integration scheme.
%Synopsis:
% =TransResp4(kk,cc,mm,fd,bcdof,nt,dt,q0,dq0)
%Variable Description:
% Input parameters
% kk, cc, mm - stiffness, damping and mass matrices
% fd - Input or forcing influence matrix
% bcdof - Boundary condition dofs vector
% nt - Number of time steps
% dt - Time step size
% q0, dq0 - Initial condition vectors
% Output parameters
% acc - Acceleration response
% vel - Velocity response
% dsp - Displacement response
%--------------------------------------------------------------------------
%(1) initial condition
%--------------------------------------------------------------------------
=size(kk);
dsp=zeros(sdof,nt); % displacement matrix
vel=zeros(sdof,nt); % velocity matrix
acc=zeros(sdof,nt); % acceleration matrix
dsp(:,1)=q0; % initial displacement
vel(:,1)=dq0; % initial velocity
theta=1.4; % select the parameters
%--------------------------------------------------------------------------
%(2) Wilson integration scheme for time integration
%--------------------------------------------------------------------------
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1));
% compute the initial acceleration (t=0)
ekk=kk+mm*(6/(theta*dt)^2)+cc*(3/(theta*dt));
% compute the effective stiffness matrix
for i=1:sdof % assign zero to dsp, vel, acc of the dofs associated with bc
if bcdof(i)==1
dsp(i,1)=0;
vel(i,1)=0;
acc(i,1)=0;
end
end
for it=1:nt % loop for each time step
cfm=dsp(:,it)*(6/(theta*dt)^2)+vel(:,it)*(6/(theta*dt))+2*acc(:,it);
cfc=dsp(:,it)*(3/(theta*dt))+2*vel(:,it)+acc(:,it)*(theta*dt/2);
efd=fd(:,it)+theta*(fd(:,it+1)-fd(:,it))+mm*cfm+cc*cfc;
%compute the effective force vector
dtheta=inv(ekk)*efd; % find the displacement at time t+ dt
acc(:,it+1)=(dtheta-dsp(:,it))*(6/(theta^3*dt^2))...
-vel(:,it)*(6/(theta^2*dt))+acc(:,it)*(1-3/theta);
% find the acceleration at time t+dt
vel(:,it+1)=vel(:,it)+acc(:,it+1)*dt/2+acc(:,it)*dt/2;
% find the velocity at time t+dt
dsp(:,it+1)=dsp(:,it)+vel(:,it)*dt...
+(acc(:,it+1)+2*acc(:,it))*(dt^2/6);
% find the displacement at time t+dt
for i=1:sdof % assign zero to acc, vel, dsp of the dofs associated with bc
if bcdof(i)==1
dsp(i,it+1)=0;
vel(i,it+1)=0;
acc(i,it+1)=0;
end
end
end
if cc(1,1)==0
disp('The transient response results of undamping system')
else
disp('The transient response results of damping system')
end
disp('The method is Wilson integration')
%--------------------------------------------------------------------------
% The end
%--------------------------------------------------------------------------
运行不了，求指导。 这个错误 程序太长没有看，说说报错部分的意思吧。plot函数绘图的时候，x方向上的间隔必须等距，简单说1,2,3,4,5可以，但是1,2,3,4,8就不行，因为4到8之间的跨度和前面不等，这个时候需要把4到8之间的数据补齐了才能用plot函数。 许伍sky 发表于 2013-1-25 13:50 static/image/common/back.gif
这个错误
没试过
time=0:dt:nt*dt; => time的size为1*1501
dsp=zeros(sdof,nt) => dsp的size为1*1500
好东西，大家一起分享~ ChaChing 发表于 2013-1-29 22:39 static/image/common/back.gif
没试过
time=0:dt:nt*dt; => time的size为1*1501
dsp=zeros(sdof,nt) => dsp的size为1*1500
那这个怎么改呢，忙烦您看一下，不甚感激 那这个怎么修改啊，高手再哪里？
页:
[1]