关于古塔变形的数学研究

2013高教社杯全国大学生数学建模竞赛

承 诺 书

我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模

竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮

件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的

成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。

如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行

公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）：

所属学校（请填写完整的全名）： 贵州电力职业技术学院 参赛队员 (打印并签名) ：

指导教师或指导教师组负责人 (打印并签名)： 江静 （论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以上内容请仔细核对，提交后将不再允许做任何修改。如填写错误，论文可能被取消评奖资格。）

日期：2013 年9 月 15 日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2013高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号（由赛区组委会评阅前进行编号）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

基于古塔变形的分析

摘要

古塔建筑，其形态千奇百态，巍巍壮观。在中国，有的塔是一个地区或一个城市的标志，如杭州雷锋塔、西安大雁塔等。古塔代表劳动人民的聪明智慧、是建筑史艺术的结晶，本文针对古塔历史久远，古代建筑以土木砖为主要材料，经历了千年的风风雨雨，古塔的整个塔身已经发生了不同程度的变形，诸如倾斜、弯曲、扭曲等。本文主要针对附件1所提供的4次观测的数据及古塔受长时间承受自重、气温、风力等外界环境因素各种作用进行综合分析，从统计数据的分析结果加深对古塔变形的认识，旨在确定古塔各层中心位置的通用方法，分析该塔倾斜、弯曲、扭曲等变形情况及该塔的变形趋势。

针对问题一，首先以古塔的第一层为研究对象，把古塔俯视图抽象成一个正八边形，其次以某点为坐标原点，在正八边形的平面上建立平面直角坐标系XOY，用平面坐标运算的方法来求古塔中心位置的坐标的通用方法。类似地，根据通用方法，可以依依地求出各层的中心位置坐标。在此建模过程中，利用CAD绘图软件绘出3个八边形在平面直角坐标系不同位置的图。还利用MATLAB软件根据附件1提供的数据通过编程绘制了古塔所有点的分布图。最后通过分析和坐标计算，得出了该古塔各层的中心位置坐标。

针对问题二，该问题主要让我们分析该塔倾斜、弯曲、扭曲等变形情况。考虑到倾斜、弯曲、扭曲这些情况的变化必将带来整个塔身各个点在空间位置上的变化。因此，我们利用问题一中的模型结果（古塔各年的中心坐标）或《表一》， 《表二》数据，通过MATLAB软件编程来绘制各层中心坐标的空间分布图像，在该图像初步分析的基础上，还利用CAD绘图软件绘制了古塔的主视图，更加直观反映出四次测量的该塔倾斜、弯曲、扭曲等变形情况。观察图像可以得出结论：2009年和2011年变化情况基本不变，1986年和1996年古塔偏移量较小。而1996到2009年偏移大小明显。

针对问题三，该问题目的是要分析该塔的变形趋势。要研究塔的变形趋势，可以根据附表1中的数据，在四次测量各层8个的位置都不变的情况下，利用MATLAB软件将各次测量各层的点1坐标绘制成函数图像，以及将各次测量的古塔塔顶点1坐标点绘制成函数图像，然后再结合这些函数图像的变化规律来研究该古塔的变形趋势。通过函数变化规律的研究，得出了如下结论：该古塔的高度有逐年下降的变化趋势，古塔逐渐向y轴的负方向上变化，古塔逐渐向x轴的正方向上变化，但古塔在各个方向上的变化是不明显的微弱变化。

建立数学模型的过程中运用了EXCLE，MATLAB，CAD，Word等相关软件。

关键词：变化趋势 平面直角坐标系 坐标运算 中心位置坐标 函数图像模型 几何模型

一．问题的重述

我国历史悠久，是世界上著名的文明古国，在历史的长河中留下了许许多多的古文明建筑，其中一些古塔记录了我国在古代建筑的辉煌成就。在我国古代，很早就有建塔的历史，而且在不同的时间和空间，古塔的风格也不一样。每一座塔的背后都会有当年的历史事件，对于研究我国不同朝代，不同地方的历史具有重要的意义。但是，对于一座已经经历了上千年风雨的古塔都发生各种的变化，包括形状，大小及各点空间的偏移。

