菜鸟编辑的Newmark法单自由度线性体系，请大侠们指正批评

xzning123

此程序参考的是《结构动力学：理论及其在地震工程中的应用》一书的Newmark法
%%%%Newmark method: linear systems,single degree freedom
clc
clear all
close all
%%%%intial values of the caculated system
tic;
m=0.2533;
c=0.1592;
k=10;
%%%%intial caculation
t0=0.6;     %the end of the excitation
total=1.0;   %totla caculation time
dt=0.1;     %time step
t1=0:dt:total;
s=(t1>0.6);
t1(s)=t1(s)*0;
t1(~s)=t1(~s);
pt=10*sin(pi*t1/0.6);
n=total/dt;
n1=t0/dt;
u1=0;
v1=0;
gama=0.5;
beta=1/6;
a1=(pt(1)-c*v1-k*u1)/m;
keq=k+gama*c/beta/dt+m/beta/dt^2;
x=m/beta/dt+gama*c/beta;
y=m/2/beta+c*dt*(gama/2/beta-1);
x0=gama/beta/dt;
x1=gama/beta;
x2=dt*(1-gama/2/beta);
x3=1/beta/dt^2;
x4=1/beta/dt;
x5=1/2/beta;
%%%%form the vector of the%%%%
u(1)=u1;
v(1)=v1;
a(1)=a1;
%%%%for every time step%%%%
for i=2:1:n+1
    dp(i-1)=pt(i)-pt(i-1);
    dpeq(i-1)=dp(i-1)+x*v(i-1)+y*a(i-1);
    du(i-1)=dpeq(i-1)/keq;
    dv(i-1)=x0*du(i-1)-x1*v(i-1)+x2*a(i-1);
    da(i-1)=x3*du(i-1)-x4*v(i-1)-x5*a(i-1);
    u(i)=u(i-1)+du(i-1);
    v(i)=v(i-1)+dv(i-1);
    a(i)=a(i-1)+da(i-1);
end
h=toc;
fprintf('the caculate time is %10.6f \n',h)
%%%%to figure out the caculation result
t=0:dt:total;
plot(t',u(:,1:n+1))   %the value corresponding to u(1) is not real