matlab 单自由度 杜哈梅积分

结构动力学单自由度体系matlab

在t ∈[0,π/w]范围内的位移时程图在Duhamel 积分、分段线性插值法、中心差分法下的比较。对表达式中的结构属性自己取值，取值结果如程序中所示。Matlab 程序如下：

clear;clf;clc;

syms t ;

m=1;wn=10;k=100;w=5;p0=100;tm=pi/250;u0=0;u0m=0;

u1=p0*(-sin(w*t)*wn+sin(wn*t)*w)/m/wn/(w^2-wn^2);

uf=subs(u1,t,0:pi/250:pi/w);

t=0:pi/250:pi/w;

grid on ;hold on ;

plot(t,uf,'r-.d' , 'MarkerSize' ,10)

A=(sin(wn*tm)-wn*tm*cos(wn*tm))/k/wn/tm;

B=(wn*tm-sin(wn*tm))/k/wn/tm;

C=cos(wn*tm);

D=sin(wn*tm)/wn;

Am=(wn*tm*sin(wn*tm)+cos(wn*tm)-1)/k/tm;

Bm=(1-cos(wn*tm))/k/tm;

Cm=-1*wn*sin(wn*tm);

Dm=cos(wn*tm);

i=0:1:50;

u=[];um=[];p=[];p(i+1)=p0*sin(w*t);

z=1:1:50;a=1;

for v=z

um(1)=0;

u(1)=0;

um(a+1)=Am*p(v)+Bm*p(v+1)+Cm*u(v)+Dm*um(v);

u(a+1)=A*p(v)+B*p(v+1)+C*u(v)+D*um(v);

a=a+1;

end

plot(t,u,'k--x' , 'MarkerSize' ,15)

u0mm=(p0-k*u0)/m;

u_1=u0-tm*u0m+tm^2*u0mm/2;

k_=m/tm^2;

a=m/tm^2;

b=k-2*m/tm^2;

p_=[];u3f1=[];

z=1:1:50;

a_=1;

for v=z

u3f1(1)=u_1;

u3f1(2)=u0;

p_(a_)=p(v)-a*u3f1(v)-b*u3f1(v+1);

u3f1(a_+2)=p_(v)/k_;

a_=a_+1;

end

l=2:1:52;

u3f=[u3f1(l)];

plot(t,u3f,'g-v' , 'MarkerSize' ,8)

legend('ÀíÂÛÖµ', '·Ö¶ÎÏßÐÔ²åÖµ·¨', 'ÖзֲîÖµ·¨', 'location' , 'northeast' )

地震波自己找

结构动力学单自由度体系matlab

在t ∈[0,π/w]范围内的位移时程图在Duhamel 积分、分段线性插值法、中心差分法下的比较。对表达式中的结构属性自己取值，取值结果如程序中所示。Matlab 程序如下：

clear;clf;clc;

syms t ;

m=1;wn=10;k=100;w=5;p0=100;tm=pi/250;u0=0;u0m=0;

u1=p0*(-sin(w*t)*wn+sin(wn*t)*w)/m/wn/(w^2-wn^2);

uf=subs(u1,t,0:pi/250:pi/w);

t=0:pi/250:pi/w;

grid on ;hold on ;

plot(t,uf,'r-.d' , 'MarkerSize' ,10)

A=(sin(wn*tm)-wn*tm*cos(wn*tm))/k/wn/tm;

B=(wn*tm-sin(wn*tm))/k/wn/tm;

C=cos(wn*tm);

D=sin(wn*tm)/wn;

Am=(wn*tm*sin(wn*tm)+cos(wn*tm)-1)/k/tm;

Bm=(1-cos(wn*tm))/k/tm;

Cm=-1*wn*sin(wn*tm);

Dm=cos(wn*tm);

i=0:1:50;

u=[];um=[];p=[];p(i+1)=p0*sin(w*t);

z=1:1:50;a=1;

for v=z

um(1)=0;

u(1)=0;

um(a+1)=Am*p(v)+Bm*p(v+1)+Cm*u(v)+Dm*um(v);

u(a+1)=A*p(v)+B*p(v+1)+C*u(v)+D*um(v);

a=a+1;

end

plot(t,u,'k--x' , 'MarkerSize' ,15)

u0mm=(p0-k*u0)/m;

u_1=u0-tm*u0m+tm^2*u0mm/2;

k_=m/tm^2;

a=m/tm^2;

b=k-2*m/tm^2;

p_=[];u3f1=[];

z=1:1:50;

a_=1;

for v=z

u3f1(1)=u_1;

u3f1(2)=u0;

p_(a_)=p(v)-a*u3f1(v)-b*u3f1(v+1);

u3f1(a_+2)=p_(v)/k_;

a_=a_+1;

end

l=2:1:52;

u3f=[u3f1(l)];

plot(t,u3f,'g-v' , 'MarkerSize' ,8)

legend('ÀíÂÛÖµ', '·Ö¶ÎÏßÐÔ²åÖµ·¨', 'ÖзֲîÖµ·¨', 'location' , 'northeast' )

地震波自己找