基于MATLAB实现对结构动力响应的几种算法的验证

1、 基于 MATLAB 实现对结构动力响应的几种算法的验证算例 首先，本文给出一算例， 结构在外力谐振荷载 P(t) = P0 sint 作用下，分别利用理论解法，杜哈梅积分， Wilson- 法求出该结构的位移时程反应。其中：m = 3.5103 kg ， P0 = 1.0104 N ， k =1.3584515107 ，=0.05 ，=52.3s1 ，=62.3s1 ， = 1-2 =62.222 ，初始位移、速度v(0) = 0 ，v(0) = 0 ；算法验证 2.1 理论解法 运动方程为： mv+cv+kv=sin由线性代数解出其理论解为:由于初始位移v(0) =0 ，v(0) =0 ；

2、则： v(t ) =1.052698986.230cos62.222t 18.106sin 62.222t +2.012808757 1146sin 52.3t 325.829cos52.3t 可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。 2.2Wilson-法 Wilson-法是Wilson于1966年基于线性加速度法的基础上提出一种无条件收敛的计算方法。该方法假定在时程步长内，体系的加速度反应按线性变化。对于地震持续时间内的每一个微小时 段 ，从第一时段开始到最后一个时段，逐一的重复以下计算步骤，即得到结构地震反应的全过程。下面以第i+1时段（）为例：2.3 杜哈梅

3、积分 杜哈梅积分在考虑阻尼的情况是： 可以用 MATLAB 进行编程分析，画图位移时程图，详细程序见附录。 3. 位移时程反应对比分析利用 MATLAB 将理论解法，杜哈梅积分， Wilson- 法求解出来的位移时程反应画在同一张图 中，进行比较分析。 从图中可以看出，以上三种方法得出来的位移时程曲线基本吻合，误差基本保持在 5%以内，所以以上几种方法在求解相关问题上都具有一定的作用效力。结论 本文通过一个简单的单自由度系统动力分析算例（仅作位移分析，其它分析雷同），基于 MATLAB，将理论解法，杜哈梅积分法，逐步积分法（本文采用 Wilson- 法）进行相互验证，从最后的位移分析图对比上，

4、可以很好的看出三种方法均能很好的彼此验证，从而说明了三种方法在相关问题上的作用效力。附录：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*10(-6)*(1146*sin(52.3*t)-325.829*cos(52.3*t); plot(t,v,k) hold on %杜哈梅积分法aa

5、=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;%计算速度

6、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

7、 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=kxindpxin; dtjsdu=6*dxtao/(ct*(ct2*dt2)-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1; jsdu=jsdu1+dtjsdu; dtsdu=(dt/2)*(jsdu+jsdu1); sdu=sdu1+dtsdu; dtwyi=dt*sdu1+(1/3)*dt2*jsdu1+(d

8、t2/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(fontsize9fontname 黑体 理论解,fontsize9fontname黑体 杜哈梅积分法,fontsize9fontname黑体 wilson-theta法)%基于双线性滞回模型的单自由度体系的地

9、震能量分析程序 %质量57041kg，阻尼36612 Ns/m,初始刚度2350000N/m,刚度折减 系数0.2，屈服位移0.01m，采用ELCENTRO波%参数替换直接在下面修改，然后运行clc format 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

10、,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表示正向弹性f

11、or ii=1:num if pd=0 %弹性阶段 k=k1; if x(ii)Xy pd=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)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/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/

12、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)-Xy pd=-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)

13、 if isreal(d(j)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/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*d

14、tx; hfl(ii)=-k1*Xy+dtf; a(ii)=(P(ii)-hfl(ii)-c*v(ii)/m; end end i f pd=1 %正向弹塑性 k=k2; if v(ii)pxmax pd=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)=1 dth=d(j); end end hp=h-dth;vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/h/2); ap=a(ii-1)+

15、(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)0 pd=-2;

16、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)=1 dth=d(j); end end hp=h-dth; xp=x(ii-1)+v(ii-1)*dth+a(ii-1)*dth2/2+(a(ii)-a(ii-1)*dth3/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)

