古塔的变形模型

1、2013高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、我网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设

2、置报名号的话）： sm063 所属学校（请填写完整的全名）： 山东现代职业学院 参赛队员 (打印并签名) ：1. 马昱轩 2. 李倩 3. 梁帅健 指导教师或指导教师组负责人 (打印并签名)： 宋祖芳 日期：2013年 9月 16 日赛区评阅编号（由赛区组委会评阅前进行编号）：sm063032013高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：古塔的变形模型摘要某古塔是我国重点保护文物,已有上千年历史。由于长时间

3、承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生诸如倾斜、弯曲、扭曲等各种变形。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。基于附件1提供的4次观测数据：对于问题1，1986年和1996年的观测数据中都缺少13层一个点的数据（因MATLAB程序中用的是循环语句，所以计算时赋予0值），其它各层均给出8个点的观测值，为使所得数据更具真实性，确定古塔各层中心位置的方法更具适用性，本文将每层所给8点构成的图形看做不规则八边形，用中垂线求

4、交点法求得古塔各层中心坐标。对于问题2，结合问题1的分析，采用垂直投影法1求古塔的倾斜度，根据所得数据，分析古塔的倾斜程度（因1986年和1996年13层赋予值后所得数据偏差较大，为使所得数据更具真实性，所以本问题起1986年和1996年13层数据予以舍弃）；弯曲是建立在二维平面上的一条曲线，通过截取古塔过轴、轴的界面，求出每层古塔的倾斜度，从而分析得出古塔塔身在轴、轴的界面的弯曲程度。同理，也可分析轴、轴的界面的弯曲程度；扭曲同样是采用垂直投影法1求古塔每层的倾斜度，建立三维立体空间，根据所得数据，分析古塔塔身的扭曲程度。对于问题3，利用问题1、2所得数据，进行合理的分析与猜想，进而分析出古

5、塔塔身的倾斜、弯曲、扭曲等变化趋势。本文的模型解决了题目给出的问题，计算过程中充分尊重观测数据，给出更符合实际的结果。本文所得结果大部分由图表给出，结合图像，较为直观地表现出古塔变形情况。结果表明，采用数学软件可得出可靠结论。关键词：古塔变形监测 垂直投影法 倾斜度一、问题的重述由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年

6、3月对该塔进行了4次观测。请你们根据附件1提供的4次观测数据，讨论以下问题：1. 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。2. 分析该塔倾斜、弯曲、扭曲等变形情况。3. 分析该塔的变形趋势。二、问题的分析2.1问题的背景某古塔是我国重点保护文物。已有上千年历史，由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生诸如倾斜、弯曲、扭曲等各种变形。对现存历史文物古塔的保护，掌握古塔的形变，显得尤为重要，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年

7、3月和2011年3月对该塔进行了4次观测。2.2问题的分析 基于附件1提供的4次观测数据： 2.2.1问题一的分析 本题建立模型是利用中垂线（即中垂线相交于一点为中点)法，求解古塔各层中心的坐标，以确保古塔中心坐标的准确性。 2.2.2问题二的分析 根据层与层之间的中心点，塔体倾斜变化采用垂直投影法1(即垂线法，利用相对位移量进行监测)；弯曲是建立在二维平面上的一条曲线，通过截取古塔过轴、轴的界面，求出每层古塔的倾斜度，从而分析得出古塔塔身在轴、轴的界面的弯曲程度。同理，也可分析轴、轴的界面的弯曲程度；扭曲同样是采用垂直投影法求古塔每层的倾斜度，建立三维立体空间，根据所得数据，分析古塔塔身的扭

8、曲程度。 2.2.3问题三的分析 我们在对古塔进行安全性监测方面要满足两个方面：塔体的倾斜变化；塔身各段的变化1。针对上述几个方面就采取对应的监测方法。3、 问题的假设1. 假设求解古塔各层中心坐标时纵坐标先不予以考虑；2. 假设古塔的地基土抗压能力强、韧性大、不易收缩，对古塔的变形不产生任何影响；3. 假设对古塔观测数据的缺失可忽略不计；4、 符号及文字说明表示建立古塔坐标轴的横坐标表示古塔各层的第(=1,2,.,8）个点的横坐标表示建立古塔坐标轴的纵坐标表示古塔各层的第（=1,2，.,8)个点的纵坐标表示建立古塔坐标轴竖坐标表示古塔各层的第(=1,2,.,8)个点的竖坐标表示古塔的第条直线

