动力学系统时域响应计算的六种方法Matlab源程序(Newmark,Houbolt法,中心差分法)

Newmark法Matlab源程序


function [acc,vel,dsp]=Newmark_2(kk,cc,mm,fd,nt,dt,q0,dq0)
%输入参数
% kk------刚度矩阵
% mm------质量矩阵
% cc------阻尼矩阵
% q0------初始位移
% dq0------初始速度
% dt------时间步长
% nt------总的计算步数，等于结束时间除以dt
%返回值
% dsp------位移
% vel------速度
% acc------加速度
[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
alpha=0.5; beta=0.5; % select the parameters
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1)); % compute the initial acceleration (t=0)
ekk=kk+mm/(alpha*dt^2)+cc*beta/(alpha*dt);
for it=1:nt % loop for each time step
cfm=dsp(:,it)/(alpha*dt^2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);
cfc=dsp(:,it)*beta/(alpha*dt)+vel(:,it)*(beta/alpha-1)...
+acc(:,it)*(0.5*beta/alpha-1)*dt;
efd=fd(:,it)+mm*cfm+cc*cfc; % compute the effective force vector
dsp(:,it+1)=inv(ekk)*efd; % find the displacement at time t+dt
acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...
-acc(:,it)*(0.5/alpha-1); % find the acceleration at time t+dt
vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt;% find the velocity at time t+dt
end


中心差分法Matlab源程序



function [acc,vel,dsp]=central_diff(kk,cc,mm,fd,nt,dt,q0,dq0)
%输入参数
% kk------刚度矩阵
% mm------质量矩阵
% cc------阻尼矩阵
% q0------初始位移
% dq0------初始速度
% dt------时间步长
% nt------总的计算步数，等于结束时间除以dt
%返回值
% dsp------位移
% vel------速度
% acc------加速度
[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
acc(:,1)=inv(mm)*(fd(:,1)-kk*dsp(:,1)-cc*vel(:,1));
dsp0=dsp(:,1)-vel(:,1)*dt+0.5*acc(:,1)*dt^2;
ekk=mm/dt^2+cc/(2*dt);
efd=fd(:,1)-(kk-2*mm/dt^2)*dsp(:,1)-(mm/dt^2-cc/(2*dt))*dsp0;
dsp(:,1+1)=inv(ekk)*efd;
for it=2:nt % loop for each time step after first step
efd=fd(:,it)-(kk-2*mm/dt^2)*dsp(:,it)-