一些古塔由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的环境因素影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。因此，保护我国的古塔建筑已经是一项重要的工作，更不容忽视的历史责任和义务。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。

某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

请你们根据附件1提供的4次观测数据，建立数学模型讨论以下问题： 问题一： 给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

问题二：分析该塔倾斜、弯曲、扭曲等变形情况。

问题三：分析该塔的变形趋势。

二．问题的分析

本文主要目标是分析和利用附件1中的数据来确定古塔各层中心位置的通用方法，并且分析该塔倾斜、弯曲、扭曲等变形情况及该塔的变形趋势。首先，通过所建立的通用方法求出4次测量的古塔各层中心坐标；其次，对各次测量的统计数据分析，利用MATLAB软件对各种情形下具有相关规律的数据绘制相应的函数图像，建立函数图像模型，在利用该模型来分析该古塔倾斜，弯曲，扭曲等变形情况。最后，结合问题一及问题二的模型和结论，分析该古塔的变形趋势。鉴于此目的，针对本文具体3个问题，可以进行如下分析：

2.1 针对问题一的分析

本问题主要是在古塔俯视图的平面上建立平面直角坐标系XOY，为了研究和建立模型的需要，首先将古塔俯视图抽象成正的八边形。这个正八边形在该平面直角坐标系上移动，移动两次产生了三种不同位置，相应地，八边形上的各点坐标也会随着整个图形的移动而发生变化。在此问题上，只要利用附表1中古塔各次测量及各层中的点1坐标，点4坐标，点8坐标进行各个状态时的相关坐标计算。

2.2 针对问题二的分析

本问题要求分析该塔倾斜，弯曲，扭曲等变形的情况，这座已经经历上千年历史的古塔，可以想象古塔全身各个部位在空间里都会发生不同程度的位置偏移，这也是该塔倾斜，弯曲，扭曲等变形的原因和反映。

为了更加直观的研究该塔的各种变形情况，可以利用附件1中四次测量的统计数据，结合问题一所得的表中各次测量的古塔各层中心坐标数据，再利用MATLAB数学绘图软件对《表一》和《表二》数据绘制出四次测量古塔各中心坐标变化规律的三维空间图像，然后对所绘制的图像观察和对比。最后，根据该中心坐标变化规律图像进行分析后，得出该塔倾斜，弯曲，扭曲等变形的情况。

2.3 针对问题三的分析

本问题的目的是为了分析该古塔的变形趋势，可以想象，古塔在空间里的变化无疑是整个古塔的塔身的各个部位在变化，本质上是整个塔的各位置的空间坐标点随时间的变化情况，即所有点在x轴，y轴，z轴三个方向的偏移。

为了利用附件1中的数据来研究塔变形趋势的问题，基于假设3：假设四次测量的各层的各个点保持不变以及基于假设8：假设各层相对应的一系列点的组合对于研究问题的作用是相同的（例如第一层到第十三层中的点1构成的点集和点2构成的点集）。因此，我们可以任意选取各年份测量的一系列点集为研究对象，再利用MATLAB软件把这4个年份的点集绘制成函数图像，通过函数图像的变化趋势及4个年份相应的函数图像的比较分析来得出古塔的变形趋势。

在本文所建立的模型中，选取了各次测量的古塔各层点1横坐标及纵坐标作为研究对象。但考虑到竖直方向上的高度变化因素，在这里可以选取各次测量的塔顶中的点1坐标作为研究对象，然后利用MATLAB软件绘制4个年份所测量的塔顶点1坐标的函数图像。最后，在结合这三个函数图像来分析古塔的变形趋势。

三、模型的假设及符号说明

3.1 模型的基本假设：

1.假设每一层塔的竖直方向的截面为梯形

2.假设计算时带来的误差忽略不计

3.假设四次测量的各层的各个点保持不变