17、+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); end end if pd=-2 %正向弹性 k=k1; if 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

18、(b); for j=1:length(d) if isreal(d(j)=1 dth=d(j); end end hp=h-dth; vp=v(ii-1)+a(ii-1)*dth+(a(ii)-a(ii-1)*dth2/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

19、; 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; end end i f ii=num break; end kd=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)=

20、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-1 EI(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

21、)=EH(jj)+Eh(jj+1); end end for j=1:num time(j)=h*(j-1); end fid=fopen(EL波弹塑性各种能.txt,wt); %输出结果到EL波弹塑性各种能.txt fprintf(fid,%g %g %g %g%gn,time;EI;ED;EH;EK); fclose(fid); plot(time,EI); hold on plot(time,ED); hold on plot(time,EK); hold on plot(time,EH); 附录资料:不需要的可以自行删除c语言典型问题处理方法小结循环问题（1）、数论问题1、求素数 fo

22、r(i=2;i1，如果它仅有平凡约数1和a，则我们称a为素数(或质数)。整数 1 被称为基数，它既不是质数也不是合数。整数 0 和所有负整数既不是素数，也不是合数。 2、求最大公约数和最小公倍数a、 if(ab) t=a; a=b; b=t; for(i=a;i=1;i-) if(a%i=0&b%i=0) break; printf(largest common divisor:%dn,i); printf(least common multiple:%dn,(a*b)/is);b、辗转相除法求解 a1=a; b1=b; while(a%b!=0) t=a%b; a=b; b=t; print

23、f(largest common divisor:%dnleast common multiple:%d,b,a1*b1/b);3、求完数 一个数如果恰好等于它的因子之和，这个数就称为“完数”。 例如：6的因子为1、2、3，而6123，因此6是“完数”。for(a=1;a=1000;a+) s=0; for(i=1;i=a) break; if(s=a) printf(%dt,a);注意S=0所放的位置 4、分解质因数 将一个整数写成几个质因数的连乘积，如： 输入36，则程序输出36=2*2*3*3 。解一、看似简单，但要自己完整地写出来还真不容易！竟然还动用了goto语句，正好可以熟悉一下g

24、oto语句的用法！main() int a,z,i; clrscr(); scanf(%d,&a);判断下一个数开始有要重新从2开始了。所以用loop语句回到for语句，这是for语句仍从2初始化。从2开始的原则不变，变的是a的值。 loop: for(z=2;z=a;z+)判断是否为质数for(i=2;i=z;i+) if(z%i=0) break;判断是否为a的质因数 if(z=i) if(a%z=0) k+; if(k=1) printf(%d=%d,a1,z);用计数器来解决每行输入不同的问题。 else printf(*%d,z); a/=z; goto loop; 解二：main(

25、) int n, k=2, isfirst=1; printf(Input n=); scanf(%d,&n); while(k=n) if(n%k=0) if(isfirst) printf(%d=%d, n, k); isfirst=0; else printf(*%d,k); n/=k; else k+; printf(n);5、从键盘输入两个整数，输出这两个整数的商的小数点后所有1000位整数 for(i=1;i=2;i-) if(fm%i=0&fz%i=0) fz/=i; fm/=i; z=fz/fm; fzx=fz%fm; if(fzx=0) printf(%d%d/%d-%d%d

26、/%d=%dn,z1,fz1,fm1,z2,fz2,fm2,z); else if(z=0) printf(%d%d/%d-%d%d/%d=%d/%dn,z1,fz1,fm1,z2,fz2,fm2,fzx,fm); else printf(%d%d/%d-%d%d/%d=%d%d/%dn,z1,fz1,fm1,z2,fz2,fm2,z,fzx,fm);（2）近似问题1、书P122习题4-6注意千万不要忘记添加#include “math.h”#include math.hmain() float x,j=1,k,s,so; int n; scanf(%f,&x); s=x; so=x+1; f

27、or(n=1;fabs(s-so)1e-6;n+) for(k=1;k1e-6) x=(x1+x2)/2; f=x*x*x+4*x*x-10;可以用/*if(f*f10) x2=x; else x1=x; printf(%fn,x);（3）枚举法（4）数列问题二、数组问题（1）排序问题1、从小到大排序main() int a10,i,j,t; for(i=0;i10;i+) scanf(%d,&ai); for(i=1;i10;i+) for(j=0;jaj+1) t=aj+1;aj+1=aj;aj=t; for(i=0;i10;i+) printf(%d ,ai); printf(n);注意

28、排序问题：1、须迅速，熟练，无差错经常插入在程序中间2、现使用最大数下沉冒泡法还可以使用最小数上浮冒泡法3、j控制前面一个数和后面一个数一一比较。由于是最大数下沉，i+1后j仍要从0开始。4、i控制这样的操作一共要做多少次5、注意i j的控制次数2、从大到小排序main()现使用最大数上浮冒泡法还可使用最小数下沉冒泡法 int a10,i,j,t; for(i=0;i10;i+) scanf(%d,&ai); for(i=1;i=i;j-) if(ajaj-1) t=aj-1; aj-1=aj; aj=t; for(i=0;i10;i+) printf(%d ,ai);（2）二维数组三、字符或

29、字符串输入输出问题（1）字符打印1、打印*此类题的溯源为书P122 4.11（1），其他题都是它的拓展 for (i=1;i=n;i+) 一共要输出的行数 for(j=1;j=i;j+) 每行要打印的*数 printf(*); printf(n); a、*解题要点：此类题关键在于找到每行要打印的个数和行数的关系。此题j=i j=n-i+1b、* for(i=1;i=n;i+) 一共要输出的行数 for(j=1;j=n-i;j+) 控制空格数 printf( ); for(k=1;k=i;k+) 每行要打印的*数 printf(*); printf(n); c、 * * *解题要点：在出现空格的

30、时候，在找到每行要打印的*个数和行数的关系后，还应找到空格和行数的关系，分不同的参数进行循环。此题k=i j=n-i j=i-1k=n-i+1d、* * for(i=1;i=n;i+) for(j=1;j=n-i;j+) printf( ); for(k=1;k=2*i-1;k+) printf(*); printf(n); * *e、 * * * for(i=1;i=n-1;i+) for(j=1;j=i;j+) printf( ); for(k=1;k=2*(n-1-i)+1;k+) printf(*); printf(n); for(i=1;i=n;i+) for(j=1;j=n-i;j

31、+) printf( ); for(k=1;k=2*i-1;k+) printf(*); printf(n); * * * * * *2、打印9*9乘法表解题要点：注意寻找行与列的规律。i*ji代表列j代表行for(i=1;i=9;i+) for(j=1;j=9;j+) printf(%-3d ,i*j); 注意输出格式的控制 printf(n); 3、九九乘法表1 2 3 4 5 6 7 8 92 4 6 8 10 12 14 16 183 6 9 12 15 18 21 24 27 9 18 27 36 45 54 63 72 814、杨晖三角形11 11 2 11 3 3 11 4 6

32、4 11 5 10 10 5 11 6 15 20 15 6 1（2）字符串打印问题for(i=1;i=7;i+) ai1=1; aii=1; for(i=3;i=7;i+) for(j=2;j=i-1;j+) aij=ai-1j-1+ai-1j; gets(a); puts(a); for(i=1;i0;j-) aj=aj-1; a0=t; for(k=0;k=a&ai=z) ai-=32; puts (a);3、逆序输出gets (a); c=strlen(a); for(i=0;i=0;i-) 蓝色部分可以简写为绿色部分coutai-1;4、如输入：ab1 3,;z 输出：ab1注意点：

33、1、= =2、while语句的使用处体会3、全面考虑问题 3,;zgets(a); while(a0= ) for(i=0;ai!=0;i+) ai=ai+1; for(i=0;ai!=0;i+)if(ai= &ai+1!= ) printf(n); else if(ai= &ai+1= ) for(k=i;ak!=0;k+) ak+1=ak+2; i-; elseprintf(%c,ai);5、输入3个字符串，按从小到大排序输出这3个字符串。 使用一个两维数组贮存多个字符串char a8181;注意：如何使用一个两维数组贮存多个字符串 int i,j; for(i=0;i3;i+) gets

34、(ai); for(i=0;i3;i+) puts(ai);注意：1、scanf(%d%s,&n,str) 其中%s为字符串格式2、逐个给字符串赋值的方法见书140页。 不可for（i=0；ai!=0;i+）3、stri=stri-A+10;4、pow函数5、任何进制转为十进制的方法6、输入一个整数n和一个字符串str，计算并输出n进制数str的值。 如输入：7 16则输出：13(16)7=(13)10如输入：16 3A则输出：58(3A)16=(58)10#include stdio.h#include math.hmain() char str81; int n,i,s=0,t; clrs

35、cr(); scanf(%d%s,&n,str); for(i=0;stri!=0;i+) if(stri=A) stri=stri-A+10; else stri=stri-0; t=strlen(str); for(i=0;stri!=0;i+) s+=strt-i-1*pow(n,i); printf(%d,s);编写程序，将一个十进制正整数转换成十六进制数。 注意类比#include main()char a20;int x,i=0,j;clrscr();scanf(%d,&x);while(x) if(x%16=10&x%16=0;j-)printf(%c,aj);printf(n)

36、;7、输入一个字符串，将其中的缩写形式展开，并输出展开后的该字符串。所谓展开缩写形式就是将其中由大小写字母或数字构成的形如a-f、U-Z、3-8 的形式展开成为 abcdef 、UVWXYZ 、345678，若出现f-a、A-7、9-5等形式则不予理睬。例如： 输入：qwe246e-hA-d$-%4-7A-Dz-xp-R4-0输出：qwr246efghA-d$-%4567ABCDz-xp-R4-0main() char a81; int i,c,s,k,t; gets(a); for(i=0;ai!=0;i+) if(ai=-) if(ai-1=A&ai+1=a&ai+1=0&ai+1i;k-

37、)ak+c-2=ak;as-1+c-2+1=0; for(;i=t;i+) ai=ai-1+1; puts(a);补充：循环：求：a+aa+aaa+.的值#includevoid main()int a,n,i=1,sn=0,tn=0;coutinput a and nan;while(i=n)tn=tn+a;sn+=tn;a*=10;i+;coutthe answer is snendl;两个乒乓球队进行比赛，各出3人。甲队为A,B,C;已对是X,Y,Z；已经抽签决定比赛名单。有人向队员大厅比赛的名单。A说他不和X比，C说他不和X,Z比。请编程序找出3对赛手的名单。#includevoid

38、main()char i,j,k;for(i=X;i=Z;i+)for(j=X;j=Z;j+)if(i!=j)for (k=X;k=Z;k+)if(i!=k&j!=k)if(i!=X&k!=X&k!=Z)coutA-i B-j C-kendl;枚举口袋中有红，黄，蓝，白，黑5种颜色的球若干。每次从口袋中任意取出3歌，问得到3种不同颜色球的可能取法，输出每种排列的情况。#include#include /在C语言中不用加这句void main()enum colorred ,yellow ,blue,white, black;color pri;int i,j,k,n=0,loop;for(i=

39、red;i=black;i+)for(j=red;j=black;j+)if(i!=j)for (k=red;k=black;k+)if(k!=i)&(k!=j)n+;coutsetw(3)n; /setw是输出格式的限定for(loop=1;loop=3;loop+)switch(loop)case 1:pri=color(i);break;case 2:pri=color(j);break;case 3:pri=color(k);break;default:break;switch(pri)case red:coutsetw(8)red;break;case yellow:coutsetw(8)yellow;break;case blue:coutsetw(8)blue;break;case white:coutsetw(8)white;break;case black:coutsetw(8)black;break;default:break;coutendl;couttotal:nendl;数组：不用strcat函数,编写程序，将两个字符串串接起来。#includeusing namespace std;voi