cendiff.m

动力学教程中的数值方法

function [X,Xd,Xdd,t] = cendiff(M,C,K,f,dt,Xi,Xdi) 
%================================================================ 
%function [X,t] = newmark(M,C,K,f,dt,gamma,beta,Xi,Xdi) 
% 
%  (c) Ajax, (all rights reserved) 
%  Central South University 
%  October 13, 2006 
% 
%  This function is intended to perform the numerical 
%  integration of a structural system subjected to an 
%  external dynamic excitation such as a wind or earthquake. 
%  The structural model is assumed to be a lumped mass shear 
%  model.The integration scheme utilized for this analysis 
%  is the central difference method. The central difference 
%  method is an implicit time steping scheme so stability of 
%  the system need not be considered.
% 
%  Input Variables: 
%  [M] = Mass Matrix (nxn) 
%  [C] = Damping Matrix (nxn) 
%  [K] = Stiffness Matrix (nxn) 
%  {f} = Excitation Vector (mx1) 
%  dt = Time Stepping Increment 
%  Xi = Initial Displacement Vector (nx1) 
%  Xdi = Initial Velocity Vector (nx1) 
% 
%  Output Variables: 
%  {t} = Time Vector (mx1) 
%  [X] = Response Matrix (mxn) 
%================================================================ 

  n = size(M,1); 
  fdimc = size(f,2); 
  fdimr = size(f,1); 
  X(:,2)=Xi;
  Xd(:,1)=Xdi;
  
%  Check Input Excitation 
%  ====================== 
   if(fdimc==n)
       f=f';
   end
   m=size(f,2);
   
%  Coefficients 
%  ============ 
        c0 = 1/(dt*dt) ; 
        c1 = 1/(2*dt) ; 
        c2 = 2*c0 ; 
        
%Initialize Mass Matricies 
% ============================== 
        Meff = c0*M + c1*C ; 
        Minv = inv(Meff) ;
        
% Initial Acceleration
% ==================== 
        R0=f(:,1);
        Xdd(:,1)= inv(M)*(f - K*Xi - C*Xdi);

% Perform First Step 
% ==================
        X(:,1) = X(:,2) - 1/(2*c1)*Xd(:,1) + 1/(2*c0)*Xdd(:,1) ; 
        f(:,1) = f(:,1) - (K - c2*M)*X(:,2) - (c0*M - c1*C)*X(:,1);
        X(:,3) = Minv*f(:,1);
        Xd(:,2) = c1*(X(:,3) - X(:,1));
        Xdd(:,2) = c0*(X(:,3) - 2*X(:,2) + X(:,1));
        
%Perform Subsequent Steps 
% ======================== 
        for i=1:size(f,2)-1; 
         X(:,i+1) = X(:,i+2) - 1/(2*c1)*Xd(:,i+1) + 1/(2*c0)*Xdd(i+1) ; 
         f(:,i+1) = f(:,i+1) - (K - c2*M)*X(:,i+2) - (c0*M - c1*C)*X(:,i+1);
         X(:,i+3) = Minv*f(:,i+1);
         Xd(:,i+2) = c1*(X(:,i+3) - X(:,i+1));
         Xdd(:,i+2) = c0*(X(:,i+3) - 2*X(:,i+2) + X(:,i+1));
        end; 

         
%  Strip Off Padded Response Zeros 
%  =============================== 
  X = X'
  Xd=Xd'
  Xdd=Xdd'
  
%  Generate the Time Vector 
%  ======================== 
  for i=0:1:size(X,1)-1 
    t(i+1,1) = i*dt; 
  end; 
  t