4.假设不受其他环境因素的影响，只受题目给定因素的影响

5.假设受地震、飓风的影响变化不大

6.假设题目缺失的数据对数据计算影响不大，可忽略不计

7.假设题目所给的数据真实可靠

8.假设各层相对应的一系列点的组合对于研究问题的作用是相同的（例如第一层到第十三层中的点1构成的点集和点2构成的点集）

3.2 模型的符号说明：

（c,d): 状态1时古塔4次测量的第一个点的平面坐标。

（e,f): 状态1时古塔4次测量的第四个点的平面坐标。

（g,h): 状态1时古塔4次测量的第八个点的平面坐标。

（n,m): 状态2时古塔4次测量的第四个点的平面坐标。

（u,v): 状态2时古塔4次测量的第八个点的平面坐标。

（a,b): 状态2时古塔4次测量的中心位置点的平面坐标。

（a1,b1): 状态3时古塔4次测量的中心位置点的平面坐标。

Zi: 表示古塔各层的第i个点的高度坐标。

Xi: 表示古塔各层的第i个点的横向坐标。

Yi: 表示古塔各层的第i个点的纵向坐标。

K: 古塔随时间变迁及各种环境因素影响产生变形的系数（当环境变化不大且时间较短时，k为常数）

→：表示某种变化趋势的负方向变化

←：表示某种变化趋势的正方向变化

四． 模型的建立与求解

4.1 针对问题一模型的建立与求解

由问题和数据的分析可知，该模型是通过建立平面直角坐标系以及相关坐标运算的初等数学几何模型。该模型主要是分析8边形在平面坐标系上做必要的移动后，产生3种不同的位置状态（状态1，状态2，状态3）的各个状态的相关坐标点变化，以及进行的坐标的运算，然后得出古塔各层中心位置的计算通用方法，最后利用该通用方法逐一求出各次测量的古塔各层中心坐标。

建立模型的过程中，主要运用了CAD绘图工具以及MATLAB软件。

利用MATLAB软件画出古塔各个点连接的立体图：

图1.

利用MATLAB软件画出古塔各个点连接的俯视图：

图2

模型的准备一：(中心位置坐标通用方法的建立）

在古塔某层的俯视图所在的平面上建立平面直角坐标系，如下图，该图可以称为正8边形移动的初始位置，即上面所说的状态1，针对该图有以下相关的坐标描述。

图3

在状态1时：

点1的坐标表示为（c,d);

点4的坐标表示为（e,f);

点8的坐标表示为（g,h);

在状态一的基础上，将八边形平移，使得点1（c,d)与坐标的原点O(0,0)重合,得出状态二，显然其他个各点的坐标也相应的改变，如下图4。

图4

在状态2时：

点1的坐标表示为(0.0)

点4的坐标表示为（n,m);

点8的坐标表示为（u,v);

其中：n=e-c m=f-d;

u=g-c v=h-d;

中心位置坐标：

横坐标 a=(n+u)/2=(e+g-2c)/2;

纵坐标 b=(m+v)/2=(f+h-2d)/2;

在图4（状态2）的基础上，已经建立了在此状态时的任一层中心坐标的通用公式。为了计算在初始状态时的任一层中心坐标，我们必须把该图经过平移到原来的位置，因此，产生了状态3，如下图所示。

图5

在状态3时，古塔某层的中心位置坐标：

横坐标 a1=a+c=(e+g-2c)/2+c=(e+g)/2;

纵坐标 b1=b+d=(f+h-2d)/2+d=(f+h)/2;

模型的准备二

针对古塔各层中心位置在竖直方向上通用方法的确定，可以利用平均值求解的方法，即在z轴方向上，将4次测量各层的8个数据相加，再求取平均值。数

学表达式如下：

z=(z1+z2+z3+z4+z5+z6+z7+z8)/8

=(zi)/8;

i8

综上所述，可以得出确定古塔各层中心位置的通用方法，该数学模型的表达式即为：

（(e+g)/2，(f+h)/2，(z1+z2+z3+z4+z5+z6+z7+z8)/8）

