function [d,v,a]=Newmark(K,M,C,f,d1,v1,a1,dt,tend)
%K---剛度
%M---質(zhì)量
%C---阻尼
%d1---初始位移
%v1---初始速度
%dt---時(shí)間步長
%tend---結(jié)束時(shí)間
%%%返回值
%d----位移
%v----速度
%a----加速度
gama=0.5;
beta=1/6;
alpha0=1/beta/dt^2;
alpha1=gama/beta/dt;
alpha2=1/beta/dt;
alpha3=1/2/beta-1;
alpha4=gama/beta-1;
alpha5=dt/2*(gama/beta-2);
alpha6=dt*(1-gama);
alpha7=gama*dt;
K1=K+alpha0*M+alpha1*C;
d=zeros(1,floor(tend/dt)+1);
v=zeros(1,floor(tend/dt)+1);
a=zeros(1,floor(tend/dt)+1);
d(1,1)=d1;
v(1,1)=v1;
a(1,1)=a1;%M\(f(1,1)-K*d1-C*v1);
for i=2:1:floor(tend/dt)+1
    t=(i-1)*dt;
    f1=f(:,i)+M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+alpha3*a(:,i-1))...
        +C*(alpha1*d(:,i-1)+alpha4*v(:,i-1)+alpha5*a(:,i-1));
    d(:,i)=K1\f1;
    v(:,i)=alpha1*(d(:,i)-d(:,i-1))-alpha4*v(:,i-1)-alpha5*a(:,i-1);
    a(:,i)=M\(f(:,i)-C*v(:,i)-K*d(:,i));
end
return