土木工程数值分析（共8篇）
班级： 姓名： 学号： 指导老师：
某路基工程施工过程数值模拟
本文首先对FLAC3D软件进行了介绍，简明阐述了其特点、应用范围及不足，然后结合具体路堤工程，采用FLAC3D软件对施工过程进行了模拟，生成了初始竖向和水平应力云图、第一次填筑及填筑结束时的沉降云图及水平位移云图;最后生成了路基中心点和坡脚节点的沉降曲线。
关键词：FLAC3D:数值模拟；应力云图；沉降云图；位移云图
1、FLAC3D的功能与特性
自R.W.Clough 1965年首次将有限元引入土石坝的稳定性分析以来，数值模拟技术在岩土工程领域获得了巨大的进步，并成功解决了许多重大工程问题。特别是个人电脑的出现及其计算性能的不断提高，使得分析人员在室内进行岩土工程数值模拟成为可能，也使得数值模拟技术逐渐成为岩土工程研究和设计的主流方法之一。数值模拟技术的优势在于有效延伸和扩展了分析人员的认知范围，为分析人员洞悉岩体、土体内部的破坏机理提供了强有力的可视化手段。FLAC系列软件的出现，为岩土工程研究工作者提供了一款功能强大的数值模拟工具。
1.1 FLAC3D主要特点
FLAC(Fast Lagrangian Analysis of Continua)是由Itasca公司研发推出的连续介质力学分析，是该公司旗下最知名的软件系统之一，FLAC目前已在全球七十多个国家得到广泛应用，在国际土木工程（尤其是岩土工程）学术界和工业界享有盛誉。
FLAC3D界面简洁明了，特点鲜明。其使用特征主要表现为：命令驱动模式、专一性、开放性。作为有限差分软件，相对于其他有限元软件，在算法上，FLAC3D有以下几个优点：采用“混合离散法”来模拟材料的塑性破坏和塑性流动，比有限元中通常采用的“离散集成法”更准确、合理；即使模拟静态系统，也采用动态运动方程进行求解，这使得FLAC3D模拟物理上的不稳定过程不存在数值上的障碍；采用显示差分法求解微分方程。采用FLAC3D进行数值模拟时，必须指定有限差分网格、本构关系和材料特性、边界和初始条件，这是FLAC3D求解的一般流程。
1.2 FLAC3D的应用范围
FLAC3D的应用范围已拓展到土木建筑、交通、水利、地质、核废料处理、石油及环境工程等领域，成为这些专业领域进行分析和设计不可或缺的工具。其研究范围主要集中在以下几个方面：
●岩体、土体的渐进破坏和崩塌现象的研究；
●岩体中断层结构的影响和加固系统（如喷锚支护、喷射混凝土等）的模拟研究；
●岩体、土体材料固结过程的模拟研究；
●岩体、土体材料流变现象的研究；
●高放射性废料的地下存储效果的研究分析；
●岩体、土体材料的变形局部化剪切带的演化模拟研究；
●岩体、土体的动力稳定性分析、土与结构的相互作用分析以及液化现象的研究等。
1.3 FLAC3D的不足
毋庸置疑，FLAC3D是十分优秀的岩土工程数值模拟软件，其实用性和专业性得到了广泛证实。但FLAC3D也存在着诸多不足，主要集中在以下几个方面：
●求解时间受网格尺寸影响很大； ●某些模式下的计算求解时间很长： ●前处理功能较弱。
●FLAC3D对于复杂三维模型的建立仍然十分困难。
2、某路基工程施工过程模拟
2.1工程概况
该路基工程地基计算深度为50m，分为两层，上部为回填土，厚度10m，下部为粘土，厚度40m:路基计算宽度为200m，填筑高度为5m，坡度为1：1.5。各土层物理、力学参数见表1。
2.2模型建立
由于几何模型具有对称性，可以采用1/2模型进行分析。首先建立坐标系，坐标系的原点O设置在地基表面与模型对称轴的交点，水平向右为X方向，竖直向上为Z方向，垂直于分析平面的方向为Y方向。网格的建立按照分区域建模的思路进行，如图2所示。由于路基坡脚的位置存在一个关键点，所以将模型划分成5个矩形区域，对每个区域按照控制点利用brick单元建立网格，并进行分组后赋值。考虑到网格尺寸的一致性，Y方向只设置一个单元，该方向单元尺寸为5m。
网格建立以后，首先设置边界条件。对底部边界节点的X、Y、Z三个方向的速度进行约束，相当于固定支座，对X两侧的边界进行水平速度约束。由于Y方向只设置一个单元长度，所以对模型中所有节点的Y方向速度均进行约束，相当于进行平面应变分析。
2.3初始应力计算
在路基施工前，需要将路基部分网格赋值为空模型，而将地基部分的网格赋值为Mohr模型。由于实例中存在null模型，不能采用solve elastic的求解方法获得初始应力，所以采用分阶段的弹塑性求解方法。先将Mohr模型的凝聚力c值和抗拉强度δt赋值为无穷大进行求解，保证在重力作用下单元不至于发生屈服，然后再将Mohr模型参数赋值为真实值，再进行求解。
图3和图4为初始应力计算结束时得到的竖向应力和水平应力云图，可以发现最大竖向应力值为85.3kPa，最大水平应力值为42.0kPa，静止侧压力系数约为0.5，与理论计算值基本一致。
2.4施工过程模拟
在进行路基施工模拟前要将初始应力计算过程中产生的节点位移和速度进行清零处理。因本工程路基高度为5m，高度方向共划分了5个单元，为了模拟路基填筑的施工过程，采用分级加载的方法激活路基单元，每次激活1m高度的单元，相当于每次填土高度为1m，分5次填筑完成，每次填土进行一次求解。
图5和圈6分别为第一次填筑结束时的沉降云图和水平位移云图，图7和图8分别为填筑结束时的沉降云图和水平位移云图，可以发现最大沉降发生在地基表面的左侧边界处，而最大水平位移发生在坡脚以下的深部地基中。
图9给出了计算过程中不平衡力的收敛过程，在初始应力计算过程中有很大的数值逐渐收敛，在后续路基填筑过程中，每一次填筑引起的不平衡力在计算过程中逐渐收敛。路基填筑计算过程中监测了路基中心点和路基坡脚处节点的沉降和水平位移。图10为监测结果，可以看出在迭代过程中节点位移随迭代步数的变化。
2.5绘制沉降曲线
在工程应用中常常需要得到某些关键点的沉降曲线，本工程中主要关心的是路基中心节点和坡脚节点的变形结果。利用FLAC3D中自带的print命令配合log文件可以得到这两个节点的变形量。用记事本打开相应的log文件，拷出沉降数据，就可在Excel中生成沉降曲线，如图11所示。
3、结论 本文结合具体路堤工程，采用FLAC3D软件对施工过程进行了模拟，生成了初始竖向和水平应力云图，得出最大竖向应力值为85.3kPa，最大水平应力值为42.0kPa，静止侧压力系数约为0.5，与理论计算值基本一致；第一次填筑及填筑结束时的沉降云图及水平位移云图，可以发现最大沉降发生在地基表面的左侧边界处，而最大水平位移发生在坡脚以下的深部地基中；最后生成了路基中心点和坡脚节点的沉降曲线。
参考文献
随着我国国民经济的发展和人们居住环境要求的提高, 高层建筑得到了快速的发展, 城市土地资源变得日趋紧张, 因而地下建设成为一个重要发展方向, 如:地下室、地下商业街以及地铁车站等, 而各种形式的地下建筑工程的建设面临的首要问题就是大规模的开挖, 因此, 深基坑开挖和支护将成为工程技术人员日益关注的话题。如何保证深基坑的开挖安全以及支护结构的稳定, 降低基坑开挖对周边建筑的影响成为工程界共同需要探究的问题。
1 围护结构水平位移监测
基坑监测中的重要一项是围护结构的水平位移监测, 这能反映出基坑开挖是否安全施工的直接表现, 围护结构变形是基坑施工过程中工程人员的关注重点。随着基坑开挖时间的增长, 围护结构水平位移也不断的增加, 当开挖到坑底时, 围护结构发生水平位移呈现缓慢的增大趋势。要减少水平位移蔓延, 要尽量减少基坑的暴露时间, 快速的将基坑垫层与底板等主体结构进行施工。在这个时间段内, 加大监测力度, 保证基坑的整体安全, 当底板浇筑结束后, 围护结构的水平位移也基本结束。
2 地层沉降监测
深基坑开挖会引起周边建筑物的沉降, 这样有可能对人们的正常生活带来安全隐患和经济损失, 为了保障周边建筑以及周围管线的安全, 必须对地层沉降进行监测, 以便于基坑开挖安全进行。
一般情况下, 在基坑开挖初期和中期周围沉降较为稳定, 在开挖和支撑过程中略有上下波动, 主要原因是在开挖过程中周围地表发生沉降。当继续开挖到后期, 可能会发现明显的沉降趋势, 这时候就要加强监测的频率, 随着基坑底板浇筑完毕之后, 基坑周围沉降的趋势就会趋于收敛状态, 这说明基坑基本处于稳定安全的状态。所以, 在施工期间对基坑周围管线沉降进行信息化监测非常必要, 有利于及时的采取安全措施来预防事故的发生, 从而保证基坑周围建筑的安全。
3 地下水位监测
直接反映基坑局部是否发生漏水及坑内降水效果的是深基坑外水位的监测, 这关系到整个基坑安全稳定。一般情况, 在基坑开挖初期, 随着基坑开挖不断的加深, 对坑内水位进行降水处理, 坑外各监测点水位变化成整体下降趋势。
4 监测数据分析与预测
获取各种监测数据后, 需要进行及时的处理, 最好当天处理, 但是难免还是会出现误差, 粗差、系统误差以及偶然误差等, 这些误差对后续变形分析和解释带来了困难, 容易导致错误的判断和预测, 这样监测便失去了意义。如何消除或者减弱误差, 提高监测精度, 就是及时采用计算机对现场监测的数据进行整理和分析。
1) 数据整理。通过现场监测采集的原始数据按照一定的数据分布情况显示出来, 呈现出数字特征值, 然后进行离群数据的取舍。
2) 插值法。依据已经监测到的数据进行分析研究, 采用函数的方法算出由测量结果得出的规律而实际上并没有监测的数据。
3) 对实际监测的数据采用统计的方法进行分析, 拟定出一种能够反应监测数据变化规律及趋势的函数关系式, 对下一阶段的施工进行预测, 预测周边地表沉降量, 预测结构物的安全性等, 据此对工程采取技术措施和对策。
5 现场监测注意事项
首先监测过程应该加强巡逻, 看监测点是否完好, 由于施工现场人员较多且复杂, 所以对监测元件也要做好监督和保护;其次要求监测人员对基坑的稳定情况有一个直观的把握, 特别是钢筋混凝土的支撑力, 地面是否开裂, 连续墙是否出水, 这些有时不一定在监测数据中体现, 但是在巡视的过程中可能会发现。再次, 为了提高监测结果的准确性, 必须坚持温度监测与日常正常监测工作同时进行, 避免温度异样的影响。再其次, 测斜观测时必须坚持慢和准, 这样做既保证测量数据的准确又可以保证使用工具的损耗, 益处多多。最后, 观测水位时, 应该在正常的天气环境下, 不要在雨后测量, 这么做能够保证数据的正确性。
6 数值模拟分析
根据工程的基坑支护结构, 将模拟分为三个施工步骤:第一工况开挖到围护桩顶部;第二工况开挖到第二道环梁底部, 第三工况开挖到基底。然后进行模拟结果分析。通常会得出以下结论:对于形状不规则的基坑, 双环梁支护体系是一种非常好的基坑支护方式, 因为其受力合理且施工方便, 对周围的环境影响也较小。而基坑开挖是的周围的土体产生了较大的水平位移, 而且位移会锁着距基坑边距离的增大而减小, 并且距离基坑越远, 减小的速率就越快。
土体及围护桩最大位移位置基本都处于基坑开挖附近, 且位移随着开挖深度的增大而增大, 而且各阶段的水平位移和深度变化规律不一样, 需要具体情况具体分析。通常模拟分析的结果和实测结果如出一辙, 但是影响位移和内力因素具有复杂性, 假定的条件和实际条件存在着差异和变异。
7 结语
实践表明, 基坑开挖过程中, 需要对整个支护结构、周边土体以及邻近结构物和管线等进行系统的监测, 才能对工程的安全情况全面了解和把握, 确保工程施工顺利进行。对深基坑工程施工进行实时监测并进行监测数据分析, 不仅可以发现问题, 还可以保证工程的安全, 通过对数值模拟分析, 对基坑工程、施工设计以及应急措施具有一定的参考和指导意义。
参考文献
[1]GB 50497-2009, 建筑基坑工程监测技术规范[S].2009.
[2]周予启.刘卫未.袁革忠, 等.超大深基坑工程支护设计与施工[J].天津建设科技, 2013 (5) .
[3]钟正雄.基坑工程监测数据库管理系统的设计及应用[J].地下空间, 1998 (S1) .
[4]罗晓辉, 白世伟.深基坑大变形耦合分析与数值模拟[J].岩土力学, 2003 (6) .
关键词：岩土工程 数值模拟有限差分有限元边界 元离散 元无界元
1.引言
近几十年来，随着计算机应用的发展，数值计算方法在岩土工程问题分析中迅速得到了广泛应用，大大推动了岩土（体）力学的发展。在岩土（体）力学中所用的数值方法主要有以下几种：有限差分法、有限元法、边界元法、加权余量法、半解析元法、刚体元法、非连续变形分析法、离散元法、无界元法和流形元法等。下面就对这些方法进行简要的介绍和分析。
2.有限差分法
有限差分法是一种比较古老且应用较广的一种数值方法。它的基本思想是将待解决问题的基本方程和边界条件近似地用差分方程来表示，这样就把求解微分方程的问题转化为求解代数方程的问题。亦即它将实际的物理过程在时间和空间上离散，分解成有限数量的有限差分量，近似假设这些差分量足够小，以致在差分量的变化范围内物体的性能和物理过程都是均匀的，并且可以用来描述物理现象的定律，只是在差分量之间发生阶跃式变化。有限差分法的原理是将实际连续的物理过程离散化，近似地置换成一连串的阶跃过程，用函数在一些特定点的有限差商代替微商，建立与原微分方程相应的差分方程，从而将微分方程转化为一组代数方程，通常采用“显式”时间步进方法来求解代数方程组。
3.有限单元法
有限元法将连续的求解域离散为有限数量单元的组合体，解析地模拟或逼近求解区域。由于单元能按各种不同的联结方式组合在一起，且单元本身又可有不同的几何形状，所以可以适应各种复杂几何形状的求解域。它的原理是利用每个单元内假设的近似函数来表示求解区域上待求的未知场函数，单元内的近似函数由未知场函数在各个单元节点上的数值以及插值函数表达。这就使未知场函数的节点值成为新未知量，把一个连续的无限自由度问题变成离散的有限自由度问题。只要解出节点未知量，便可以确定单元组合体上的场函数，随着单元数目的增加，近似解收敛于精确解。按所选未知量的类型，有限元法可分为位移型、平衡型和混合型有限元法。位移型有限元法在计算机上更易实现，且易推广到非线性和动力效应等方面，故比其他类型的有限元法应用广泛。
4.边界元法
边界元法出现在20 世纪60 年代，是一种求解边值问题的数值方法。它是以Betti 互等定理为基础，有直接法与间接法两种。直接边界元法是以互等定理为基础建立起来的，而间接边界元法是以叠加原理为基础建立起来的。边界元法原理是把边值问题归结为求解边界积分方程的问题，在边界上划分单元，求边界积分方程的数值解，进而求出区域内任意点的场变量，故又称为边界积分方程法。边界元法只需对边界进行离散和积分，与有限元法相比，具有降低维数、输入数据较简单、计算工作量少、精度高等优点。比较适合于在无限域或半无限域问题的求解，尤其是等效均质围岩地下工程问题。边界元法的基本解本身就有奇异性，可以比较方便地处理所谓奇异性问题，故目前边界元法得到研究人员的青睐。
5.加权余量法
加权余量法也是一种求解微分方程的数值法，它在流体力学、热传导以及化学工程等方面应用较广。它具有两个方面的优点：①由于加权余量法是直接从控制方程出发去求解问题，理论简单，不需要复杂的数学处理，且它的应用与问题的能量泛函是否存在无关，因而它的应用范围较广，利用加权余量法这一优点去建立有限单元的刚度矩阵，可以大大扩展有限元法的应用范围；②加权余量法的计算程序简单，要求解的代数方程组阶数较低，对计算机内存容量要求不高，计算所需要的原始数据较少，这样就大大减轻了准备工作量。
6.半解析元法
半解析元法是Y. K. Cheung 于1968 年提出来的，同有限元法一样，它也是基于变分原理的。不同点是半解析元法根据结构的类型和特点，利用部分已有的解析结果，选择一定的位移函数，使解沿某些方向直接引入已知解析函数系列，而不再离散为数值计算点，因此自由度和计算量大大降低。这几年半解析法发展很快，种类很多，主要包括有限条法、有限层法、有限厚条法、有限壳条法、样条有限元法以及无限元法等。这类方法适用于求解高维、无限域及动力场等较复杂的问题。
7.无界元法
无界元法是P. Bettess 于1977 年提出来的，用于解决用有限元法求解无限域问题时，人们常会遇到的“计算范围和边界条件不易确定”的问题，是有限元法的推广。其基本思想是适当地选取形函数和位移函数，使得当局部坐标趋近于1 时，整体坐标趋于无穷大而位移为零，从而满足计算范围无限大和无限远处位移为零的条件。它与有限元法等数值方法耦合对于解决岩土（体）力学问题也是一种有效方法。上述介绍的几种数值法都是针对连续介质的，只能获得某一荷载或边界条件下的稳定解。
8.离散单元法
离散单元法随着非连续岩石力学的发展而不断进步，与现有的连续介质力学方法相比，还有以下问题需要研究：
（1）刚体离散单元法是基于非连续岩石力学的，更适合于低应力状态下具有明显发育构造面的坚硬岩体的变形失稳分析。对于软弱破碎、节理裂隙非常发育和高应力状态下的岩体变形失稳分析，则不适合。
（2）岩体介质种类繁多，性质非常复杂。在通常情况下，节理岩体或颗粒体表现为非均质和各向异性，并且常表现有很强的非线性，所处的地质环境不尽相同，这就使得岩土工程计算有很多不确定性因素。离散元的主要计算参数（如阻尼参数、刚度系数），影响到岩土工程稳定过程的正确模拟以及最终结果的可靠性，尤其是离散元计算中的参数选取，没有统一和完善的确定方法。
（3）计算时步的确定。现在的选取原则是出于满足数学方程趋于收敛的条件，与实际工程问题中的“时间”概念如何联系起来，合理地考虑时间效应，是今后需要研究的问题。
（4）迭代运算的时间较长。用计算机进行离散元计算时，CPU 占用时间较多，特别是在考虑岩块变形的情况下，模型划分单元数受到限制，对迭代方法需做进一步的改进。
9.刚体节理元法
刚体节理元法是Asai 在1981 年提出的，它是在Cundall 刚体离散元间夹有Goodman 节理单元的组合单元，但此节理单元有一定厚度而使离散元间不能“叠合”。刚体节理元法也可考虑不含节理单元的情况，即所谓的单一三角形刚体元非连续变形分析法，是石根华博士和古德曼教授于1984 年首次提出的一种新型数值分析方法，至1988 年该方法已形成了一种较为完整的数值计算方法体系。非连续变形分析方法以严格遵循经典力学规则为基础，是一种平行于有限元法的数值计算方法。
10.流形元方法
学 号：
学 院：
机 电 学 院 日 期：
2015 年 X 月X 日 目 录 实验一 函数插值方法 1 实验二 函数逼近与曲线拟合 5 实验三 数值积分与数值微分 7 实验四 线方程组的直接解法 9 实验五 解线性方程组的迭代法 15 实验六 非线性方程求根 19 实验七 矩阵特征值问题计算 21 实验八 常微分方程初值问题数值解法 24 实验一 函数插值方法 一、问题提出 对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。
数据如下：
（1）0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,）（2）1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，计算的,值。（提示：结果为,）二、要求 1、利用Lagrange插值公式 编写出插值多项式程序；
2、给出插值多项式或分段三次插值多项式的表达式；
3、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何；
4、对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下：
其中：
三、目的和意义 1、学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；
2、明确插值多项式和分段插值多项式各自的优缺点；
3、熟悉插值方法的程序编制；
4、如果绘出插值函数的曲线，观察其光滑性。
四、实验步骤（1）0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,）第一步：先在matlab中定义lagran的M文件为拉格朗日函数代码为：
function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1;for j=1:n+1 if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));end end l(k,:)=v;end c=y*l;end 第二步：然后在matlab命令窗口输入：
>>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382];>> lagran(x,y)回车得到：
ans =121.6264-422.7503 572.5667-377.2549 121.9718-15.0845 由此得出所求拉格朗日多项式为 p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步：在编辑窗口输入如下命令：
>> x=[0.4 0.55 0.65 0.80,0.95 1.05];>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845;>> plot(x,y)命令执行后得到如下图所示图形，然后 >> x=0.596;>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.6262 得到f（0.596）=0.6262 同理得到f（0.99）=1.0547（2）1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，和分段三次插值多项式，计算的,值。（提示：结果为,）实验步骤：
第一步定义 function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1;for j=1:n+1 if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));end end l(k,:)=v;end c=y*l;end 定义完拉格朗日M文件 第二步：然后在matlab命令窗口输入：
>>>> x=[1 2 3 4 5 6 7];y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001];>> lagran(x,y)回车得到：
ans =0.0001-0.0016 0.0186-0.1175 0.4419-0.9683 0.9950 由此得出所求拉格朗日多项式为 p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950 第三步：在编辑窗口输入如下命令：
>> x=[1 2 3 4 5 6 7];>> y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950;>> plot(x,y)命令执行后得到如下图所示图形，然后 >> x=1.8;>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.1650 得到f（0.596）=0.6262 同理得到f（6.15）=2.3644 五、实验结论 插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。
实验二 函数逼近与曲线拟合 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。
t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 二、要求 1、用最小二乘法进行曲线拟合；
2、近似解析表达式为；
3、打印出拟合函数，并打印出与的误差，；
4、另外选取一个近似表达式，尝试拟合效果的比较；
5、* 绘制出曲线拟合图。
三、目的和意义 1、掌握曲线拟合的最小二乘法；
2、最小二乘法亦可用于解超定线代数方程组；
3、探索拟合函数的选择与拟合精度间的关系 四、实验步骤：
第一步先写出线性最小二乘法的M文件 function c=lspoly(x,y,m)n=length(x);b=zeros(1:m+1);f=zeros(n,m+1);for k=1:m+1 f(:,k)=x.^(k-1);end a=f＇*f;b=f＇*y＇;c=a＼b;c=flipud(c);第二步在命令窗口输入：
>>lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)回车得到：
ans =-0.0024 0.2037 0.2305 即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305 在编辑窗口输入如下命令：
>> x=[0,5,10,15,20,25,30,35,40,45,50,55];>> y=-0.0024*x.^2+0.2037*x+0.2305;>> plot(x,y)命令执行得到如下图 五、实验结论 分析复杂实验数据时,常采用分段曲线拟合方法。利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。
实验三 数值积分与数值微分 一、问题提出 选用复合梯形公式，复合Simpson公式，Romberg算法，计算（1）（2）（3）（4）二、要求 1、编制数值积分算法的程序；
2、分别用两种算法计算同一个积分，并比较其结果；
3、分别取不同步长，试比较计算结果（如n = 10, 20等）；
4、给定精度要求ε，试用变步长算法，确定最佳步长。
三、目的和意义 1、深刻认识数值积分法的意义；
2、明确数值积分精度与步长的关系；
3、根据定积分的计算方法，可以考虑二重积分的计算问题。
四、实验步骤 第一步：编写各种积分的程序 复合梯形程序如下：
function I=TX(x,y)n=length(x);m=length(y);if n~=m error(＇The lengths of X and Y must be equal＇);return;end h=(x(n)-x(1))/(n-1);a=[1 2*ones(1,n-2)1];I=h/2*sum(a.*y);复合Simpson程序如下：
function s = simpr1(f,a,b,n)h=(b-a)/(2*n);s1=0;s2=0;for k=1:10 x=a+h*(2*k-1);s1=s1+feval(f,x);end for k=1:(10-1)x=a+h*2*k;s2=s2+feval(f,x);end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;end Romberg程序如下：
function I = Romber_yang(fun,a,b,ep)if nargin<4 ep=1e-5;end;m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while 1 N=2^(m-1);h=h/2;I=I/2;for i=1:N I=I+h*feval(fun,a+(2*i-1)*h);end T(m+1,1)=I;M=2*N;k=1;while M>1;T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;end if abs(T(k,k)-T(k-1,k-1)) 2、对于积分Ι=01sinXXdx，f(0)=1,梯形积分T=0.94607307，辛普森积分S=0.94607308，Romberg积分R=0.94607307。 3、对于积分Ι=01eX4+X2dx，梯形积分T=0.39081248，辛普森积分S=0.39081185，Romberg积分R=0.39081885。 4、对于积分Ι=01ln1+X1+X2dx，梯形积分T=0.27218912，辛普森积分S=0.27219844，Romberg积分R=0.27219827。 五、实验结论，通过本实验学会复合梯形公式，复合Simpson公式，Romberg公式的编程与应用，掌握MATLAB提供的计算积分的各种函数的使用方法。 实验四 线方程组的直接解法 一、问题提出 给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、设线性方程组 2、设对称正定阵系数阵线方程组 3、三对角形线性方程组 二、要求 1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法； 平方根法与改进平方根法； 追赶法求解（选择其一）； 2、应用结构程序设计编出通用程序； 3、比较计算结果，分析数值解误差的原因； 4、尽可能利用相应模块输出系数矩阵的三角分解式。 三、目的和意义 1、通过该课题的实验，体会模块化结构程序设计方法的优点； 2、运用所学的计算方法，解决各类线性方程组的直接算法； 3、提高分析和解决问题的能力，做到学以致用； 4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。 四、实验步骤： 列主元高斯消去法的matlab的M文件程序 function [x,det,index]=Gauss(A,b)% 求线形方程组的列主元Gauss消去法，其中，% A为方程组的系数矩阵； % b为方程组的右端项； % x为方程组的解； % det为系数矩阵A的行列式的值； % index为指标变量，index=0表示计算失败，index=1表示计算成功。 [n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。 if n~=m error(＇The rows and columns of matrix A must be equal!＇);return;end % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m~=nb error(＇The columns of A must be equal the length of b!＇);return;end % 开始计算，先赋初值 index=1;det=1;x=zeros(n,1);for k=1:n-1 % 选主元 a_max=0;for i=k:n if abs(A(i,k))>a_max a_max=abs(A(i,k));r=i;end end if a_max<1e-10 index=0;return;end % 交换两行 if r>k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;end z=b(k);b(k)=b(r);b(r)=z;det=-det;end % 消元过程 for i=k+1:n m=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j);end b(i)=b(i)-m*b(k);end det=det*A(k,k);end det=det*A(n,n);% 回代过程 if abs(A(n,n))<1e-10 index=0;return;end for k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j);end x(k)=b(k)/A(k,k);end 然后在命令窗口输入 >> A=[4 2-3-1 2 1 0 0 0 0;8 6-5-3 6 5 0 1 0 0;4 2-2-1 3 2-1 0 3 1;0-2 1 5-1 3-1 1 9 4;-4 2 6-1 6 7-3 3 2 3;8 6-8 5 7 17 2 6-3 5;0 2-1 3-4 2 5 3 0 1;16 10-11-9 17 34 2-1 2 2;4 6 2-7 13 9 2 0 12 4;0 0-1 8-3-24-8 6 3-1];>> b=[5 12 3 2 3 46 13 38 19-21];>> gauss(A,b)ans = 1.0000-1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000-1.0000 2.0000 高斯-约当消去法maltab的M文件程序 function [x,flag]=Gau_Jor(A,b)% 求线形方程组的列主元Gauss-约当法消去法，其中，% A为方程组的系数矩阵； % b为方程组的右端项； % x为方程组的解； [n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。 if n~=m error(＇The rows and columns of matrix A must be equal!＇);return;end % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m~=nb error(＇The columns of A must be equal the length of b!＇);return;end % 开始计算，先赋初值 flag=＇ok＇;x=zeros(n,1);for k=1:n % 选主元 max1=0;for i=k:n if abs(A(i,k))>max1 max1=abs(A(i,k));r=i;end end if max1<1e-10 falg=＇failure＇;return;end % 交换两行 if r>k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;end z=b(k);b(k)=b(r);b(r)=z;end % 消元过程 b(k)=b(k)/A(k,k);for j=k+1:n A(k,j)=A(k,j)/A(k,k);end for i=1:n if i~=k for j=k+1:n A(i,j)=A(i,j)-A(i,k)*A(k,j);end b(i)=b(i)-A(i,k)*b(k);end end end % 输出x for i=1:n x(i)=b(i);end 然后保存后在命令窗口输入： >> A=[4 2-4 0 2 4 0 0;2 2-1-2 1 3 2 0;-4-1 14 1-8-3 5 6;0-2 1 6-1-4-3 3;2 1-8-1 22 4-10-3;4 3-3-4 4 11 1-4;0 2 5-3-10 1 14 2;0 0 6 3-3-4 2 19];>> b=[0-6 20 23 9-22-15 45];>> Gau_Jor(A,b)ans = 121.1481-140.1127 29.7515-60.1528 10.9120-26.7963 5.4259-2.0185 五、实验结论 用LU法，调用matlab中的函数lu中，L往往不是一个下三角，但可以直接计算不用它的结果来计算，不用进行行变换。如果进行行变b也要变，这样会很麻烦。 实验五 解线性方程组的迭代法 一、问题提出 对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。 二、要求 1、体会迭代法求解线性方程组，并能与消去法做以比较； 2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢； 3、对方程组2，3使用SOR方法时，选取松弛因子ω=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者； 4、给出各种算法的设计程序和计算结果。 三、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较； 2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序； 3、体会上机计算时，终止步骤或k >（给予的迭代次数），对迭代法敛散性的意义； 4、体会初始解，松弛因子的选取，对计算结果的影响。 四、实验步骤 第一步编写实验所需的Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序。 Jacobi迭代法： function [x,k,index]=J(A,b,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 for i=1:n y(i)=b(i);for j=1:n if j~=i y(i)=y(i)-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end y(i)=y(i)/A(i,i);end if norm(y-x,inf) function [x,k,index]=G(A,b,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 y=x;for i=1:n z=b(i);for j=1:n if j~=i z=z-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end z=z/A(i,i);x(i)=z;end if norm(y-x,inf) function [x,k,index]=SOR(A,b,ep,w,itmax)if nargin<5 itmax=100;end if nargin<4 w=1;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 y=x;for i=1:n z=b(i);for j=1:n if j~=i z=z-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end z=z/A(i,i);x(i)=(1-w)*x(i)+w*z;end if norm(y-x,inf) 1、设线性方程组 2、设对称正定阵系数阵线方程组 3、三对角形线性方程组 五、实验结论 迭代法是解线性方程组的一个重要的实用方法，特别适用于求解在实际中大量出现的，系数矩阵为稀疏阵的大型线性方程组。通过此次实验学会了Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序编写，并掌握了它们各自的优缺点及其适用条件。 实验六 非线性方程求根 一、问题提出 设方程有三个实根 现采用下面六种不同计算格式，求 f(x)=0的根或 1、2、3、4、5、6、二、要求 1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况； 2、用事后误差估计来控制迭代次数，并且打印出迭代的次数； 3、初始值的选取对迭代收敛有何影响； 4、分析迭代收敛和发散的原因。 三、目的和意义 1、通过实验进一步了解方程求根的算法； 2、认识选择计算格式的重要性； 3、掌握迭代算法和精度控制； 4、明确迭代收敛性与初值选取的关系。 四、实验步骤 第一步：编写实验所需的程序。 function [x_star,index,it]=DD(fun,x,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end index=0;k=1;while k 1、，x1=0，x2=0。 2、，x1=无穷大，x2=-0.3473。 3、，x1=1.8794，x2=1.8794。 4、，x1=-0.3473，x2=-0.3473.。 5、，x1=1.8794，x2=1.8794。 6、，x1=1.8794，x2=-0.3473。 五、实验结论 对于非线性方程，求它的解析解有时候是很困难的，但采用数值方法可以很容易地求它的近似解。此次实验就是采用迭代法求非线性方程的根。对于一个非线性方程，选用不同的迭代形式，因为其收敛程度不一样，造成其效率与精确度有很大的差别。 实验七 矩阵特征值问题计算 一、问题提出 利用冪法或反冪法，求方阵的按模最大或按模最小特征值及其对应的特征向量。 设矩阵A的特征分布为： 且 试求下列矩阵之一（1）求，及 取 结果（2）求及 取 结果： （3）求及 取 结果(4)取 这是一个收敛很慢的例子,迭代次才达到 结果（5）有一个近似特征值,试用幂法求对应的特征向量，并改进特征值（原点平移法）。 取 结果 二、要求 1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计； 2、会用原点平移法改进算法，加速收敛； 对矩阵B=A-PI取不同的P值，试求其效果； 3、试取不同的初始向量，观察对结果的影响； 4、对矩阵特征值的其它分布，如且如何计算。 三、目的和意义 1、求矩阵的部分特征值问题具有重要实际意义，如求矩阵谱半径，稳定性问题往往归于求矩阵按模最小特征值； 2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧； 3、问题中的题（5），反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征向量。 四、实验步骤 第一步：写出实验所需的幂法求最大特征值及反幂法求最小特征值的程序。 幂法程序： function [m,u,index]=TZ(A,ep,itmax)if nargin<3 itmax=100;end if nargin<2 ep=1e-5;end n=length(A);u=ones(n,1);index=0;k=0;m1=0;while k<=itmax v=A*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;if abs(m-m1) function [m,u,index]=FTZ(A,ep,itmax)if nargin<3 itmax=100;end if nargin<2 ep=1e-5;end n=length(A);u=ones(n,1);index=0;k=0;m1=0;invA=inv(A);while k<=itmax v=invA*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;if abs(m-m1) λ3=3.4723，x3=（1.0000 0.5229 0.2422）T。，λ1=21.3053，X1=（0.8724 0.5401 0.9973 0.5644 0.4972 1.0000）T； λ6=1.6214。 五、实验结论 求n阶方阵A的特征值和特征向量，也是实际中常常碰到的问题。通过此次实验掌握了用幂法和反幂法求一个方阵的最大特征值和特征向量，绝对值最小的特征值和特征向量。 实验八 常微分方程初值问题数值解法 一、问题提出 科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题： （1）分别取h=0.1,0.2,0.4时数值解。 初值问题的精确解。 （2）用r=3的Adams显式和预-校式求解 取步长h=0.1，用四阶标准R-K方法求值。 （3）用改进Euler法或四阶标准R-K方法求解 取步长0.01,计算数值解，参考结果。 （4）利用四阶标准R-K方法求二阶方程初值问题的数值解（I）(II)(III)(IV) 二、要求 1、根据初值问题数值算法，分别选择二个初值问题编程计算； 2、试分别取不同步长，考察某节点处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。 三、目的和意义 1、熟悉各种初值问题的算法，编出算法程序； 2、明确各种算法的精度与所选步长有密切关系； 3、通过计算更加了解各种算法的优越性。 四、实验步骤 function [x,y]=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;end h=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n;x(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end 实验程序及分析（Ⅰ）（1）、算法程序 function E =Euler_1(fun,x0,y0,xN,N)% Euler向前公式，其中 % fun为一阶微分方程的函数 % x0,y0为初始条件 % xN为取值范围的一个端点 % h为区间步长 % N为区间个数 % x为Xn构成的向量 % y为yn构成的向量 x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));end T=[x＇,y＇] function z=f(x,y)z=4*x/y-x*y;(2)、运行程序 >> Euler_1(＇f＇,0,3,2,20)结果 ： 一、实验3.1 题目： 考虑线性方程组，，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss消去过程。 （1）取矩阵，则方程有解。取计算矩阵的条件数。分别用顺序Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解，结果如何？ （2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。 （3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。 1.算法介绍 首先，分析各种算法消去过程的计算公式，顺序高斯消去法： 第k步消去中，设增广矩阵中的元素（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k行以下各行计算，分别用乘以增广矩阵的第行并加到第行，则可将增广矩阵中第列中以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即； 列主元高斯消去法： 第k步消去中，在增广矩阵中的子方阵中，选取使得，当时，对中第行与第行交换，然后按照和顺序消去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即； 完全主元高斯消去法： 第k步消去中，在增广矩阵中对应的子方阵中，选取使得，若或，则对中第行与第行、第列与第列交换，然后按照和顺序消去法相同的步骤进行即可。重复此方法，从第1步进行到第n-1步，就可以得到最终的增广矩阵，即； 接下来，分析回代过程求解的公式，容易看出，对上述任一种消元法，均有以下计算公式： 2.实验程序的设计 一、输入实验要求及初始条件； 二、计算系数矩阵A的条件数及方程组的理论解； 三、对各不同方法编程计算，并输出最终计算结果。 3.计算结果及分析 （1） 先计算系数矩阵的条件数，结果如下，可知系数矩阵的条件数较大，故此问题属于病态问题，b或A的扰动都可能引起解的较大误差； 采用顺序高斯消去法，计算结果为： 最终解为x=(1.***,1.***,1.***,1.***,0.***,1.***,0.***,1.***,0.***,1.***)T 使用无穷范数衡量误差，得到=2.842***1e-14，可以发现，采用顺序高斯消元法求得的解与精确解之间误差较小。通过进一步观察，可以发现，按照顺序高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此顺序高斯消去法方式并没有对结果造成特别大的影响。 若采用列主元高斯消元法，则结果为： 最终解为x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T 同样使用无穷范数衡量误差，有=0； 若使用完全主元高斯消元法，则结果为 最终解x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T 同样使用无穷范数衡量误差，有=0； （2） 若每步都选取模最小或尽可能小的元素为主元，则计算结果为 最终解x=(1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.***)T 使用无穷范数衡量误差，有为2.842***1e-14；而完全主元消去法的误差为=0。 从（1）和（2）的实验结果可以发现，列主元消去法和完全主元消去法都得到了精确解，而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了一定的影响，但由于方程组的维度较低，并且元素之间相差不大，所以误差仍比较小。 为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下： 第n次消元 顺序 列主元 完全主元 模最小 6.*** 6.*** 4.*** 4.*** 4.*** 4.*** 4.***3333 4.***3333 4.*** 4.*** 4.*** 4.*** 4.0***063 4.0***063 4.*** 4.*** 4.0039*** 4.0039*** 4.*** 0.0***469 0.0***469 4.*** 从上表可以发现，对这个方程组而言，顺序高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小，所以计算过程中的累积误差也较小，最终4种方法的输出结果均较为精确。 在这里，具体解释一下顺序法与模最小法的计算结果完全一致的原因。该矩阵在消元过程中，每次选取主元的一列只有两个非零元素，对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。在这样的情况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一致的。 理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特殊性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。 （3） 时，重复上述实验过程，各种方法的计算结果如下所示，在这里，仍采用无穷范数衡量绝对误差。 顺序高斯消去法 列主元高斯消去 完全主元高斯消去 选取模最小或尽可能小元素作为主元消去 X 1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 2.***e-11 0 0 2.***e-11 可以看出，此时列主元和完全主元的计算结果仍为精确值，而顺序高斯消去和模尽可能小方法仍然产生了一定的误差，并且两者的误差一致。与n=10时候的误差比相比，n=20时的误差增长了大约1000倍，这是由于计算过程中舍入误差的不断累积所致。所以，如果进一步增加矩阵的维数，应该可以看出更明显的现象。 （4） 不同矩阵维度下的误差如下，在这里，为方便起见，选取2-条件数对不同维度的系数矩阵进行比较。 维度 条件数 顺序消去 列主元 完全主元 模尽量小 1.7e+3 2.84e-14 0 0 2.84e-14 1.8e+6 2.91e-11 0 0 2.91e-11 5.7e+7 9.31e-10 0 0 9.31e-10 1.8e+9 2.98e-08 0 0 2.98e-08 1.9e+12 3.05e-05 0 0 3.05e-05 3.8e+16 3.28e+04 3.88e-12 3.88e-12 3.28e+04 8.5e+16 3.52e+13 4.2e-3 4.2e-3 3.52e+13 从上表可以看出，随着维度的增加，不同方法对计算误差的影响逐渐体现，并且增长较快，这是由于舍入误差逐步累计而造成的。不过，方法二与方法三在维度小于40的情况下都得到了精确解，这两种方法的累计误差远比方法一和方法四慢；同样地，出于与前面相同的原因，方法一与方法四的计算结果保持一致，方法二与方法三的计算结果保持一致。 4.结论 本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此，不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显，四种方法都足够精确。 对比四种方法，可以发现采用列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确。不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。 另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发现列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。 附录：程序代码 clear clc; format long; %方法选择 n=input(＇矩阵A阶数:n=＇); disp(＇选取求解方式＇); disp(＇1 顺序Gauss消元法，2 列主元Gauss消元法，3 完全选主元Gauss消元法，4 模最小或近可能小的元素作为主元＇); a=input(＇求解方式序号：＇); %赋值A和b A=zeros(n,n); b=zeros(n,1); for i=1:n A(i,i)=6; if i>1 A(i,i-1)=8; end if i A(i,i+1)=1; end end for i=1:n for j=1:n b(i)=b(i)+A(i,j); end end disp(＇给定系数矩阵为：＇); A disp(＇右端向量为：＇); b %求条件数及理论解 disp(＇线性方程组的精确解：＇); X=(A＼b)＇ fprintf(＇矩阵A的1-条件数: %f ＼n＇,cond(A,1)); fprintf(＇矩阵A的2-条件数: %f ＼n＇,cond(A)); fprintf(＇矩阵A的无穷-条件数: %f ＼n＇,cond(A,inf)); %顺序Gauss消元法 if a==1 A1=A;b1=b; for k=1:n if A1(k,k)==0 disp(＇主元为零，顺序Gauss消元法无法进行＇); break end fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A1(k,k)) %disp(＇此次消元后系数矩阵为：＇); %A1 for p=k+1:n l=A1(p,k)/A1(k,k); A1(p,k:n)=A1(p,k:n)-l*A1(k,k:n); b1(p)=b1(p)-l*b1(k); end end x1(n)=b1(n)/A1(n,n); for k=n-1:-1:1 for w=k+1:n b1(k)=b1(k)-A1(k,w)*x1(w); end x1(k)=b1(k)/A1(k,k); end disp(＇顺序Gauss消元法解为:＇); disp(x1); disp(＇所求解与精确解之差的无穷-范数为＇); norm(x1-X,inf) end %列主元Gauss消元法 if a==2 A2=A;b2=b; for k=1:n [max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k)))); if max_i~=k A2_change=A2(k,:); A2(k,:)=A2(max_i,:); A2(max_i,:)=A2_change; b2_change=b2(k); b2(k)=b2(max_i); b2(max_i)=b2_change; end if A2(k,k)==0 disp(＇主元为零，列主元Gauss消元法无法进行＇); break end fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A2(k,k)) %disp(＇此次消元后系数矩阵为：＇); %A2 for p=k+1:n l=A2(p,k)/A2(k,k); A2(p,k:n)=A2(p,k:n)-l*A2(k,k:n); b2(p)=b2(p)-l*b2(k); end end x2(n)=b2(n)/A2(n,n); for k=n-1:-1:1 for w=k+1:n b2(k)=b2(k)-A2(k,w)*x2(w); end x2(k)=b2(k)/A2(k,k); end disp(＇列主元Gauss消元法解为:＇); disp(x2); disp(＇所求解与精确解之差的无穷-范数为＇); norm(x2-X,inf) end %完全选主元Gauss消元法 if a==3 A3=A;b3=b; for k=1:n VV=eye(n); [max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n))))); if numel(max_i)==0 [max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n))))); end W=eye(n); W(max_i(1)+k-1,max_i(1)+k-1)=0; W(k,k)=0; W(max_i(1)+k-1,k)=1; W(k,max_i(1)+k-1)=1; V=eye(n); V(k,k)=0; V(max_j(1)+k-1,max_j(1)+k-1)=0; V(k,max_j(1)+k-1)=1; V(max_j(1)+k-1,k)=1; A3=W*A3*V; b3=W*b3; VV=VV*V; if A3(k,k)==0 disp(＇主元为零，完全选主元Gauss消元法无法进行＇); break end fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A3(k,k)) %disp(＇此次消元后系数矩阵为：＇); %A3 for p=k+1:n l=A3(p,k)/A3(k,k); A3(p,k:n)=A3(p,k:n)-l*A3(k,k:n); b3(p)=b3(p)-l*b3(k); end end x3(n)=b3(n)/A3(n,n); for k=n-1:-1:1 for w=k+1:n b3(k)=b3(k)-A3(k,w)*x3(w); end x3(k)=b3(k)/A3(k,k); end disp(＇完全选主元Gauss消元法解为:＇); disp(x3); disp(＇所求解与精确解之差的无穷-范数为＇); norm(x3-X,inf) end %模最小或近可能小的元素作为主元 if a==4 A4=A;b4=b; for k=1:n AA=A4; AA(AA==0)=NaN; [min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k)))); if numel(min_i)==0 [min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n)))); end W=eye(n); W(min_i(1)+k-1,min_i(1)+k-1)=0; W(k,k)=0; W(min_i(1)+k-1,k)=1; W(k,min_i(1)+k-1)=1; A4=W*A4; b4=W*b4; if A4(k,k)==0 disp(＇主元为零，模最小Gauss消元法无法进行＇); break end fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A4(k,k)) %A4 for p=k+1:n l=A4(p,k)/A4(k,k); A4(p,k:n)=A4(p,k:n)-l*A4(k,k:n); b4(p)=b4(p)-l*b4(k); end end x4(n)=b4(n)/A4(n,n); for k=n-1:-1:1 for w=k+1:n b4(k)=b4(k)-A4(k,w)*x4(w); end x4(k)=b4(k)/A4(k,k); end disp(＇模最小Gauss消元法解为:＇); disp(x4); disp(＇所求解与精确解之差的无穷-范数为＇); norm(x4-X,inf) end 二、实验3.3 题目： 考虑方程组的解，其中系数矩阵H为Hilbert矩阵： 这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端的办法给出确定的问题。 （1）选择问题的维数为6，分别用Gauss消去法（即LU分解）、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何。 （2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何？计算的结果说明的什么？ （3）讨论病态问题求解的算法。 1.算法设计 对任意线性方程组，分析各种方法的计算公式如下，（1）Gauss消去法： 首先对系数矩阵进行LU分解，有，则原方程转化为，令，则原方程可以分为两步回代求解： 具体方法这里不再赘述。 （2）J迭代法： 首先分解，再构造迭代矩阵，其中，进行迭代计算，直到误差满足要求。 （3）GS迭代法： 首先分解，再构造迭代矩阵，其中，进行迭代计算，直到误差满足要求。 （4）SOR迭代法： 首先分解，再构造迭代矩阵，其中，进行迭代计算，直到误差满足要求。 2.实验过程 一、根据维度n确定矩阵H的各个元素和b的各个分量值； 二、选择计算方法（Gauss消去法，J迭代法，GS迭代法，SOR迭代法），对迭代法设定初值，此外SOR方法还需要设定松弛因子； 三、进行计算，直至满足误差要求（对迭代法，设定相邻两次迭代结果之差的无穷范数小于0.0001； 3.计算结果及分析 （1）时，问题可以具体定义为 计算结果如下，Gauss消去法 第1次消元所选取的主元是：1 第2次消元所选取的主元是：0.0833333 第3次消元所选取的主元是：0.00555556 第4次消元所选取的主元是：0.000357143 第5次消元所选取的主元是：2.26757e-05 第6次消元所选取的主元是：1.43155e-06 解得X=(0.*** 1.*** 0.*** 1.*** 0.*** 1.***)T 使用无穷范数衡量误差，可得=4.***e-10； J迭代法 设定迭代初值为零，计算得到 J法的迭代矩阵B的谱半径为4.30853＞1，所以J法不收敛； GS迭代法 设定迭代初值为零，计算得到GS法的迭代矩阵G的谱半径为：0.999998＜1，故GS法收敛，经过541次迭代计算后，结果为X=(1.001***6 0.*** 0.*** 1.*** 1.*** 0.***)T 使用无穷范数衡量误差，有=0.***； SOR迭代法 设定迭代初值为零向量，并设定，计算得到SOR法迭代矩阵谱半径为0.***，经过100次迭代后的计算结果为 X=(1.*** 0.*** 1.03***59 1.06***81 1.*** 0.9***527)T； 使用无穷范数衡量误差，有=0.***； 对SOR方法，可变，改变值，计算结果可以列表如下 迭代次数 迭代矩阵的谱半径 0.*** 0.*** 0.*** 0.*** X 1.*** 0.*** 1.01***40 1.*** 1.0***681 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.*** 1.*** 0.*** 0.*** 1.05***66 0.*** 1.*** 0.*** 1.*** 0.*** 0.*** 0.*** 0.*** 0.*** 可以发现，松弛因子的取值对迭代速度造成了不同的影响，上述四种方法中，松弛因子=0.5时，收敛相对较快。 综上，四种算法的结果列表如下： 算法 Gauss消去法 Jacobi法 GS法 SOR法（取） 迭代次数 -- 不收敛 541 迭代矩阵的谱半径 -- 4.30853 0.999998 0.*** X 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** -- 1.001***6 0.*** 0.*** 1.*** 1.*** 0.*** 1.*** 0.*** 1.03***59 1.06***81 1.*** 0.9***527 4.***e-10 -- 0.*** 0.*** 计算可得，矩阵H的条件数为>>1，所以这是一个病态问题。由上表可以看出，四种方法的求解都存在一定的误差。下面分析误差的来源： LU分解方法的误差存在主要是由于Hilbert矩阵各元素由分数形式转换为小数形式时，不能除尽情况下会出现舍入误差，在进行LU分解时也存在这个问题，所以最后得到的结果不是方程的精确解，但结果显示该方法的误差非常小； Jacobi迭代矩阵的谱半径为4.30853，故此迭代法不收敛； GS迭代法在迭代次数为541次时得到了方程的近似解，其误差约为0.05，比较大。GS迭代矩阵的谱半径为0.999998，很接近1，所以GS迭代法收敛速度较慢； SOR迭代法在迭代次数为100次时误差约为0.08，误差较大。SOR迭代矩阵的谱半径为0.999999，也很接近1，所以时SOR迭代法收敛速度不是很快，但是相比于GS法，在迭代速度方面已经有了明显的提高；另外，对不同的，SOR方法的迭代速度会相应有变化，如果选用最佳松弛因子，可以实现更快的收敛； （2） 考虑不同维度的情况，时，算法 Gauss消去 J法 GS法 SOR法（w=0.5） 计算结果 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** -- 0.*** 1.*** 0.*** 1.*** 1.*** 1.*** 0.9968*** 0.*** 1.*** 0.9397*** 0.*** 1.*** 1.*** 1.*** 0.*** 0.*** 迭代次数 -- -- 356 谱半径 -- 6.04213 0.*** -- 时，算法 Gauss消去法 Jacobi法 GS法 SOR法（w=0.5） 计算结果 0.*** 1.*** 0.*** 1.000***1 0.*** 1.*** 0.*** 1.*** 0.*** 1.*** 0.*** -- 0.*** 1.*** 0.*** 0.*** 0.*** 1.02***91 1.*** 1.*** 1.*** 0.*** 0.947***7 1.0***572 0.*** 0.*** 0.*** 1.*** 1.*** 1.*** 1.*** 0.*** 0.*** 0.*** 迭代次数 -- -- 1019 谱半径 -- 8.64964 0.*** -- 时 算法 Gauss消去法 Jacobi法 GS法 SOR法（w=0.5） 计算结果 0.*** 1.*** 0.*** 0.*** 1.*** 0.*** 2.*** -2.*** 7.*** -7.*** 7.*** -1.*** 0.*** 1.*** 0.*** -- 不收敛 1.*** 1.*** 0.907***9 0.*** 0.*** 1.*** 1.09***64 1.*** 1.*** 1.*** 1.0385*** 0.*** 0.942***3 0.*** 0.*** 迭代次数 -- -- 262 谱半径 -- 6.04213 ＞1 1.*** 8.*** -- -- 0.*** 分析以上结果可以发现，随着n值的增加，Gauss消去法误差逐渐增大，而且误差增大的速度很快，在维数小于等于10情况下，Gauss消去法得到的结果误差较小；但当维数达到15时，计算结果误差已经达到精确解的很多倍； J法迭代不收敛，无论n如何取值，其谱半径始终大于1，因而J法不收敛，所以J迭代法不能用于Hilbert矩阵的求解； 对于GS迭代法和SOR迭代法，两种方法均收敛，GS迭代法是SOR迭代法松弛因子取值为1的特例，SOR方法受到取值的影响，会有不同的收敛情况。可以得出GS迭代矩阵的谱半径小于1但是很接近1，收敛速度很慢。虽然随着维数的增大，所需迭代的次数逐渐减少，但是当维数达到15的时候，GS法已经不再收敛。因此可以得出结论，GS迭代方法在Hilbert矩阵维数较低时，能够在一定程度上满足迭代求解的需求，不过迭代的速度很慢。另外，随着矩阵维数的增加，SOR法的误差水平基本稳定，而且误差在可以接受的范围之内。 经过比较可以得出结论，如果求解较低维度的Hibert矩阵问题，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的结果精确度较高；如果需要求解较高维度的Hibert矩阵问题，只有采用SOR迭代法。 （3） 系数矩阵的条件数较大时，为病态方程。由实验可知，Gauss法在解上述方程时，结果存在很大的误差。而对于收敛的迭代法，可以通过选取最优松弛因子的方法来求解，虽然迭代次数相对较多，但是结果较为精确。 总体来看，对于一般病态方程组的求解，可以采用以下方式： 1.低维度下采用Gauss消去法直接求解是可行的； Jacobi迭代方法不适宜于求解病态问题； GS迭代方法可以解决维数较低的病态问题，但其谱半径非常趋近于1，导致迭代算法收敛速度很慢，维数较大的时候，GS法也不再收敛； SOR方法较适合于求解病态问题，特别是矩阵维数较高的时候，其优势更为明显。 2.采用高精度的运算，如选用双倍或更多倍字长的运算，可以提高收敛速度； 3.可以对原方程组作某些预处理，从而有效降低系数矩阵的条件数。 4.实验结论 （1）对Hibert矩阵问题，其条件数会随着维度的增加迅速增加，病态性会越来越明显；在维度较低的时候，Gauss消去法、GS迭代法和SOR迭代法均可使用，且可以优先使用Gauss消去法；如果需要求解较高维度的Hibert矩阵问题，只有SOR迭代法能够求解。 （2）SOR方法比较适合于求解病态问题，特别是矩阵维数较高的时候，其优点更为明显。从本次实验可以看出，随着矩阵维数的增大，SOR方法所需的迭代次数减少，而且误差基本稳定，是解决病态问题的适宜方法。 附录：程序代码 clear all clc; format long; %矩阵赋值 n=input(＇矩阵H的阶数:n=＇); for i=1:n for j=1:n H(i,j)=1/(i+j-1); end end b=H*ones(n,1); disp(＇H矩阵为：＇); H disp(＇向量b：＇); b %方法选择 disp(＇选取求解方式＇); disp(＇1 Gauss消去法，2 J迭代法，3 GS迭代法，4 SOR迭代法＇); a=input(＇求解方式序号：＇); %Gauss消去法 if a==1; H1=H;b1=b; for k=1:n if H1(k,k)==0 disp(＇主元为零，Gauss消去法无法进行＇); break end fprintf(＇第%d次消元所选取的主元是：%g＼n＇,k,H1(k,k)) for p=k+1:n m5=-H1(p,k)/H1(k,k); H1(p,k:n)=H1(p,k:n)+m5*H1(k,k:n); b1(p)=b1(p)+m5*b1(k); end end x1(n)=b1(n)/H1(n,n); for k=n-1:-1:1 for v=k+1:n b1(k)=b1(k)-H1(k,v)*x1(v); end x1(k)=b1(k)/H1(k,k); end disp(＇Gauss消去法解为:＇); disp(x1); disp(＇解与精确解之差的无穷范数＇); norm((x1-a),inf) end D=diag(diag(H)); L=-tril(H,-1); U=-triu(H,1); %J迭代法 if a==2; %给定初始x0 ini=input(＇初始值设定:x0=＇); x0(:,1)=ini*diag(ones(n)); disp(＇初始解向量为:＇); x0 xj(:,1)=x0(:,1); B=(D^(-1))*(L+U); f=(D^(-1))*b; fprintf(＇(J法B矩阵谱半径为：%g＼n＇,vrho(B)); if vrho(B)<1; for m2=1:5000 xj(:,m2+1)=B*xj(:,m2)+fj; if norm((xj(:,m2+1)-xj(:,m2)),inf)<0.0001 break end end disp(＇J法计算结果为:＇); xj(:,m2+1) disp(＇解与精确解之差的无穷范数＇); norm((xj(:,m2+1)-diag(ones(n))),inf) disp(＇J迭代法迭代次数:＇); m2 else disp(＇由于B矩阵谱半径大于1，因而J法不收敛＇); end end %GS迭代法 if a==3; %给定初始x0 ini=input(＇初始值设定:x0=＇); x0(:,1)=ini*diag(ones(n)); disp(＇初始解向量为:＇); x0 xG(:,1)=x0(:,1); G=inv(D-L)*U; fG=inv(D-L)*b; fprintf(＇GS法G矩阵谱半径为：%g＼n＇,vrho(G)); if vrho(G)<1 for m3=1:5000 xG(:,m3+1)=G*xG(:,m3)+fG; if norm((xG(:,m3+1)-xG(:,m3)),inf)<0.0001 break; end end disp(＇GS迭代法计算结果:＇); xG(:,m3+1) disp(＇解与精确解之差的无穷范数＇); norm((xG(:,m3+1)-diag(ones(n))),inf) disp(＇GS迭代法迭代次数:＇); m3 else disp(＇由于G矩阵谱半径大于1，因而GS法不收敛＇); end end %SOR迭代法 if a==4; %给定初始x0 ini=input(＇初始值设定:x0=＇); x0(:,1)=ini*diag(ones(n)); disp(＇初始解向量为:＇); x0 A=H; for i=1:n b(i)=sum(A(i,:)); end x_star=ones(n,1); format long w=input(＇松弛因子:w=＇); Lw=inv(D-w*L)*((1-w)*D+w*U); f=w*inv(D-w*L)*b; disp(＇迭代矩阵的谱半径：＇) p=vrho(Lw) time_max=100;%迭代次数 x=zeros(n,1);%迭代初值 for i=1:time_max x=Lw*x+f; end disp(＇SOR迭代法得到的解为＇); x disp(＇解与精确解之差的无穷范数＇); norm((x_star-x),inf) end pause 三、实验4.1 题目： 对牛顿法和拟牛顿法。进行非线性方程组的数值求解 （1）用上述两种方法，分别计算下面的两个例子。在达到精度相同的前提下，比较其迭代次数、CPU时间等。 （2）取其他初值，结果又如何？反复选取不同的初值，比较其结果。 （3）总结归纳你的实验结果，试说明各种方法适用的问题。 1.算法设计 对需要求解的非线性方程组而言，牛顿法和拟牛顿法的迭代公式如下，（1）牛顿法： 牛顿法为单步迭代法，需要取一个初值。 （2）拟牛顿法：（Broyden秩1法） 其中，拟牛顿法不需要求解的导数，因此节省了大量的运算时间，但需要给定矩阵的初值，取为。 2.实验过程 一、输入初值； 二、根据误差要求，按公式进行迭代计算； 三、输出数据； 3.计算结果及分析 （1）首先求解方程组（1），在这里，设定精度要求为，方法 牛顿法 拟牛顿法 初始值 计算结果X x1 0.*** 0.*** x2 1.*** 1.0852*** x3 0.*** 0.*** 迭代次数 CPU计算时间/s 3.777815 2.739349 可以看出，在初始值相同情况下，牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果基本相同，但牛顿法的迭代次数明显要少一些，但是，由于每次迭代都需要求解矩阵的逆，所以牛顿法每次迭代的CPU计算时间更长。 之后求解方程组（2），同样设定精度要求为 方法 牛顿法 拟牛顿法 初始值 计算结果X x1 0.*** 0.*** x2 0.*** 0.*** x3 -0.*** -0.*** 迭代次数 CPU计算时间/s 2.722437 3.920195 同样地，可以看出，在初始值相同情况下，牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果是基本相同的，但牛顿法的迭代次数明显要少，但同样的，由于每次迭代中有求解矩阵的逆的运算，牛顿法每次迭代的CPU计算时间较长。 （2）对方程组（1），取其他初值，计算结果列表如下，同样设定精度要求为 初始值 方法 牛顿法 拟牛顿法 计算结果 0.*** 1.*** 0.*** 9.21***94 -5.*** 18.1***205 迭代次数 CPU计算时间/s 3.907164 4.818019 计算结果 0.*** 1.*** 0.*** 9.21***91 -5.*** 18.1***807 迭代次数 2735 CPU计算时间/s 8.127286 5.626023 计算结果 0.*** 1.*** 0.*** 0.*** 1.0852*** 0.*** 迭代次数 CPU计算时间/s 3.777815 2.739349 计算结果 0.*** 1.*** 0.*** 0.*** 1.*** 0.*** 迭代次数 188 CPU计算时间/s 3.835697 2.879070 计算结果 9.21***22 -5.*** 18.1***605 Matlab警告矩阵接近奇异值，程序进入长期循环计算中 迭代次数 -- CPU计算时间/s 4.033868 -- 计算结果 0.*** 1.*** 0.*** Matlab警告矩阵接近奇异值，程序进入长期循环计算中 迭代次数 -- CPU计算时间/s 12.243263 -- 从上表可以发现，方程组（1）存在另一个在(9.2,-5.6,18.1)T附近的不动点，初值的选取会直接影响到牛顿法和拟牛顿法最后的收敛点。 总的来说，设定的初值离不动点越远，需要的迭代次数越多，因而初始值的选取非常重要，合适的初值可以更快地收敛，如果初始值偏离精确解较远，会出现迭代次数增加直至无法收敛的情况； 由于拟牛顿法是一种近似方法，拟牛顿法需要的的迭代次数明显更多，而且收敛情况不如牛顿法好（初值不够接近时，甚至会出现奇异矩阵的情况），但由于牛顿法的求解比较复杂，计算时间较长； 同样的，对方程组（2），取其他初值，计算结果列表如下，同样设定精度要求为 初始值 方法 牛顿法 拟牛顿法 计算结果 0.*** 0.*** -0.*** 0.*** 0.*** -0.*** 迭代次数 CPU计算时间/s 2.722437 3.920195 计算结果 0.*** 0.*** -0.*** 0.*** -0.*** 76.*** 迭代次数 CPU计算时间/s 5.047111 5.619752 计算结果 0.*** 0.*** -0.*** 1.0e+02 * -0.*** -0.000***6 1.754***3 迭代次数 CPU计算时间/s 3.540668 3.387829 计算结果 0.*** 0.*** -0.*** 1.0e+04 * 0.*** -0.*** 1.*** 迭代次数 CPU计算时间/s 2.200571 2.640901 计算结果 0.*** 0.*** -0.*** 矩阵为奇异值，无法输出准确结果 迭代次数 -- CPU计算时间/s 1.719072 -- 计算结果 0.*** 0.*** -0.*** 矩阵为奇异值，无法输出准确结果 迭代次数 149 -- CPU计算时间/s 2.797116 -- 计算结果 矩阵为奇异值，无法输出准确结果 矩阵为奇异值，无法输出准确结果 迭代次数 -- -- CPU计算时间/s -- -- 在这里，与前文类似的发现不再赘述。 从这里看出，牛顿法可以在更大的区间上实现压缩映射原理，可以在更大的范围上选取初值并最终收敛到精确解附近； 在初始值较接近于不动点时，牛顿法和拟牛顿法计算所得到的结果是基本相同的，虽然迭代次数有所差别，但计算总的所需时间相近。 （3） 牛顿法在迭代过程中用到了矩阵的求逆，其迭代收敛的充分条件是迭代满足区间上的映内性，对于矩阵的求逆过程比较简单，所以在较大区间内满足映内性的问题适合应用牛顿法进行计算。一般而言，对于函数单调或者具有单值特性的函数适合应用牛顿法，其对初始值敏感程度较低，算法具有很好的收敛性。 另外，需要说明的是，每次计算给出的CPU时间与计算机当时的运行状态有关，同时，不同代码的运行时间也不一定一致，所以这个数据并不具有很大的参考价值。 4.实验结论 对牛顿法和拟牛顿法，都存在初始值越接近精确解，所需的迭代次数越小的现象； 在应用上，牛顿法和拟牛顿法各有优势。就迭代次数来说，牛顿法由于更加精确，所需的迭代次数更少；但就单次迭代来说，牛顿法由于计算步骤更多，且计算更加复杂，因而每次迭代所需的时间更长，而拟牛顿法由于采用了简化的近似公式，其每次迭代更加迅速。当非线性方程组求逆过程比较简单时，如方程组1的情况时，拟牛顿法不具有明显的优势；而当非线性方程组求逆过程比较复杂时，如方程组2的情况，拟牛顿法就可以体现出优势，虽然循环次数有所增加，但是CPU耗时反而更少。 另外，就方程组压缩映射区间来说，一般而言，对于在区间内函数呈现单调或者具有单值特性的函数适合应用牛顿法，其对初始值敏感程度较低，使算法具有很好的收敛性；而拟牛顿法由于不需要在迭代过程中对矩阵求逆，而是利用差商替代了对矩阵的求导，所以即使初始误差较大时，其倒数矩阵与差商偏差也较小，所以对初始值的敏感程度较小。 附录：程序代码 %方程1，牛顿法 tic; format long; %%初值 disp(＇请输入初值＇); a=input(＇第1个分量为：＇); b=input(＇第2个分量为：＇); c=input(＇第3个分量为：＇); disp(＇所选定初值为＇); x=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; while e>E F=[12*x(1)-x(2)^2-4*x(3)-7;x(1)^2+10*x(2)-x(3)-11;x(2)^3+10*x(3)-8]; f=[12,-2*x(2),-4;2*x(1),10,-1;0,3*x(2)^2,10]; det_x=((f)^(-1))*(-F); x=x+det_x; e=max(norm(det_x)); i=i+1; end disp(＇迭代次数＇); i disp(＇迭代次数＇); x toc; %方程1，拟牛顿法 tic; format long; %%初值 %%初值 disp(＇请输入初值＇); a=input(＇第1个分量为：＇); b=input(＇第2个分量为：＇); c=input(＇第3个分量为：＇); disp(＇所选定初值为＇); x0=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; A0=eye(3); while e>E F0=[12*x0(1)-x0(2)^2-4*x0(3)-7;x0(1)^2+10*x0(2)-x0(3)-11;x0(2)^3+10*x0(3)-8]; x1=x0-A0^(-1)*F0; s=x1-x0; F1=[12*x1(1)-x1(2)^2-4*x1(3)-7;x1(1)^2+10*x1(2)-x1(3)-11;x1(2)^3+10*x1(3)-8]; y=F1-F0; A1=A0+(y-A0*s)*s＇/(s＇*s); x0=x1; A0=A1; e=max(norm(s)); i=i+1; end disp(＇迭代次数＇); i disp(＇迭代次数＇); x0 toc; %方程2，牛顿法 tic; format long; %%初值 disp(＇请输入初值＇); a=input(＇第1个分量为：＇); b=input(＇第2个分量为：＇); c=input(＇第3个分量为：＇); disp(＇所选定初值为＇); x=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; while e>E F=[3*x(1)-cos(x(2)*x(3))-0.5;x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;exp(1)^(-x(1)*x(2))+20*x(3)+(10*pi-3)/3]; f=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3));2*x(1),-162*x(2)-81/5,cos(x(3));-x(2)*exp(1)^(-x(1)*x(2)),-x(1)*exp(1)^(-x(1)*x(2)),20]; det_x=((f)^(-1))*(-F); x=x+det_x; e=max(norm(det_x)); i=i+1; end disp(＇迭代次数＇); i disp(＇迭代次数＇); x toc; %方程2，拟牛顿法 tic; format long; %%初值 %%初值 disp(＇请输入初值＇); a=input(＇第1个分量为：＇); b=input(＇第2个分量为：＇); c=input(＇第3个分量为：＇); disp(＇所选定初值为＇); x0=[a;b;c] %%误差要求 E=0.0001; %%迭代 i=0; e=2*E; A0=eye(3); while e>E F0=[3*x0(1)-cos(x0(2)*x0(3))-0.5;x0(1)^2-81*(x0(2)+0.1)^2+sin(x0(3))+1.06;exp(1)^(-x0(1)*x0(2))+20*x0(3)+(10*pi-3)/3]; x1=x0-A0^(-1)*F0; s=x1-x0; F1=[3*x1(1)-cos(x1(2)*x1(3))-0.5;x1(1)^2-81*(x1(2)+0.1)^2+sin(x1(3))+1.06;exp(1)^(-x1(1)*x1(2))+20*x1(3)+(10*pi-3)/3]; y=F1-F0; A1=A0+(y-A0*s)*s＇/(s＇*s); x0=x1; A0=A1; e=max(norm(s)); i=i+1; end disp(＇迭代次数＇); i disp(＇迭代次数＇); x0 1．应提交一份完整的实习报告。具体要求如下： （1）要有封面，封面上要标明姓名、学号、专业和联系电话； （2）要有序言，说明所用语言及简要优、特点，说明选用的考量； （3）要有目录，指明题目、程序、计算结果，图标和分析等内容所在位置，作到 信息简明而完全； （4）要有总结，全方位总结机编程计算的心得体会； （5）尽量使报告清晰明了，一般可将计算结果、图表及对比分析放在前面，程序 清单作为附录放在后面，程序中关键部分要有中文说明或标注，指明该部分的功能和作用。 2．程序需完好保存到期末考试后的一个星期，以便老师索取用于验证、询问或质疑部分内容。 3．认真完成实验内容，可以达到既学习计算方法又提高计算能力的目的，还可以切身体会书本内容之精妙所在，期间可以得到很多乐趣。 4．拷贝或抄袭他人结果是不良行为，将视为不合格。 5．报告打印后按要求的时间提交给任课老师。 数值分析上机试题 10（选择其中两个题目） 1． 给定三个n阶线性方程组Ax=b,其中A的元素aij(i,j,=1,…,n)与阶数分别为 （1）aij=(i+j-1)2, n=3,4,5,…,9; （2）aij=(i+j)2, n=3,4,5 （3）aij=1/(i+j-1), n=3,4,5,6; b的元素bi=ai1+ai2+…,ain,(i=1,2,…,n) 已知其准确解为x=(1,1,...,1)T, 用列主元素法分别求解上列方程组。输出各步主元，解释所遇到的现象。 2．用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。 (1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]； b1=[-3,2,4]T, b2=[100,-200,345]T,(2)A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1] T, b2=[5,0,-10]T,1 (3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T,3.松弛因子对SOR法收敛速度的影响。 用SOR法求解方程组Ax=b,其中 41-3141-2-2...,B ....141-2-314 要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2,...,1.9进行迭代，记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w0或w2会有什么影响？ 4.某实际问题中两次过程记录测得(x,y)数据点如下：其中yi=f(xi),i=1,…,21 第一次的数据为： x=[-5.0000-4.5000-4.0000-3.5000-3.0000-2.5000-2.0000 -1.5000-1.0000-0.500000.50001.00001.5000 2.00002.50003.00003.50004.00004.50005.0000]; y =[-0.0000-0.0001-0.0002-0.0003-0.0004-0.0048-0.0366 -0.1581-0.3679-0.389400.38940.36790.1581 0.03660.00480.00040.00030.00020.00010.0000] 第二次的数据为： x=[-5.0000-4.7000-4.4000-4.1000-3.8000-3.5000-3.2000 -2.9000-2.6000-2.3000-2.0000-1.7000-1.4000-1.1000 -0.8000-0.5000-0.20000.10000.40000.70001.0000 1.30001.60001.90002.20002.50002.80003.1000 3.40003.70004.00004.30004.60004.9000] y=[-0.0000-0.0000-0.0000-0.0000-0.0000-0.0000-0.0001 -0.0006-0.0030-0.0116-0.0366-0.0945-0.1972-0.3280 -0.4218-0.3894-0.19220.09900.34090.42880.3679 0.23990.12370.05140.01740.00480.00110.0002 0.00000.00000.00000.00000.00000.0000] 现考虑分段低次或样条插值，或直接插值，总之用不超过7次的分段或不分段插值多项式近似该过程的对应规律y=f(x)，请给出一种解决方案使得误差尽量小表示尽量简单。 5.用Runge-Kutta 4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。 1.1背景 在海洋平台的建造过程中,由于吊架比塔架节省工期,大跳因其抗弯刚度大,变形小,故应用广泛。 1.2有限元模型建立 1.2.1模型基本参数及模型单元选取 钢管直径d=48mm,厚度θ=3.5mm,钢板密度ρ=7.85x103Kg/m3,重力加速度g=9.8m/s2,钢材泊松比μ=0.3,弹性模量E=2.06x105MPa,屈服强度σy=235MPa,屈服后硬化强度E=2.06x103MP。脚手架钢管采用8节点等参壳单元。 1.2.2边界条件与加载 立杆和大跳的连接通过刚性连接协调变形的方法,采用分布线荷载的方法进行加载,有限元模型如图1所示。 2理论验证 将有限元计算结果与理论[2]计算结果比较,见表1。得出有限元模拟与规范算法之间的误差在1.2%~9.3%之间,满足需要。 3大跳数值模拟分析 已知:吊架吊点间隔4m,采用6m大跳连接,立杆直接挂在6m大跳上,间距1.4m,这种搭设方案缺乏理论依据,本节通过有限元模拟对其进行分析。 强度分析: 3.1有效应力分析。大跳有效应力分布如图2所示,有效应力最大值132.0MPa,位置出现在大跳小立杆和上部水平杆交界处,距大跳左边缘0.75m。 3.2弯曲应力分析。大跳弯曲应力分布如图3所示,最大值106.6MPa,位置出现在大跳小立杆和上部的水平杆交界处,距大跳左边缘0.75m。 3.3挠度分析。大跳竖向挠度分布如图4示,挠度最大值5.61mm,位置出现在大跳下部水平杆中部。 经分析,大跳的有效应力及弯曲应力最大值均小于允许应力156.7 MPa,挠度最大值为5.61mm,小于《建筑施工扣件式脚手架安全技术规程》规定的挠度(1400/150=9.33mm),可见大跳的强度满足要求。因此,通过6m大跳来提供4m间隔的这种搭设形式,在每个工作面荷载为312Kg/m2时,安全性能符合规范要求,这种吊架搭设方案可推广。 4结论 4.1 ADINA有限元软件分析计算的脚手架横向水平杆、纵向水平强度和挠度与规范算法之间的误差在1.2%~9.3%之间,误差比较小。 4.2通过6m大跳来提供4m间隔的这种吊架搭设形式,在每个工作面荷载为312Kg/m2时,安全性能符合规范,此搭设方案可推广使用。 摘要：本文选取以海洋平台钢结构上吊架式脚手架的6m大跳在吊架中新的搭设形式为切入点,应用ADINA有限元软件,考虑扣件的半刚性连接及材料的弹塑性力学性能,采用逐步加载的方法,建立了大跳的有限元模型,对吊架中应用6m大跳来连接4m间隔搭设方案的安全性能进行评价。 关键词：吊架,大跳,有限元分析 参考文献 [1]《建筑施工扣件式脚手架安全技术规程》JGJ130-2001[S].2001.13~22. 【关键词】城市隧道；数值模拟；地表变形；CRD工法 0.引言 目前，随着岩土工程数值方法和计算机技术的快速发展，复杂定解条件问题的处理才能成为可能，数值方法被人们广泛用来进行求解隧道施工过程中遇到的各种问题的最有效的通用方法。 目前对于CRD法在城市浅埋隧道施工中的运用及其对地表变形和衬砌受力的内在因素和机理方面的研究较少。本文结合龙岩至厦门高速铁路（简称龙厦铁路）石桥头隧道进口段工程， 对CRD法施工进行数值模拟，重点研究隧道施工引起的地表变形沉降规律、衬砌受力特点、衬砌受力与围岩变形的关系，最后通过现场监测结果验证数值模拟的合理性。 1.工程概况 龙厦高速铁路石桥头隧道位于福建省龙岩市城区，隧道整体由西北折向东南方向，为电气化双线隧道，线间距4~4.4m，行车速度120km/h，位于R=1000m的右偏曲线上，隧道净高11.2m，净宽12.8m，全长1586m，是龙厦铁路重点控制性工程。地质调查和钻探揭示，进口段DK2+450~DK3+021长571m，为Ⅴ级围岩。且该段隧道埋深浅，开挖跨度大，围岩稳定条件差，因此，控制围岩变形及地表沉降是该段隧道施工中的关键技术问题，因此本文以此段为依据进行数值模拟分析，来研究隧道施工引起的地表变形规律。 2.隧道施工数值模型的建立 由于进口段地质条件较差，选取隧道埋深约17m建立计算模型。隧道净空断面为（12.8m×11.2m），根据隧道开挖的影响范围，参考已有的计算理论和经验，模型计算范围为：中心线两侧各取50m，竖向选取地表以下70m，为了节约时间，隧道纵向取三个开挖步长即9m。在本文分析计算中，假定隧道围岩按弹塑性材料考虑，破坏准则采用摩尔-库伦准则。结合该工程，在模拟隧道施工过程中，超前管棚预加固，采用改善围岩参数的等效方法来模拟；钢拱架加喷射混凝土初期支护，采用壳体结构单元模拟；围岩、加固圈和二次衬砌，采用实体单元模拟。整个地层数值模型计算区域为100m×9m×70m，共划分网格单元13620个，节点23360个。 模型计算边界条件为：左右边界施加水平方向的约束，底部边界施加固定约束，顶部边界为自由面。模拟计算中，土层力学参数和围岩支护参数，按石桥头隧道地质勘查报告和规范选取，其值见表1。 表1模型參数表 3.CRD工法施工过程的数值模拟 建立三维网格模型→初始条件模拟分析(自重应力和顶面的均布荷载) →加固圈加固处理(超前管棚加固) →开挖左上导洞模拟计算、初期支护模拟计算→开挖右上导洞模拟计算、初期支护模拟计算→开挖左下导洞模拟计算、初期支护模拟计算→开挖右下导洞模拟计算、初期支护模拟计算→二次衬砌施加模拟计算→最后对模拟计算结果进行分析。 注意的是在开挖前必须采取一些加固措施(如超前管棚等)进行预加固地层，每部开挖之后要及时施作初期支护，最后拆除临时支护并施作二次衬砌。 通过对CRD工法在城市浅埋隧道施工进行三维数值模拟，得到各部开挖支护和二次支护的竖向位移及影响范围分布(结合tecplot10.3后期处理软件得图1)、围岩应力分布(结合tecplot10.3后期处理软件得图2)、衬砌受力分布云图(如图3)、开挖结束后的变形矢量及塑性状态分布(如图4、图5)。 第一部开挖支护第二部开挖支护 第三部开挖支护 第四部开挖支护 二次支护 图1各部开挖支护和二次支护竖向位移及影响范围 表2模拟开挖地表下沉统计/mm 由图3和表2模拟计算结果表明： 3.1采用CRD法对隧道进行施工，隧道上方出现了变形沉降，沉降影响范围随着开挖位置变化而变化，沉降影响范围首先产生并偏向于开挖一则，随着开挖的继续，影响范围也随之变化，最后影响区域几乎对称位于隧道中心线上方两侧各26m左右，最大沉降值为16.2mm，位于隧道拱顶正上方。因此在对地表沉降进行监控量测时，基准点要埋设远离隧道中心线26m以外。 3.2隧道各部开挖进行支护后，沉降虽然继续增加，但数值不大，并且主要影响范围明显缩小，由此知道在隧道开挖过程中，隧道毛洞开挖结束时，初期支护和临时支护要及时跟进，对稳定隧道开挖面极为重要。 3.3隧道不同部位的开挖对地表沉降的影响各异，左上导洞开挖沉降比例较大，因土体开挖打破了原有土体三向应力平衡，地应力第一次将重新分布，必然引起周围土体的位移和变形，由于第一步开挖非常关键，在进行开挖前有必要采取相应措施（如超前管棚加固），并及时配合临时支护，尽量控制因开挖引起的变形位移。上部导洞开挖沉降量占总沉降量多达70%，由于隧道上部是主要支撑拱顶的关键部分，而这部分土体被开挖，上部围岩此时处于临空状态，内部应力大量释放，使得拱部土体受到了较大的拉应力作用，其上部土体岩层变形较大，导致地表沉降量也较大，因此在该阶段施工过程中支护要及时，并且拱顶部位钢拱架喷混凝土的固定及连接要严密，超前小导管注浆及时跟进，以保证开挖掌子面的整体稳定。随后进行的下部开挖，地表沉降进一步增大，但相对上部而言，该部位开挖对地表沉降影响明显较小，但在施工过程中也要严格按照设计方案施工。最后二次衬砌的沉降量占总沉降4%左右，在隧道施工过程中，要求初次衬砌结束后，并且在隧道地表沉降和洞内收敛基本稳定时，再进行二次衬砌施工。 图2a围岩水平应力分布图图2b围岩竖向应力分布图 图3a衬砌竖向应力云图图3b衬砌水平应力云图 图4 变形矢量图图5 塑性状态图 由图2到图5模拟计算结果表明： 3.4衬砌周围土层应力主要为压应力，最大压应力分布在隧道顶部和底部区域，且顶部土层压应力分布与地表沉降类似为抛物线，仅在隧道两侧的局部区域有较小的拉应力。 3.5在隧道周围土体和地表荷载作用下，衬砌应力以压应力为主，局部为拉应力。在垂直方向上受力最大值达4.06MPa，主要位于隧道拱腰处，且应力比较集中，因此在施工过程中，拱腰处的钢支撑接头需重点处理，做好锁脚加强支护，防止应力集中对衬砌造成变形开裂，降低衬砌结构的整体刚度，从而影响整个围岩的稳定。隧道两腰上半部分区域承受拉应力，最大值约为1.28MPa，拉应力的产生与地层应力分布、隧道断面形状及拱顶变形等诸多因素有关，且有扩大趋势，并向腰部，顶部集中。拱顶部位受到水平应力集中影响，最大值达0.72MPa，这主要是由于拱顶部位在开挖过程中产生较大沉降有关，且该部位在开挖时作为主要受力部位承受垂直和水平两个方向的地应力。因此施工时，要保证拱顶处的固定连接，接头的连续性会对衬砌变形和结构整体沉降产生较大影响，再配合必要的加固措施以保证开挖面的整体稳定。 3.6隧道顶部受力复杂，应力集中，顶部以上区域变形较大，是隧道开挖造成地面沉降主要区域；在隧道坑壁及向外延伸线上，主要受压缩剪切，是剪切破坏的主要区域，如果长期不支护，可能会出现流动现象，发生剪切破坏，造成隧道失稳而影响地面变形过大。 4.变形监控量测 对于城市浅埋隧道施工而言，监控量测具有重要地位，它是信息化施工的基础，是确保隧道施工和周围环境安全的重要手段，是检验设计参数、地面稳定性、评价施工方法的主要依据。对于石桥头隧道工程，地表沉降使用高可靠性磁阻尼摆式补偿(补偿范围±15、、)的索佳SDL30电子水准仪，配合RAB码玻璃钢水准尺(标准差±1.0mm)进行量测；由于隧道上方多为民居低层建筑物，且建筑物周围空地较窄，因而采用通过测量建筑物基础相对沉陷的方法来确定建筑物的倾斜，因此量测仪器与地表沉降测量仪器相同，而对于一些高层建筑物，则采用可实现高精度，远距离无协作目标测距，反射片测距及棱镜测距的索佳SET230RK全站仪。由图6看出，隧道地表最终沉降的实测值与模拟值相比，误差在1.8~3.9mm之间，对于相同里程，拱顶最终沉降量一般比地表最终沉降量稍大，但最终沉降量都处于规范允许之內。地表房屋得到安全保障，隧道顺利贯通，工程工期明显缩短，不仅为之后铺轨作业提供了充裕时间，且为业主单位节省了工程经费，达到预期经济和社会效益。 图6实测地表和拱顶沉降曲线 5.结论 ５．１本文所建立的三维有限差分模型，可有效的分析CRD法在城市浅埋隧道施工引起的地表沉降规律和衬砌受力特点。 ５．２采用CRD法对城市浅埋隧道施工，沉降影响范围随着开挖位置变化而变化，沉降影响范围首先产生并偏向于开挖一则，随着开挖的继续，影响范围也随之变化，最后影响范围几乎对称的位于隧道中心线上方。 ５．３采用CRD法对城市浅埋隧道施工，初期支护对地表沉降影响范围有明显控制效果。上部导坑开挖对整个地表变形起到关键作用，沉降量占总沉降量的70%左右。 ５．４衬砌周围的土层应力主要为压应力，且顶部土层应力的分布与地表沉降类似，为抛物线，仅在隧道两侧的局部区域有较小的拉应力。而衬砌应力同样以压应力为主，局部区域为拉应力，高应力区集中在顶部和腰部。 ５．５隧道顶部受力复杂，应力集中，顶部以上区域变形较大，是隧道开挖造成地表沉降主要区域。在隧道坑壁及向外延伸线上，主要受压缩剪切，是剪切破坏的主要区域，如果长期不支护，可能会出现流动现象，发生剪切破坏，造成隧道失稳而影响地面变形过大。 【参考文献】 ［１］严绍洋，李亮辉，张春宇.公路隧道开挖渗流场的有限差分法分析[J].中外公路，2007，27：120-123． ［２］杨玉胜.红砂岩隧道开挖稳定性的三维数值分析[J].路基工程，2007，67-69．清华大学数值分析实验报告 篇5
西南交大数值分析上机实习报告 篇6
土木工程数值分析 篇7
城市隧道施工的三维数值模拟分析 篇8
