基于 MATLAB 实现对结构动力响应的几种算法的验证1. 算例首先，本文给出一算例， 结构在外力谐振荷载 P (t ) = P 0 sin θt 作用下，分别利用理论解法，杜哈梅积分， Wilson-θ 法求出该结构的位移时程反应。
其中：m = 3.5×103 kg ， P 0 = 1.0×104 N ， k =1.3584515×107 ，ξ=0.05 ，θ=52.3s −1 ，ω=62.3s −1 ，⋅D ω= ω 1-ξ2=62.222 ，初始位移、速度v (0) = 0 ，v (0) = 0 ；2. 算法验证2.1 理论解法运动方程为： mv+cv+kv=0P sin t ϑ由线性代数解出其理论解为:]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(t]sin )0()0(cos )0([2222222202222222222t t m P t m P ev v t v e v D DD tD DD t θξωθθθωθωξθωωωθωωξωξωθωξθωθωωξωωξωξω-++-+⋅++++=--由于初始位移v(0) =0 ，v(0) =0 ；则：]cos 2)sin -[(]4)-[(t]sin )-(2cos [2]4)-[(22222222022222222220t t m P t m P e v D DD tθξωθθθωθωξθωωωθωωξωξωθωξθωθξω-++-+⋅+=-v(t ) =-3.115te⋅ 1.05269898×410-⋅[6.230cos62.222t −18.106sin 62.222t]+2.012808757 610-⨯⋅[1146sin 52.3t −325.829cos52.3t]可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。
2.2Wilson-θ法Wilson-θ法是Wilson 于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。
该方法假定在θt ∆时程步长内，体系的加速度反应按线性变化。
对于地震持续时间内的每一个微小时 段t ∆ ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。
下面以第i+1时段（时刻时刻到1+i i t t ）为例：{}{}[]{}[][][][]{}[]{}[]{}{}[]{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}[][]{}[][]{}附录。
移时程图，详细程序见进行编程分析，画图位可以用；时刻的相对加速度反应，求得第代入结构运动微分方程和最后将和速度增量移反应时刻的各质点的相对位得到第，和速度增量，分别加上位移增量时刻的地震反应为基础再以各质点在第；和速度增量内位移增量正常时段再求解第；内的加速度增量正常时段计算出第先利用位移增量的加长时段内各质点的计算MATLAB x xC xx x x x x x xx x x xx x x t xt x t x t x xxt x x x i xxx xx x xx C x x xP C K P K x x i i g i i i i i i i i i i i i i i i i i i i i i g iii )(-t )5(t )4(631)(2t 1)3()2(6t 1i )2()23()36(136;t )1(1i 11i 11i 1i 11i 1i 1i 11i 1i 1122122211,+-+-+++++++++++++***-*K M +M +=∆+=∆+=∆∆∆+∆+∆=∆+∆=∆∆∆∆+--∆=∆∆∆+∆+++M +∆M -=∆+M +K =∆=∆∆∆=+ τττττθτττττθτττττττ2.3 杜哈梅积分杜哈梅积分在考虑阻尼的情况是：ττωτωωωξωωτξωξωd t p e m v v t v e v D tt DD DD t )(sin )(1t]sin )0()0(cos )0([0)(-+++=⎰---可以用 MA TLAB 进行编程分析，画图位移时程图，详细程序见附录。
3. 位移时程反应对比分析利用MATLAB将理论解法，杜哈梅积分，Wilson-θ法求解出来的位移时程反应画在同一张图中，进行比较分析。
从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。
4.结论本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用Wilson-θ法）进行相互验证，从最后的位移分析图对比上，可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。
附录：MA TLAB 源程序%理论解，杜哈梅积分，Wilson-θ法程序clc;clearh1=figure(8);set( h1, 'color','w')%理论解法t=0:0.01:1;v=1.05269898*10^(-4)*exp(-3.115*t).*(6.230*cos (62.222*t)-18.106*sin(62.222*t))+2.012808757*1 0^(-6)*(1146*sin(52.3*t)-325.829*cos(52.3*t)); plot(t,v,'k')hold on%杜哈梅积分法aa=1;%输入时间长度bb=0.01;%输入精度t=bb:bb:aa;t1=t;theta=52.3;%输入荷载频率w=62.3;%输入自振频率m=3500;%输入质量p0=10000;%输入荷载幅值p0=p0*ones(1,aa/bb);p=p0.*sin(theta*t);%荷载函数for i=1:(aa/bb)for j=1:icanshu1(j)=p(j)/(m*w)*bb*sin(w*(t(i)-t1(j)));%杜哈梅积分中的被积函数endy(i)=sum(canshu1);%%位移值endfor i=1:aa/bb-1v1(i)=(y(i+1)-y(i))/bb;%计算速度endfor i=1:(aa/bb-2)a(i)=(v1(i+1)-v1(i))/bb;%计算加速度endhold onplot(t,y,'r')%画位移图hold on%Wilson-θ法dt=0.01;ct=1.4;ndzh=100;k=13584515;c=21805;t=0:dt:ndzh*dt;ag=10000*sin(52.3*t);ag1=ag(1:ndzh);ag2=ag(2:ndzh+1);agtao=ct*(ag2-ag1);wyi1=0;sdu1=0;jsdu1=0;wyimt=0;sdumt=0;jsdumt=0;for i=1:ndzhkxin=k+(3/(ct*dt))*c+(6/(ct*ct*dt*dt))*m;%kxin为新的刚度dpxin=-m*agtao(i)+m*(6/(ct*dt)*sdu1+3*jsdu1)+c*(3*sdu1+ct*dt/2*jsdu1); %新的力增量dxtao=kxin\dpxin;dtjsdu=6*dxtao/(ct*(ct^2*dt^2))-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1;jsdu=jsdu1+dtjsdu;dtsdu=(dt/2)*(jsdu+jsdu1);sdu=sdu1+dtsdu;dtwyi=dt*sdu1+(1/3)*dt^2*jsdu1+(dt^2/6)*jsdu;wyi=wyi1+dtwyi;jsdu=-m\(m*ag2(i)+c*sdu+k*wyi); % 调整加速度wyi1=wyi;sdu1=sdu;jsdu1=jsdu;wyimt=[wyimt wyi];sdumt=[sdumt sdu];jsdumt=[jsdumt jsdu];endplot(t,-wyimt/3000,'b');hold offlegend('\fontsize{9}\fontname{ 黑体} 理论解','\fontsize{9}\fontname{黑体} 杜哈梅积分法','\fontsize{9}\fontname{黑体} wilson-{\theta}法')%基于双线性滞回模型的单自由度体系的地震能量分析程序%质量57041kg，阻尼36612 N·s/m,初始刚度2350000N/m,刚度折减系数0.2，屈服位移0.01m，采用ELCENTRO波%参数替换直接在下面修改，然后运行clcformat long;m=57041;%质量ug=importdata('ELCENTRO.txt');%地震波txt 文件ug=ug/100;P=-m*ug;num=size(P,1);c=36612;%阻尼k1=2350000;%初始刚度k2=k1*0.2;a=zeros(1,num);v=zeros(1,num);x=zeros(1,num);%加速度速度位移Ei=zeros(1,num);Ek=zeros(1,num);Ed=zeros(1,num);Eh=zeros(1,num);%输入能动能阻尼耗能滞回耗能EI=zeros(1,num);EK=zeros(1,num);ED=zeros(1,num);EH=zeros(1,num);%累积的各种能量time=zeros(1,num);a(1)=P(1)/m;hfl=zeros(1,num);h=0.02;%地震波采样间隔Xy=0.01;%屈服位移pxmax=0;nxmax=0;pd=0;%双线型滞回模型折线的标识，0表示弹性，1表示正向弹塑性，2表示反向弹性，-1表示反向弹塑性，-2表示正向弹性for ii=1:numif pd==0 %弹性阶段k=k1;if x(ii)>Xypd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+k1*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+dtf;a(ii)=(P(ii)-hfl(ii+1)-c*v(ii))/m;elseif x(ii)<-Xypd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)+Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp-k1*Xy;kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=-Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=-k1*Xy+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendi f pd==1 %正向弹塑性k=k2;if v(ii)<0pd=2;b=[(a(ii)-a(ii-1))/2/ha(ii-1) v(ii-1)];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth^2/2+(a(ii)-a(ii-1))*dth^3/h/6;ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*0+k1*Xy+k2*(xp-Xy);kd=k1+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(ii)=xp+dtx;v(ii)=0+dtv;dtf=k1*dtx;hfl(ii)=k1*Xy+k2*(xp-Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;pxmax=xp;endendif pd==2 %反向弹性k=k1;if x(ii)>pxmaxpd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-pxmax];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+k1*Xy+k2*(pxmax-Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=pxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+k2*(pxmax-Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;elseif x(ii)<(pxmax-2*Xy)pd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-pxmax+2*Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+(pxmax*k2-k1*Xy-k2*Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=pxmax-2*Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=pxmax*k2-k1*Xy-k2*Xy+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendif pd==-1 %反向弹塑性k=k2;if v(ii)>0pd=-2;b=[(a(ii)-a(ii-1))/2/h a(ii-1)v(ii-1)];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth^2/2+(a(ii)-a(ii-1))*dth^3/h/6;ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*0-k1*Xy+k2*(xp+Xy);kd=k1+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*0/hp+3*ap)+c*(3*0+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*0-hp*ap/2;x(ii)=xp+dtx;v(ii)=0+dtv;dtf=k1*dtx;hfl(ii)=-k1*Xy+k2*(xp+Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;nxmax=x(ii);endendif pd==-2 %正向弹性k=k1;if x(ii)<nxmaxpd=-1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-nxmax];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp-k1*Xy+k2*(nxmax+Xy);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=nxmax+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=-k1*Xy+k2*(nxmax+Xy)+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;elseif x(ii)>(nxmax+2*Xy)pd=1;b=[(a(ii)-a(ii-1))/6/h a(ii-1)/2 v(ii-1)x(ii-1)-nxmax-2*Xy];% 拐点处理d=roots(b);for j=1:length(d)if isreal(d(j))==1dth=d(j);endendhp=h-dth;vp=v(ii-1)+a(ii-1)*dth+((a(ii)-a(ii-1))*dth^2/h/2);ap=a(ii-1)+(a(ii)-a(ii-1))*dth/h;pp=m*ap+c*vp+(k1*Xy+nxmax*k2+Xy*k2);kd=k2+3*c/hp+6*m/hp/hp;dtpd=P(ii+1)-pp+m*(6*vp/hp+3*ap)+c*(3*vp+hp*ap/2);dtx=dtpd/kd;dtv=3*dtx/hp-3*vp-hp*ap/2;x(ii)=nxmax+2*Xy+dtx;v(ii)=vp+dtv;dtf=k2*dtx;hfl(ii)=k1*Xy+nxmax*k2+Xy*k2+dtf;a(ii)=(P(ii)-hfl(ii)-c*v(ii))/m;endendi f ii>=numbreak;endkd=k+3*c/h+6*m/h/h;dtpd=P(ii+1)-P(ii)+m*(6*v(ii)/h+3*a(ii))+c*(3*v(ii)+h*a(ii)/2);dtx=dtpd/kd;dtv=3*dtx/h-3*v(ii)-h*a(ii)/2;x(ii+1)=x(ii)+dtx;v(ii+1)=v(ii)+dtv;dtf=k*dtx;hfl(ii+1)=hfl(ii)+dtf;a(ii+1)=(P(ii+1)-hfl(ii+1)-c*v(ii+1))/m;Ei(ii+1)=-m*ug(ii+1)*v(ii+1)*h;Ed(ii+1)=c*v(ii+1)^2*h;Ek(ii+1)=m*a(ii+1)*v(ii+1)*h;Eh(ii+1)=hfl(ii+1)*v(ii+1)*h;for jj=1:num-1EI(jj+1)=EI(jj)+Ei(jj+1);ED(jj+1)=ED(jj)+Ed(jj+1);EK(jj+1)=EK(jj)+Ek(jj+1);EH(jj+1)=EH(jj)+Eh(jj+1);endendfor j=1:numtime(j)=h*(j-1);endfid=fopen('EL波弹塑性各种能.txt','wt'); %输出结果到EL波弹塑性各种能.txt fprintf(fid,'%g %g %g %g%g\n',[time;EI;ED;EH;EK]);fclose(fid);plot(time,EI);hold onplot(time,ED);hold onplot(time,EK);hold onplot(time,EH);。