9、的斜率表示两个点之间的中点表示两个点之间的中点表示古塔中心点到水平面的距离表示古塔中心点的位移量表示古塔的倾斜度表示古塔的第(=1,2,.,13)层数（MATLAB程序中）表示古塔各层的第(=1,2,.,8）个点（MATLAB程序中）表示关于横坐标的矩阵（MATLAB程序中）表示关于纵坐标的矩阵（MATLAB程序中）表示各层每两点之间中垂线的斜率（MATLAB程序中）表示每层8条中垂线两两相交所得交点的关于横坐标的矩阵（MATLAB程序中）表示每层8条中垂线两两相交所得交点的关于纵坐标的矩阵（MATLAB程序中）五、模型的建立与求解（1） 问题一的模型建立与求解1. 通过中垂线的方法对古塔各层

10、中心坐标进行求解，根据附件1表格数据，1986年和1996年的观测数据中都缺少13层一个点的数据计算时赋予0值。然后观察数据，古塔的每层八个点的高度变化不大，假设把古塔的每一层的点看作是在一个平面上，求解古塔各层中心坐标时纵坐标先不予以考虑，之后用平均值的方法计算出。 分析过程如下图1-1：（用MATLAB软件作图（附件1，程序1） 1-1 根据上图，做出任意两边的中垂线，交与一点，如下图所示1-2（在上图上用画图工具作图）： 1-2由图可知，点坐标为，点坐标。根据公式，得: ，；， 根据直线的斜率公式，求出直线的斜率。同理，求出直线的斜率。 再根据中垂线定理，求出中垂线，的斜率，； 设中垂线

11、的直线方程为： 设中垂线的直线方程为： 通过方程、建立方程组: 解得： ； 再把、式中所求的数据代入中，即可解出两中垂线交点坐标。同理，通过数学软件MATLAB编程求解（附录1，程序2）可得出每层八条中垂线两两相交的八个交点坐标，得出数据并对这八个中垂线的交点求平均值，因1986年和1996年的观测数据中都缺少13层一个点的数据计算时赋予0值得出的数据偏差太大，此时予以舍弃。塔顶数据是塔顶的四个点求平均值而得。各层中心坐标如下表：表一： 古塔层次各层中心坐标1566.522.2566.522.3566.522.4566.522.5566.522.6566.522.7566.522.8566.5

12、22.9567.522.10567.522.11567.522.12567.522.13-塔顶567.24725522.24375 而每层古塔都有一定的高度差，因此，假设每层古塔的中心高度为每层八点的平均高度，即：，因此，古塔每层的中心坐标是：表二：1986年古塔层次古塔各层中心坐标1566.522.1.2566.522.7.320253566.522.12.755254566.522.17.078255566.522.21.72056566.522.26.7566.522.29.8566.522.33.9567.522.36.10567.522.40.11567.522.44.12567.5

13、22.48.13-塔顶567.24725522.2437555.表三：1996年古塔层次古塔各层中心坐标1566.522.1.7832566.522.7.3566.522.12.750754566.522.17.5566.522.21.766566.522.26.22957566.522.29.8566.522.33.9567.522.36.10567.522.4760940.11567.522.44.12567.522.48.13-塔顶567.25435522.2366555.11975表四：2009年古塔层次古塔各层中心坐标1566.522.1.76452566.522.7.3093566

14、.522.12.732254566.522.17.069755566.522.21.6566.522.26.2117566.991522.29.8567.522.33.9567.522.36.8437510567.522.40.11567.522.44.12567.522.48.6997513567.522.52.塔顶567.336522.214855.091表五：2011年古塔层次古塔各层中心坐标1566.522.1.763252566.522.7.29053566.81276522.12.4566.522.17.0255566.522.21.6566.522.26.20457566.991

15、8522.29.8178567.522.33.9567.522.36.8222510567.522.40.11567.522.44.12567.522.48.13567.522.52.塔顶567.3375522.213555.087（2） 问题二的模型建立与求解1，根据问题一求出的每层的中心坐标数据，通过分析古塔的倾斜、弯曲、扭曲及层与层之间的中心点，确定了塔体倾斜变化采用垂直投影法1(即垂线法，利用相对位移量进行监测);弯曲、扭曲即为塔身各段的变化采用求倾斜度的方法。用MATLAB数学软件编程（附件1，程序3）和画图工具画图，如下图2-1： 2-1 通过对上图的观察与分析，每层的中心大致在一

16、条倾斜的直线上，采用垂直投影的方法，对古塔的倾斜角进行求解，如下图通过MATLAB数学软件编程和画图工具并用作图（附件1，程序4），建立的数学模型，如下图2-2； 2-2 通过对图形的观察，点为古塔塔顶的中心坐标，点为古塔第一层的中心坐标，点是点在点所在水平面上的投影， 根据倾斜度公式,3计算出的值： 根据以上公式，得出以下数据，如表六：倾斜度198619962009201171.70.69.69.比较数据： 分析数据，的值越小，倾斜程度越严重，古塔随着年代的增加而越来越倾斜。2. 对于弯曲的情况，弯曲是建立在二维平面上的一条曲线，在解决问题的过程中，为了对弯曲便以分析，用轴、轴所在的平面建立

17、二维平面的数学模型，（由于古塔的弯曲方向具有不确定性，此方法具有一定的局限性）通过画图工具画出简易二维图形2-3； 2-3 根据倾斜度公式：3，求出每层的倾斜度，每层塔的的高度： 中心点的位移距离： 根据上述公式计算出1986、1996、2009、2011年的每层古塔的倾斜度，由于1986、1996年的第十三层丢失数据，因此舍去第十三层的斜率，用MATLAB数学软件编程（附件1，程序5），计算结果，如下表：层次年代1986199620092011倾斜度197.95.159.159.297.96.164.164.398.97.161.159.498.98.165.166.592.90.49.49

18、.6108.107.107.107.7106.105.66.66.8101.99.66.65.9101.99.59.59.1075.75.97.97.1176.76.97.97.1276.75.82.82.13-49.49. 根据表格中所得的数据，可以分析出，每层古塔的倾斜度不同，塔身的弯曲程度也不同，倾斜度越小，塔身的弯曲程度越严重，通过这四个年代每层对应的倾斜度的比较，得出古塔的塔身的弯曲程度每年都在加深。3. 扭曲等情况，需要按着上述方法计算出每层古塔的倾斜角，通过比较每层的倾 斜角，可得出古塔的扭曲情况。下图2-4为古塔每年的中心坐标图（通过MATLAB数学软件作图（附件1，程序6），

19、旋转选取适当的视角）： 2-4 通过对上图的观察，每层都有一定的扭曲程度，采用第1小题中的方法求出每层古塔的倾斜度，如图（用画图工具作简易图）： 根据倾斜度公式：，求出每层的倾斜度，每层塔的的高度： 中心点的位移距离： 根据上述公式计算出1986、1996、2009、2011年的每层古塔的倾斜度，由于1986、1996年的第十三层丢失数据，因此舍去第十三层的斜率，用MATLAB数学软件编程（附件1，程序7），计算结果，如下表：层次年代1986199620092011倾斜度181.80.123.122.281.80.127.128.382.81.123.120.482.81.127.127.57

20、8.76.44.44.696.94.90.90.794.93.57.57.891.89.57.56.993.92.42.42.1061.61.71.71.1162.61.71.71.1235.35.60.60.13-26.26.观察图表中数据，（3） 问题三的模型建立与求解 古塔存在时间久远，经过长时间大自然的侵袭和破坏，并且还受到人们生产生活的影响。塔的内部和外部都受到一定的影响。根据建筑力学得知，当发生地震时，地基结构发生了改变使得古塔发生不均匀沉降，导致古塔发生了变形。经过四次观测数据得出古塔每次都发生了很微小的变形，如建筑物超过自重，会使得抗外界干扰能力及承载力降低从而造成破坏。就风压

21、而言，相当于是对古塔产生额外的与重力相同的外力作用，会增加内部古塔柱梁板等构件的内力，可能造成构件内力超过其承载能力，从而造成古塔变形。就风吸而言，建筑上一般不考虑其对古塔减轻重力的有利一面，而且需要古塔在风的吸力作用下是否会造成风吸力大于建筑物屋面重力的情况，而造成屋面被掀起。这样的例子，主要表现在建筑物外围雨棚的设计中，对于一般屋面，其外形一般产生的只是风压。因此会出现古塔的倾、弯曲、扭曲等各种变形。根据以上问题得出的数据，分析得出的古塔的倾斜、弯曲、扭曲等变化情况，如果不受外力影响（除重力外），塔身会随着时间的变化会继续倾斜、弯曲、扭曲等变化，如图2-4，根据四次测量的数据及分析，猜想古

22、塔会继续向某些方向倾斜、弯曲、扭曲等变化；如果考虑到外力对塔身的影响，塔身会发生不规则倾斜、弯曲、扭曲等变化。如果古塔到了一定的年代，塔身的变形程度过大，古塔重心不足以维持塔身的平衡，古塔将会出现倒塌一部分，或者完全倒塌。 6、 模型的评价与改进方向6.1模型的优缺点（1） 本模型综合利用了、中垂线法、垂直投影法等数学方法，数学推导严谨，理论性强，所得数据真实可靠。（2）模型具有一般性，在一定条件下适用于其他古建筑物的变形情况监测。（3） 由于数据较多，计算量过大，对于问题一，算得八条垂直平分线的八个交点后，采用极限思想的计算会使数据更准确一些。（4） 问题二中的弯曲，扭曲在已得的数据上用拟合

23、方法会更优化一些。6.2模型的改进放向我们的模型通过一些合理的假设，使得问题的数学描述较为简单直观，但实际上还有很多因素与要讨论的问题密切相关，而且是应该加以考虑的。对于问题一的求解方法，如果多次计算每次求出的数据，根据极限思想，也就是利用求极限方法求各层中心坐标，最终的结果更无限接近于真实值。如问题二，在解决弯曲变化时，应多方位考虑弯曲方向等问题，这些都是在研究最优策略时应该认真考虑的问题，因为时间关系，我们只讨论了比基础、较简单的建模方法，还有待于完善和改进。7、 参考文献1梁海奎.古塔变形测量方法探讨.城市勘探2011年第03期2软件编程；清华大学出版社3韩煊,李宁.地铁施工引起的建筑物

24、扭曲变形分析.土木工程学报2010年1月8、 附录附件1程序1：x=565.454562.058561.39563.782567.941571.255571.938569.5565.454;y=528.012525.544521.447518.108517.407519.857523.953527.356528.012;plot(x,y)程序2：（1986年）xa=load(xa.txt);ya=load(ya.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 fo

25、r i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j

26、,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data1.xls,xx,1)xlswrite(data2.xls,yy,1)（1996年）xa=load(xb.txt);ya=load(yb.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)

27、-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data3.xls,xx,1)xlswrite

