【求助，matlab动力学编程出错】

许伍sky

%--------------------------------------------------------
%  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=[0 0;4 0;4 3;8 0;8 3;12 0];            %node_coordinate
en=[1 2;1 3;2 3;2 4;3 4;3 5;4 5;4 6;5 6];   % 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*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; % 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 size  with  k
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=[cos(theta(loopi)) sin(theta(loopi)); -sin(theta(loopi)) cos(theta(loopi))];
    t=[lmd zeros(2); zeros(2) lmd];
    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
[acc,vel,dsp]=TransResp4(k,c,m,fd,bcdof,nt,dt,q0,dq0);
plot(time, dsp(8,:))
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
——————————————————————————————————————————————
用到的函数 wilson
——————————————————————————————————————————————


function [acc,vel,dsp]=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:
%     [acc, vel, dsp]=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
%--------------------------------------------------------------------------
[sdof,n2]=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
%--------------------------------------------------------------------------

许伍sky

运行不了，求指导。

许伍sky

这个错误

alexqxp

程序太长没有看，说说报错部分的意思吧。plot函数绘图的时候，x方向上的间隔必须等距，简单说1,2,3,4,5可以，但是1,2,3,4,8就不行，因为4到8之间的跨度和前面不等，这个时候需要把4到8之间的数据补齐了才能用plot函数。

ChaChing

许伍sky 发表于 2013-1-25 13:50
这个错误

没试过
time=0:dt:nt*dt; => time的size为1*1501
dsp=zeros(sdof,nt) => dsp的size为1*1500

laocl07

好东西，大家一起分享~

许伍sky

ChaChing 发表于 2013-1-29 22:39
没试过
time=0:dt:nt*dt; => time的size为1*1501
dsp=zeros(sdof,nt) => dsp的size为1*1500

那这个怎么改呢，忙烦您看一下，不甚感激

许伍sky

那这个怎么修改啊，高手再哪里？