结构动力学作业

结构动力学作业

用El Centro地震波记录数据计算地震速

度与位移

学生姓名

学号

指导教师

学院土木工程学院

完成时间 2012年9月

El Centro 地震波记录数据计算地震速度与位移

EL Centro 地震波是世界上第一条成功记录全过程数据的地震波，对于人类地震的研究有着重大的意义。 EL Centro 地震波记录了地震过程中加速度随时间的变化情况，现需要研究根据已知的记录数据，计算出地震过程速度和位移随时间的变化情况，并绘制出时程图。

该数据每隔0.02s 记录一次地震加速度的数值，在地震发生的54s 的时间内，共记录下2676个数据。用MATLAB 根据数据绘制的加速度时程图如下：

图1 EL Centro 地震波加速度时程图

加速度经过一次积分可以获得t 时刻的速度，速度经过一次积分可以获得t 时刻的位移。所以通过加速度求速度和位移就是对加速度关于时间的函数对时间求积分的问题，在几何上即是求加速度图像与坐标轴围成的图形的几何面积。

由于加速度数据是以一系列离散数据的形式给出，难以求出加速度关于时间的函数的解析表达式，于是考虑采用数值积分的方法，构造合适的算法，求出积分的近似值。

设:

1()b

a I f x dx =⎰ (1) 2()b

a I g x dx =⎰ (2) 在高等数学中知道，当|()()|f x g x ε-<时，12||()I I

b a ε-<-。这说明，当ε充分小的时候，可以用2I 近似地代替1I 。所以，求任意函数()f x 在[a ,b ]上的定积分时，如果难以使用解析的方法求出()f x ，则可以寻找一个在[a ,b ]与()f x 逼近，但形式上却简单且易于求积分的函数()g x ，用()g x 在[a ,b ]上的积分值近似地代替()f x 的积分值。一般选择被积函数的插值多项式充当这样的替代函数。

选择的插值多项式的次数不同，就形成了不同的数值积分公式。选择一次多项式时，称为梯型公式，选择二次多项式时，称为辛普森(Simpson)公式。[1]

通过插值节点0x a =和1x b =作线性插值函数1()g x ，利用1()()f x g x ≈，得梯型公式：

1()()[()()]b b b a a a x b x a f x dx g x dx f a f b dx a b b a

--≈=+--⎰⎰⎰ (3) 若()f x 用通过节点0x a =，1x b =，2x c =的二次多项式2()g x 代替，得辛普森公式：

0201122012010210122021()()()()()()()()[()()()]()()()()()()

b

b b a a a x x x x x x x x x x x x f x dx g x dx f x f x f x dx x x x x x x x x x x x x ------≈=++------⎰⎰⎰ (4)

使用MATLAB 编制计算程序，分别用梯型公式和辛普森公式进行数值积分，结果如下：

梯型公式：t=54s 时刻（地震结束时刻）速度v=40.48mm/s ，位移s=2510.1mm,时程图为：

图2 EL Centro 地震波速度时程图（梯型公式）

图3 EL Centro地震波位移时程图（梯型公式）

辛普森公式：t=54s时刻（地震结束时候）速度v=73.17mm/s，位移

s=3772mm,时程图为：

图4 EL Centro地震波速度时程图（辛普森公式）

图5 EL Centro地震波位移时程图（辛普森公式）

根据数值分析代数精度相关研究可知[2]，梯型公式和辛普森公式的代数精度分别为1次和3次，对于充分光滑的函数，辛普森公式的计算精度要比复合梯型公式高很多。

观察可知，梯型公式和辛普森公式得出的时程图走势基本相同，结果相差不大，证明了积分方法和结果的正确性。综合考虑代数精度，辛普森公式得出的结果更加精确可信，可以用来指导工程实践。

参考文献

[1]刘卫国，MATLAB程序设计与应用. 北京：高等教育出版社，2006

[2]同济大学计算数学教研室，现代数值计算. 北京：人民邮电出版社，2003

[3]RoyR.Craig,Jr著，常岭，李振邦翻译，结构动力学. 北京：人民交通出版社，1995.12

附录

附录一梯型公式MATLAB计算程序：

clear

load AB A B

v=[];s=[]; %v速度向量s距离向量

for m=1:2676

if m==1

I=0;

elseif m==2

I=0.02*(0.5*B(1)+0.5*B(2));

else

I=0.02*(0.5*B(1)+sum(B(2:m-1))+0.5*B(m)); end

v=[v I];

end

C=v;

for m=1:2676

if m==1

J=0;

elseif m==2

J=0.02*(0.5*C(1)+0.5*C(2));

else

J=0.02*(0.5*C(1)+sum(C(2:m-1))+0.5*C(m)); end

s=[s J];

end

figure;

plot(A,B);

title('加速度时程图');

legend('加速度0.1cm/s/s');

xlabel('t');

ylabel('a');

grid;

figure;

plot(A,v);

title('速度时程图');

legend('速度0.1cm/s');

xlabel('t');

ylabel('ｖ');

grid;