基于 MATLAB 实现对结构动力响应的几种算法的验证
1. 算例
首先，本文给出一算例， 结构在外力谐振荷载 P(t) = P0 sinθt 作用下，分别利用理论解法，杜哈梅积分， Wilson-θ 法求出该结构的位移时程反应。其中：
m = 3.5×103 kg ， P0 = 1.0×104 N ， k =1.3584515×107 ，ξ=0.05 ，θ=52.3s?1 ，ω=62.3s?1 ，
?D= ω 1-ξ =62.222 ，初始位移、速度v(0) = 0 ，v(0) = 0 ；
2
?
2. 算法验证
2.1 理论解法
运动方程为： mv+cv+kv=P0sin?t由线性代数解出其理论解为:
v?e???t[v(0)cos?Dt??e???t2v(0)???v(0)?D222sin?Dt]2?2?2?(?2-?2)sin?Dt]
?P0m[(?-?)?4???]22?[2??cos?Dt??DP022?[(?-?)sin?t?2???cos?t]222222m[(?-?)?4???]
由于初始位移v(0) =0 ，v(0) =0 ；则：
v?e????t?P0m[(?2-?2)2?4?2?2?2]?[2??cos?Dt?2?2?2?(?2-?2)?Dsin?Dt]
P0[(?2-?2)sin?t?2???cos?t]222222m[(?-?)?4???]-3.115t v(t ) =e?1.05269898×10?4?[6.230cos62.222t ?18.106sin 62.222t]
?6+2.012808757 ?10?[1146sin 52.3t ?325.829cos52.3t]
可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。
2.2Wilson-?法
Wilson-?法是Wilson于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。该方法假定在??t时程步长内，体系的加速度反应按线性变化。对于地震持续时间内
的每一个微小时 段?t ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。下面以第i+1时段（ti时刻到ti?1时刻）为例：
(1)计算????t的加长时段内各质点的位移增量??x??;??x????K???P?i2?????1?3?K???????6??????C??i??P????g?i????(??????1???x6????i?3????i)??C?i(3?x??i?????i)?xxx2??i；???(2)先利用??x??计算出第i?1正常时段?t内的加速度增量x??i????x6??2??i?(??x?????x?22??i)??x??(3)再求解第i?1正常时段?t内位移增量??x?和速度增量??x；????t??x??2(??xi?1??i)??x12?t2??i??t????i???i?1???x??t?xxx36????x?和速度增量??x(4)再以各质点在第ti时刻的地震反应为基础，分别加上位移增量，??i?1得到第ti?1时刻的各质点的相对位移反应?x?i?1和速度增量?x?x?i?1??x?i???x???i?1??x??i???x???x??i?1代入结构运动微分方程??i?1；??(5)最后将?x?i?1和?x，求得第ti?1时刻的相对加速度反应x??i?1?-(??g??xx?,i?1??i?1???????i?x?i?1)?????C?i?x?1?1可以用MATLAB进行编程分析，画图位移时程图，详细程序见附录。
2.3 杜哈梅积分
杜哈梅积分在考虑阻尼的情况是：
v?e???t[v(0)cos?Dt?v(0)???v(0)?Dsin?Dt]?1m?D?t0e???(t??)p(?)sin?D(t??)d?
可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。
3. 位移时程反应对比分析
利用 MATLAB 将理论解法，杜哈梅积分， Wilson-θ 法求解出来的位移时程反应画在同一张图 中，进行比较分析。
从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在 5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。
4. 结论
本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于 MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用 Wilson-θ 法）进行相互验证，从最后的位移分析图对比上，可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。
附录：MATLAB 源程序
%理论解，杜哈梅积分，Wilson-θ法程序 clc; clear
h1=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:i
canshu1(j)=p(j)/(m*w)*bb*sin(w*(t(i)-t1(j)));%杜 哈梅积分中的被积函数 end
y(i)=sum(canshu1);%%位移值 end
for i=1:aa/bb-1
v1(i)=(y(i+1)-y(i))/bb;%计算速度 end
for i=1:(aa/bb-2)
a(i)=(v1(i+1)-v1(i))/bb;%计算加速度 end hold on
plot(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:ndzh
kxin=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]; end
plot(t,-wyimt/3000,'b'); hold off
legend('\\fontsize{9}\\fontname{ 黑体} 理论解','\\fontsize{9}\\fontname{黑体} 杜哈梅积分法 ','\\fontsize{9}\\fontname{黑体} wilson-{\\theta}法')
%基于双线性滞回模型的单自由度体系的