28、(data4.xls,yy,1)（2009年）xa=load(xc.txt);ya=load(yc.txt);X=xa;Y=ya;K=zeros(13,8);M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1

29、);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data5.xls,xx,1)xlswrite(data6.xls,yy,1)（2011年）xa=load(xd.txt);ya=load(yd.txt);X=xa;Y=ya;K=zeros(13,8);

30、M=zeros(13,8);N=zeros(13,8);x=X,X(:,1);y=Y,Y(:,1);for j=1:13 for i=1:8 M(j,i)=(x(j,i)+x(j,i+1)/2; N(j,i)=(y(j,i)+y(j,i+1)/2; K(j,i)=-(x(j,i+1)-x(j,i)./(y(j,i+1)-y(j,i); endendM N K xx=zeros(13,8);yy=zeros(13,8);N1=N,N(:,1);M1=M,M(:,1);K1=K,K(:,1);for j=1:13 for i=1:8 xx(j,i)=(N1(j,i+1)-N1(j,i)+K1(j,

31、i)*M1(j,i)-K1(j,i+1)*M1(j,i+1)/(K1(j,i)-K1(j,i+1); yy(j,i)=K1(j,i)*xx(j,i)+N1(j,i)-K1(j,i)*M1(j,i); endendxxyyxlswrite(data7.xls,xx,1)xlswrite(data8.xls,yy,1)程序3：x=566.566.566.566.566.566.566.566.567.567.567.567.567.24725;y=522.522.522.522.522.522.522.522.522.522.522.522.522.24375;z=1.7.3202512.7552

32、517.0782521.720526.29.33.36.40.44.48.55.;plot3(x,y,z,r.)程序4：（1）x=565.454562.058561.39563.782567.941571.255571.938569.5565.48562.238561.663564.001567.995571.165571.801569.414565.506562.415561.931564.216568.048571.076571.666569.33565.526562.555562.144564.387568.091571.005571.558569.263565.548562.70656

33、2.373564.571568.136570.929571.443569.191565.57562.854562.6564.752568.18570.857571.333569.121565.671563.132562.883564.949568.172570.679571.094568.994565.77563.403563.158565.141568.164570.506570.862568.87565.868563.674563.433565.333568.156570.333570.63568.747565.961563.927563.693565.516568.148570.1715

34、70.408568.631566.078564.193563.958565.649568.094570.013570.236568.615566.195564.459564.224565.782568.039569.854570.063568.598566.308564.716564.481565.91569.701569.897568.582567.255567.235567.247567.252;y=528.012525.544521.447518.108517.407519.857523.953527.356527.764525.364521.42518.226517.563519.96

35、1523.908527.141527.52525.188521.394518.343517.716520.063523.864526.93527.327525.047521.373518.435517.838520.144523.829526.762527.119524.896521.351518.534517.969520.232523.791526.581526.915524.748521.329518.632518.095520.315523.755526.406526.652524.585521.356518.846518.346520.441523.672526.167526.397

36、524.427521.382519.055518.59520.564523.591525.933526.141524.268521.408519.263518.834520.686523.51525.701525.9524.12521.433519.462519.068520.801523.433525.482525.628523.95521.463519.607519.242520.885523.35525.259525.355523.78521.492519.753519.415520.969523.268525.037525.092523.616521.521519.893521.055

37、23.188524.822522.238522.242522.251522.244;z=1.7921.8181.7831.7691.7721.771.7941.8017.3267.3517.3147.3017.3067.3047.3247.33612.76112.78612.74912.73612.74112.7412.75812.77117.08417.10917.07217.05917.06417.06317.08117.09421.72621.75121.71421.70121.70521.70821.72321.73626.26726.30926.30826.26426.18926.1

38、3626.16426.24429.86929.91129.9129.86629.79129.73729.76529.84633.38333.42533.42433.3833.30533.25133.27933.3636.88736.92936.92836.88436.80936.75536.78336.86440.20140.21440.24440.22340.17140.03840.12940.15744.47244.48544.50544.48644.44244.30944.444.42848.74348.75648.77648.75748.71348.5848.67148.69952.8

39、6652.87852.89752.8852.70352.79452.82255.12855.10855.12855.129;plot3(x,y,z,b.)（2）z=-1:1;x=565.454562.058561.39563.782567.941571.255571.938569.5565.454;y=528.012525.544521.447518.108517.407519.857523.953527.356528.012;z=000000000; plot3(x,y,z,r)程序5：（1986年）i=1:13;xa=566.566.566.566.566.566.566.566.567.

40、567.567.567.567.24725;za=1.7.3202512.7552517.0782521.720526.29.33.36.40.44.48.55.;x=xa;z=za;Q=zeros(1,13);for i=1:12 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,1)（1996年）i=1:13;xa=566.566.566.566.566.566.566.566.567.567.567.567.567.25435;za=1.7837.12.7507517.21.7626.229529.33.36.40.44.48.

41、55.11975;x=xa;z=za;Q=zeros(1,13);for i=1:12 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,2)（2009年）i=1:13;xa=566.566.566.566.566.566.566.991567.567.567.567.567.567.567.336;za=1.76457.30912.7322517.0697521.26.21129.33.36.8437540.44.48.6997552.55.091;x=xa;z=za;Q=zeros(1,13);for i=1:13 Q(i)=(z

42、(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,3)（2011年）i=1:13;xa=566.566.566.81276566.566.566.566.9918567.567.567.567.567.567.567.3375;za=1.763257.290512.17.02521.26.204529.81733.36.8222540.44.48.52.55.087;x=xa;z=za;Q=zeros(1,13);for i=1:13 Q(i)=(z(i+1)-z(i)/(x(i+1)-x(i);endQ=Qxlswrite(bb.xls,Q,4)

43、程序6：（1986年）x=566.566.566.566.566.566.566.566.567.567.567.567.567.24725;y=522.522.522.522.522.522.522.522.522.522.522.522.522.24375;z=1.7.3202512.7552517.0782521.720526.29.33.36.40.44.48.55.;plot3(x,y,z)（1996年）x=566.566.566.566.566.566.566.566.567.567.567.567.567.25435;y=522.522.522.522.522.522.522.5

44、22.522.522.47609522.522.522.23665;z=1.7837.12.7507517.21.7626.229529.33.36.40.44.48.55.11975;plot3(x,y,z)（2009年）x=566.566.566.566.566.566.566.991567.567.567.567.567.567.567.336;y=522.522.522.522.522.522.522.522.522.522.522.522.522.522.2148;z=1.76457.30912.7322517.0697521.26.21129.33.36.8437540.44.48.6997552.55.091;plot3(x,y,z)（2011年）x=566.566.566.81276566.566.566.566.9918567.567.567.567.567.567.567.3375;y=522.522.522.522.522.52