纽马克地震时程分析

西风独自凉

function newmark1(t_max,dt)
%t_max为持续时间，dt为时间步长
%t_max=10
%dt=0.01
M=[1];
C=[0.2*pi];
K=[4*pi*pi];
u=[0.2]';
v=[0]';
a=[0]';
t(1)=0;               %时间
x(:,1)=u;             %位移
x1(:,1)=v;            %速度
x2(:,1)=a;            %加速度
%newmark参数
gama=0.5;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt;
%等效刚度矩阵
K1=K+a0*M+a1*C;
i=1;
Xg=load('C:\Documents and Settings\Administrator.WWW-ADAB10B430E\桌面\acc_ElCentro_0.34g_0.02s.txt');
for t=dt:dt:t_max;
      q(:,i+1)=-Xg(i)+M*(a0*x(:,i)+a2*x1(:,i)+a3*x2(:,i))+C*(a1*x(:,i)+a4*x1(:,i)+a5*x2(:,i));
      x(:,i+1)=inv(K1)*q(:,i+1);
      x2(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*x1(:,i)-a3*x2(:,i);
      x1(:,i+1)=a1*(x(:,i+1)-x(:,i))-a4*x1(:,i)-a5*x2(:,i);
      i=i+1;
end  
x1=x(1:t_max/dt)';
t=dt:dt:t_max;
plot(t,x1,'-')
xlabel('t[sec]')
ylabel('x(m)')



输入：newmark1(10,0.01)

电脑克星LK

（1）这是计算单自由度体系的反应，最好在标题中说明一下，因为新手可能不清楚。而多自由度体系的计算比这个复杂得多。
（2）给的加速度文件是一个正弦波，而不是任何一种El Centro地震波，与标题“地震时程分析”不太相符。地震波一般给的是加速度，这个正弦波是不是位移波？位移波和加速度波计算是有区别的。
（3） x(:,i+1)=inv(K1)*q(:,i+1); 这个计算式中，把inv（K1)的计算放到循环外面比较好。如果K1很大的话，每次循环都对它求逆还是比较费时间的，当然，在这个例子中影响几乎为零。

hongtianwenduo

正在学习，谢谢!

小星星

正在学习，谢谢!

无奈和孤独的人

向高手致敬，学习了。

scan11

致敬，学习了