动力学系统时域响应计算的六种方法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)-