Newmark单自由度，自己写的，大家批评指正

xzning123

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