模型一：

利用模型的准备所建立的通用方法，根据1986年所观察的数据，可以建立该古塔第一层中心位置坐标。

图6

图6可视为第一层时的初始位置（状态一），相关的坐标如下：

点4的坐标（563.782,518.108）；

点1的坐标（565.454,528.012）；

点8的坐标（569.5,527.356）；

由状态1变为状态2如下图所示：

图7

图7可视为第一层时的坐标平移，点1坐标平移到坐标系原点O时（状态二），各个相关点新坐标如下：

利用通用方法，可以求出点4的坐标（-0.672，-9.904）； 点1的坐标（0,0）；

点8的坐标（4.046，-0.656）； 再根据状态2的点4和点8的坐标，利用直线中点公式，可以求出中心坐标 (x.y)=(1.714,-5.28)

由状态2变为状态3如下图所示：

图7

根据通用公式可以得复原后：横坐标x=566.66475 纵坐标y=522.7105

竖坐标z=1.787375

根据1986年所观察的数据得出古塔第一层中心坐标为： （x, y，z)=(566.66475,522.7105,1.787375）

因此，综合模型准备一和二的中心位置的通用方法以及模型一的求解过程，再根据1986年，1996年，2009年和2011年四次测量的统计数据可以得出各年古塔各层的中心坐标，列表如下：

4.2 针对问题2的模型建立与求解

4.2.1 问题的背景分析

中国是世界上古塔最多的国家之一。隋唐之前多为木塔，但因易燃易蚀，山西应县木塔是硕果仅存的一座。其余保存下来的多为砖、石、铁塔及明清后期的琉璃塔。“木塔的框架结构都用榫卯（实木建筑的连接部件）连接，不用钉子，因此允许一定的变形，不容易倾斜。而除此之外的其他塔，受到地基的不均匀沉降、地震、大风等影响，都会有可能倾斜。[1]

倾斜是指基础两端点倾斜方向的沉降差与其距离的比值，建筑中心线或其墙、柱等，在不同高度的点对其相应底部点的偏移现象。对于一座古塔的倾斜亦是如此。

弯曲，即受到力的作用而造成形变，这种力的作用是合力最终形成的结果，纸板在制造过程中有多种作用力的存在。如三层板有浪面和贴合面纸，五层有A浪、B浪贴合与贴合面纸，这几种工艺形成产品过程中，其最终所形成的平衡作用力，即形变力。

4.2.2 问题的重要性分析

据文物部门统计，在全国3000多座古塔中，斜塔有近20座。比如西安大雁塔向西北倾斜1.01米；洛阳道士塔倾斜15度；上海护珠塔向东偏移2.28米；辽宁前卫镇塔向东北偏移1.7米等。[2]

4.2.3 模型的建立

模型一：

利用MATLAB软件所绘制的1986年，1996年，2009年以及2011年古塔各层中心坐标在空间的分布图，如下图：（函数图像模型）

图 8. 模型一的求解：

根据如上图8，可以得出初步结论：1986年和1996年的两年该古塔中心坐标在空间的分布是基本吻合的，变化不大，因此，可以认为从1986年到1996年这十年间该塔倾斜，弯曲，扭曲等变形的情况基本没有变化。2009年到2011年的两三年间，该古塔中心坐标在空间的分布也是基本吻合的，变化也是不大，可以认为没有发生倾斜，弯曲，扭曲等变形。然而，从1996年到2009年这13年间古塔中心坐标在空间的分布却发生了比较大变化，可以认为该古塔在这段时间内已经发生一定程度的倾斜，弯曲，扭曲等变形。究其原因，我们推测这段时间很大可能受到地震或者飓风等自然因素的影响导致其发生各种变形。

模型二：（几何模型）

根据上图的图像分析，我们利用CAD绘图工具绘制了该古塔在这四次测量数据下的古塔偏移情况的主视图。在如下图中，可以更直观的发现古塔的倾斜，弯曲，扭曲等变形情况。

图 9. 模型二的求解：

图9中‘0年中心轴线’表示为该古塔未发生倾斜，弯曲，扭曲等变化时各层的中心坐标连线。‘塔在x,y轴上的偏移’可以作为为扭曲程度的评价指标。

在图8的基础上，在结合上图9，我们可以直观地看出该古塔在这四年的变形情况。图9中，2009年和2011年变化情况基本不变，1986年和1996年古塔偏移量较小。而1996到2009年偏移大小明显。

模型准备一：

下面的图10表示四次测量年份的古塔在x轴方向上的倾斜情况图；图11表示为四次测量年份的古塔在y轴方向上的倾斜情况图；图12表示为古塔由于塔身发生一定倾斜后古塔的高度在各个测量年份的变化情况。

图10

图11

图

12

组合图13

组合图13的目的是为了研究该古塔的四次测量年份分别在坐标系上的x轴，y轴，z轴方向上各个年份之间相对偏移的角度，其中各个相对角度具体的计算过程如下所述：

在x轴方向上有

1986年与1996年的相对偏移角度为 α1=arccos(567.2473/567.2544)=0.0050; 1996年与2009年的相对偏移角度为 α2=arccos(567.2544/567.336)=0.0170; 2009年与2011年的相对偏移角度为 α3=arccos(567.336/567.3376)=0.0024;

在y轴方向上有

1986年与1996年的相对偏移角度为 α4=arccos(522.2387/522.2438)=0.0044; 1996年与2006年的相对偏移角度为 α5=arccos(522.2148/522.2438)=0.0105; 2009年与2011年的相对偏移角度为 α6=arccos(522.2135/522.2148)=0.0022;

在z轴方向上有

1986年与1996年的相对偏移角度为 α7=arccos(55.11975/55.12325)=0.0113; 1986年与1996年的相对偏移角度为 α8=arccos(55.091/55.12325)=0.0342; 1986年与1996年的相对偏移角度为 α9=arccos(55.087/55.09)=0.0104;

弯曲度=0.01863 倾斜度=0.01863

X轴上的扭曲度=0.00813 Y轴上的扭曲度=0.0057 Z轴上的扭曲度=0.01863

4.3 针对问题三的模型建立与求解

模型一：

要分析该古塔的变形趋势，首先选用1986年，1996年，2009年以及2011年各层上的点1的横坐标x，利用MATLAB软件绘制每一年的x坐标在平面的分布情况的函数图像和各个年份之间x坐标的变化情况，建立函数图像模型，如下图10所示：

图 14.

根据图 10，可以直观的得出一些初步的结论：1986年的变化曲线与1996年以及2009年的变化曲线与2011年的变化曲线基本重合，这一点可以表明1986年到1996年间以及2009年到2011年间该古塔基本没有发生任何变形。

但是，通过观察该图像很明显可以看出由两条分开的曲线，可以表明1996

年到2009年这13年可能受到某些不利因素的影响，比如环境因素（台风，地震等）使该古塔突然发生较大的变形。

对于1986年和1996年的曲线可以看出，x值的波动较小，基本没有发生任何变形；但对于2009年和2011年的曲线可以看出，x值的波动较大，且有向着x轴正方向偏移的趋势，因此此时该古塔已经发生了一定的倾斜，弯曲或扭转的变形情况。

模型二： 同理，要分析该古塔的变形趋势，首先选用1986年，1996年，2009年以及2011年各层上的点1的纵坐标y，利用MATLAB软件绘制每一年的y坐标在平面的分布情况的函数图像和各个年份之间x坐标变化情况，建立第二个函数图像模型，如下图11所示：

图 15.

同理，根据图 11.，可以直观的得出一些初步的结论：1986年的y值变化曲线与1996年的以及2009年的y值变化曲线与2011年的变化曲线基本重合，这一点可以表明1986年到1996年间以及2009年到2011年间该古塔基本没有发生任何变形。

但是，通过观察该图像很明显可以看出由两条分开的曲线，这表明1996年到2009年这13年可能受到某些不利因素的影响，比如环境因素（台风，地震等）使该古塔突然发生较大的变形。

由1986年和1996年的曲线可以看出，y值的波动较大，基本存在一定程度的变形；但由2009年和2011年的曲线可以看出，y值的波动较小，且有向着y轴负方向偏移的趋势，因此此时该古塔在y轴方向已经有慢慢复原的变化趋势。

模型三：

在分析该古塔变形趋势时，考虑到高度方向（z轴方向）的影响，在这里可以选取各次测量的塔顶中的点1坐标作为研究对象，然后利用MATLAB软件绘制4个年份所测量的塔顶点1坐标的函数图像如图 12.所示：

图 16

根据图12，该古塔的顶点位置相对于地面是呈逐年下降的趋势，只表明该古塔经历千年历史，经历了无数的风雨和阳光暴晒等各种环境因素的相互作用，塔身材料慢慢松软，导致塔身有下降的趋势。

在x轴方向上，塔身的各个部位也慢慢地发生了偏移，偏移方向为x轴的正方向。

在y轴方向上，塔身各部位向着y轴的负方向偏移。 从整体上看，塔身的各个部位在x,y,z方向的偏移量都比较小，这表明该古塔存在各种微妙的变形趋势。

变迁的影响，产生形变，影响系数k的作用下变为坐标Xi。向X轴正的方向变化。

Yi→Yk：该数学表达式表示古塔从建成开始的Y坐标因受到环境以及时间变迁的影响，产生形变，影响系数k的作用下变为坐标Yi。向Y轴负的方向变化。

Zi→Zk：该数学表达式表示古塔从建成开始的Z坐标因受到环境以及时间变迁的影响，产生形变，影响系数k的作用下变为坐标Zi。向Z轴负的方向变化。

五．模型的分析

6.1 假设的合理性分析

本题涉及到三个最为关键假设：3.假设四次测量的各层的各个点保持不变。

7.假设题目所给的数据真实可靠。8.假设各层相对应的一系列点的组合对于研究问题的作用是相同的（例如第一层到第十三层中的点1构成的点集和点2构成的点集）。

对于假设3，只有假设四次测量的各层的各个点保持不变，才能对问题三建立一个关于一系列相关点在这四年的变化规律的函数模型，可见，该假设是建立模型的基础，该假设合理。

对于假设7，为了建立能更好研究该古塔问题相关的数学模型，都必须是利用每年测量统计下的真实可靠的数据，因此，该假设合理。

对于假设8,只有当各层相对应的一系列点的组合对于研究问题的作用是相同的话，才能在各层8点中任意取一点作为研究对象，建立相关模型，也是建立模型的基础，该假设合理。

六．模型的推广

其实这个模型和综合分析的模型相似，我们运用了坐标法、函数图像去分别比较古塔倾斜和其他因素的关系，这个模型可以推广到更多和更复杂的相关问题中和评估问题中，比如对水域沉桩的测量，对土方的测量计算，对深基坑的水平位移监测，曲线桥墩的中心位置评估等。

19

七．模型的评价与优化

1．在建立模型时，我们利用统计数据来找绘图相关的点，再利用MATLAB软件绘制出函数图像，利用直观的函数图像模型来研究古塔的各种变化情况。

2．通过对统计数据的处理，找出问题的规律从而找出函数，这样的模型能运用于以后的相似问题中去解决问题。

3．我们评估的数据可能会带来误差，影响到计算的结果，但是模型还是比较合理的。

4．多取几组数据比较，具有代表性。

5．该模型未能体现到数据间的关系，应该尽可能的做到精确，以便更直观更充分的体现出函数模型，并应该进行更多数据的检验。

参考文献：

[1]：罗哲文：中国古建筑专家

[2]：http://it.sohu.com/20081218/n261294145.shtml 2013-9-14

[3]：作者：颜文勇，书名：数学建模

[4]：作者：同济大学数学系，书名：工程数学线形代数第五版 出版社：高等教育出版社

20