function [uu,vv,aa,ttt]=NewmarkBeta1(t,dt,delta,beta,M,C,K,acc_x,Constr)
n=length(t);
beta=0.25;
gamma=0.5;
a0=1/(beta*dt^2);
a1=gamma/(beta*dt);
a2=1/(beta*dt);
a3=1/(2*beta)-1;
a4=gamma/beta-1;
a5=(dt/2)*(gamma/beta-2);
a6=dt*(1-gamma);
a7=dt*gamma;
uu(:,1)=zeros(size(M,1),1);;%位移初始值
vv(:,1)=zeros(size(M,1),1);%速度初始值
aa(:,1)=zeros(size(M,1),1);%加速度初始值
P(:,1)=zeros(size(M,1),1); %荷载初始值
P1(:,1)=zeros(size(M,1),1);%等效荷载初始值
% I=ones(size(M,1),1);; %6*1单位阵
I=zeros(size(M,1),1);; %6*1单位列向量
for i=1:size(K,1)/3
  
    I(3*(i-1)+1)=1;%确定惯性力施加的自由度
end
K1=K+a0*M+a1*C; %形成等效刚度矩阵
for i=1:n-1

    P(:,i+1)=-M*I*acc_x(i+1);%t(i+1)时刻结构的荷载
    P1(:,i+1)=P(:,i+1)+M*(a0*uu(:,i)+a2*vv(:,i)+a3*aa(:,i))+C*(a1*uu(:,i)+a4*vv(:,i)+a5*aa(:,i));%t(i+1)时刻结构的等效荷载
%     % %施加约束--乘大数法;当代入结构整体的刚度矩阵、质量矩阵、阻尼矩阵需要施加边界约束
%     for iConstr=1:size(Constr,1)
%            for j=2:4
%                   if ~isnan(Constr(iConstr,j))
%                      K1(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4)=...
%                      1e12*K1(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4);
%                      P1(3*Constr(iConstr,1)+j-4)=Constr(iConstr,j)*...
%                      K1(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1)+j-4);
%                   end
%            end
%     end  
    uu(:,i+1)=(K1)\P1(:,i+1);%t(i+1)时刻结构的位移
    aa(:,i+1)=a0*(uu(:,i+1)-uu(:,i))-a2*vv(:,i)-a3*aa(:,i);%t(i+1)时刻结构的加速度
    vv(:,i+1)=vv(:,i)+a6*aa(:,i)+a7*aa(:,i+1);%t(i+1)时刻结构的速度
    ttt(i+1)=dt*(i-1);
end