地震数值模拟实验报告（共7篇）
地震数值模拟实验报告 篇1
地震模拟实验报告
一、实验目的
感受地震灾害发生时带来的感受，了解地震灾害发生时的现象，明白地震的成因、带来的危害及防灾自救措施。
二、实验内容
观看介绍地震相关知识的影片，到地震屋感受模拟地震，观看地震知识展板，聆听老师讲解地震知识。
三、实验原理阐述
（一）地震的分类及其成因 根据地震的成因，可分为构造地震、火山地震、塌陷地震、人工地震、诱发地震五类。其中，构造地震、火山地震和塌陷地震是由自然原因引起的地震，称为天然地震。而人工地震和诱发地震则是由人类活动引起的地震，称为人为地震。
1、构造地震
构造地震是由于地下深处岩层张裂、扩大而造成的地震。这类地震发生的次数最多，破坏力也最大，约占全世界地震总数的90%以上。在众多构造地震形成机制的解释中，断层说是最受认可、最重要的一种解释，它较为合理科学地解释了构造地震的形成机制。
2、火山地震
火山地震是由于火山作用，如岩浆活动、气体爆炸等引起的地震。火山地震在火山喷发前后都会发生。在火山喷发前，地下岩浆的冲击或热力膨胀的作用会使岩层断裂造成地震。在火山喷发后，由于大量岩浆喷出地表，地壳内部压力减小，造成岩层断裂错动而发生地震。火山地震只可能发生在火山活动区。这类地震发生次数不多，仅占全世界地震总数的7%左右。
3、塌陷地震
塌陷地震是由地下岩洞或矿井顶部塌陷而引起的地震。这类地震的规模比较小，发生的次数也很少。如果发生，往往发生在溶洞密布的石灰岩地区或大规模地下开采的矿区。
4、人工地震
人工地震是指由地下核爆炸、炸药爆破等人为活动引起的地面震动。人工地震是由人为活动引起的地震，如工业爆破、地下核爆炸造成的震动；在深井中进行高压注水以及大水库蓄水后增加了地壳的压力，有时也会诱发地震。
5、诱发地震
诱发地震是由于水库蓄水、油田注水等人类活动引发的地震。这类地震仅发生在特定的水库库区或油田地区。
图一
构造地震
图二
火山、塌陷、人工、诱发地震
（二）地震波
地震波是一种由地震震源发出并在地球内部传播的机械波。地震波可分为纵波、横波和面波三种。纵波又称P波，是一种推进波，传播速度在三种地震波中最快，它使地面发生上下震动，破坏性较弱。横波又称S波，是一种剪切波，它的传播速度慢于纵波，快于面波，它使地面发生前后、左右抖动，破坏性较强。面波又称L波，是由纵波与横波在地表相遇后激发产生的混合波，它只能沿地表面传播，是造成建筑物强烈破坏的主要因素。
（三）地震的分布
世界上的地震主要发生在地震带地区，主要包括环太平洋地震带、地中海—喜马拉雅地震带、大洋中脊地震带和东非裂谷地震带。其中，以环太平洋地震带、地中海—喜马拉雅地震带和东非裂谷地震带造成的危害较大。而大洋中脊地震带因位于大洋深处，且地震活动性较弱，对人类造成的危害较小。
我国位于环太平洋地震带和地中海—喜马拉雅地震带之间，是地震较多的国家之一。我国的地震活动主要分布在五个地区的二十三条地震带上，这五个地区分别是：台湾地区、西南地区、西北地区、华北地区和东南沿海地区。
图三
世界主要地震带分布图
图四
中国地震带分布图
（四）地震的危害
地震发生时往往会伴随着一系列的危害，并导致人员伤亡、建筑物崩塌、交通破坏、农田毁坏等后果。由地震引发的次生灾害包括：火灾、海啸、瘟疫、滑坡和崩塌、水灾、地面塌陷等。
（五）震后自救
首先，无论是在面对突如其来、无法避免的地震袭击时，还是在震后被埋在废墟中时，都要保持镇静。地震来时，迅速寻找路线撤离，如若来不及撤离，则寻找狭小的、牢固的、有利于避震的地区躲避。地震后，如果被困在废墟中，则应保持冷静，设法改变自己的不利处境。在不利处境中，首要的是保护呼吸顺畅，如果闻到有毒气体，就用湿衣服捂住口鼻。在确保自身生命安全的前提下，可以通过用石块敲击能发出响声的物体来向外发出求救信号。如果救援人员长时间未到，则要设法维持自己的生命，尽可能寻找水和食品。
四、实验过程
实验过程可主要分成三个步骤：第一步骤是观看介绍地震成因、现象、危害、自救措施及监测预报的影片；第二步骤是到科技馆可模拟地震的地震屋体验七级地震；第三步骤就是观看地震知识展板，并聆听老师补充的地震的相关知识。
在实验的第二个步骤中，我们进入地震屋亲身感受了七级地震的震感。地震屋是一个装饰得像一个平常的小房间一样的小屋。在屋内，有电视、壁画、风铃、防盗窗和小门。此外，在屋子的四个角落还放置着四个软座，提供给体验者座位；屋子中间摆放着一张固定的圆桌，以给体验者提供攀拉的地方。我们八人一组有序地进入地震屋坐好，双手都紧紧地攀拉着圆桌的钢制外环。在语音提示过后，纵波首先袭来，我们明显地感受到屋子在上下颠簸，自己所坐的位置有上下震动的感觉。然后是横波带来的震感。我们先是感到屋 2
子在缓缓地向右移动，自己的身体也有向右滑动的倾向；然而紧接着，屋子突然加速左移，并猛地有激烈碰撞带来的强烈震动感，我们不得不加大力量握住桌子的外环。在感受横波阶段，其中的一个同学将一本书放到圆桌上，观察书是否会被震落。不一会儿，书就被甩到了地上。从这一现象中，我们可以想象得到地震中书架翻倒、杯子等生活物品坠落的混乱情景。在体验地震的过程中，还伴随着由风铃模拟出的玻璃碎裂的声音，仿佛真的置身于地震现场。
图五
体验地震
图六
地震屋1
图七
地震屋2
五、实验结果及分析
在实验的第一步骤中，我们通过观看影片目睹了地震带给人类的巨大灾难与伤害，我们的心灵受到了深深的震撼，无不感叹于自然力量的强大和人类力量的微弱。人类在自然灾难面前是如此渺小和不堪一击。
在实验的第二步骤中，我们亲身体验了七级地震带来的感受。或许是意识中有了体验地震的心理准备，而不是遭遇突袭的真正地震，我觉得七级的地震并不如想象中的来得强烈和恐怖。但是，我知道，当真正遭受七级地震时，崩塌的房屋、碎裂的玻璃、无法逃避的死亡感，都会让我们的心灵受到无比的震惊与恐惧，甚至丧失了逃生意识。
在实验的第三步骤，通过老师的讲解和地震知识展板上的系统介绍，我们对地震的成因、世界地震带的分布、地震的震级、预测地震的仪器——候风地动仪、地震带来的危害和地震发生后的自救知识都有了较为详细的了解。
图八
候风地震仪
图九
震后废墟
参考文献：
[1]实验教学指导书——走进科学
（一）.[2]陈颙,史培军.自然灾害[M].北京:北京师范大学出版社,2007.
地震数值模拟实验报告 篇2
我国抗震设计规范规定在6、7、8度区, 混凝土小砌块结构分别可以建7、6、5层, 并且开间较小.适量注芯的砌块结构更经济, 节能效果更好, 但是国内外对适量注芯的大开间结构研究很少。且这些成果主要集中在试验研究阶段, 并且都采用被动抗震方式.应用ANSYS有限元软件对八层大开间、注芯率为19%的砌块模型进行模拟, 并对其隔震前后进行模态分析, 并与模型结构试验数据进行对比[1]。
1三维有限元模型
1.1基本假定
为了使模型的反应更接近真实结构, 计算模型采用三维有限元模型。同时为了减少计算量, 作了如下基本假设[2,3]:
(1) 砌体是一种各向同性材料;
(2) 楼板在其自身平面内无限刚;
(3) 模型底板与振动台的连接为刚性连接。
1.2模型单元
模型采用大型通用软件ANSYS来建模, 砌体采shell单元, 钢筋采用link8单元, 底板采用mass21元加在结构底层的最下边的节点上, 而夹层橡胶用combin14单元来模拟, 此单元为带有阻尼的弹单元[4]。
2结构在地震波激励下的地震反应
模型有两个开间:1.515m+1.515m, 进深:2.5m, 层数:8, 层高:0.7m, 总高度:5.6m, 4个橡胶垫的性能参数见表1。其中γ为夹层橡胶垫的剪切变形。按照图1建立隔震模型三维立体有限元模型。输入Northridge和ElCentro两条地震波, 通过调整峰值以满足7度设防对应的小震、中震和大震要求, 沿X方向输入地震波。设计输入峰值的大小分别为El-0.062g, 0.082g, 0.25g和N-0.05g, 0.13g, 0.33g。隔震前后对比地震加速度有限元结果记录见表2。
从表2可看出, 隔震结构上部0层、4层、8层峰值加速度均比台面加速度降低, 上部结构类似于整体平动。屋顶水平加速度反应, 隔震结构只相当于非隔震结构的0.36～0.62, 夹层橡胶垫隔震结构对水平地震作用的隔震效果很明显, 各工况下结构总水平地震作用降低至1/4～1/12, 相当于降低地震烈度约2度。
由图2、图3可以看出, 无论是隔震结构, 还是传统抗震结构, 层间剪力的最大值都集中在第一层和第二层上。因此在进行隔震结构设计时, 同样应该注意结构底下两层的设计。且隔震结构的层间剪力小于非隔震结构, 因此相同工况下进行结构设计时, 可以大大的减少隔震结构的梁柱的截面尺寸和配筋量, 从而可以取得很大的经济效益。
3 结论
(1) 夹层橡胶垫对水平地震作用的隔震效果很明显, 能够对隔振结构降低烈度设计。
(2) 可减小隔震结构的梁柱的截面尺寸和配筋量, 取得很大的经济效益。
参考文献
[1]周抚生, 郭迅, 郑志华, 等.大开间小型混凝土砌块十层模型房屋抗震性能试验研究 (Ⅰ) .地震工程与工程振动, 2002;22 (6) :58—64
[2]张敏政.地震模拟试验中相似率应用的若干问题.地震工程与工程震动, 1997;17 (2) :68—71
[3]GB50011-2001, 建筑抗震设计规范.北京:中国建筑工业出版社, 2001
地震数值模拟实验报告 篇3
关键词 电磁搅拌 耦合场 实验研究
中图分类号：TG27 文献标识码：A
0引言
利用电磁搅拌铝合金半固态浆料的过程为例，进行分析物理量耦合的数字模型，采用有限差分和有限元结合，来进行分析电磁搅拌下温度场、流场和电磁场的耦合模拟，研究流场和温度场参数受电磁搅拌工艺影响的规律，并实验进行验证。根据实验结果显示：铝合金熔体在受水平旋转的电磁力场影响下均匀分布，交变电磁场的旋转方向与其方向一致，大小从边部到中心则因集肤效应递减；电磁力与熔体内速度场分布相似，但随着延长搅拌时间，熔体温度和速度不断降低；熔体凝固速度在电磁搅拌条件下加快，熔体温度梯度在心部和边部的变小。
1电磁搅拌过程模型的建立
1.1 电磁搅拌器结构及模型
电磁搅拌装置中搅拌器磁轭和铁芯系统使用高磁导率的硅钢片制成，搅拌线圈采用在铁芯上安装中空的铜导线作为绕组，一对极的电磁搅拌器制备形成。通过电容使线圈的两对线圈相位角相差90€啊F涔ぷ髟砘救缦拢喊严辔唤浅？0€暗慕涣鞯缌魍ㄈ氲礁杏ο呷χ校桓鼋槐湫懦【突嵩诮涟枸巅鲋胁彼懈盥梁辖鹑厶迨备杏ο呷χ谢岵杏Φ缌鳌4懦『驮亓髀梁辖鹑厶逑嗷プ饔茫梁辖鹑厶迳喜叫绱帕Γ诘绱帕η侣梁辖鹑厶蹇荚硕6杂诟呔侗群艽蟮慕涟枸巅觯朔奖慵扑悖梢杂枚侍饫醇蚧淼绱沤涟枳爸眉扑愕恼瞿P汀5绱沤涟柚平爸糜上呷Α③巅觥⒋砰睢⒖掌敖鹗羧厶宓炔糠肿槌桑治龅闹氐阄梁辖鹑厶濉？
1.2 数字模型的建立
为便于数字模型的建立，假设流动的铝合金熔体是不可压缩的黏性流体。电磁场的搅拌强度在试验过程中并不是很大，相对比较小的熔体流动速度，雷诺数小于2100，因此可以采用层流模型进行建模。除熔体的黏度随温度升降而发生改变，假设均为各向同性材料的不锈钢坩埚、磁轭、铝合金熔体和铜线圈，磁导率均为常数；温度变化对密度的影响不作考虑，除焦耳热外无其他内热源，即加热过程的热膨胀被忽略，密度一直被设定为固定值；产生的位移电流忽略不计；不考虑电磁场受铝合金熔体流动时的影响，与其它场的耦合通过时均电磁力代替时变电磁力来实现；本研究不考虑固液两相流动，而把液态金属熔体的流动状态作为主要考虑的问题。
2 模拟结果和分析
2.1磁场分布分析
在频率为10 Hz、电流为30A的搅拌电流下，计算得到铝合金熔体中电磁力分布和磁力线。不难得出，铝合金熔体被线圈产生的磁力线完全穿过。在铝合金熔体内均匀分布着水平旋转的电磁力场，磁场的旋转方向与电磁力的作用方向一致，使铝合金熔体沿着水平方向上向上旋转流动。由于感应电流的集肤效应，电磁力在合金熔体中由外向内依次减小。说明铝合金熔体均在搅拌区域内受到电磁力的作用，因此熔体的整个区域都能获得均匀的搅拌处理。
2.2流场分布分析
把电流频率设定为10Hz，电流为10 A，计算铝合金熔体内流场在不同时间下的分布，可以看出，在5s的电磁搅拌时间下，还处于完全液态的铝合金熔体，受到旋转电磁力的作用，在液态铝水平面上呈旋转运动，基本实现速度场均匀分布。由于集肤效应，由外向内合金熔体中的电磁力依次减小，这也导致在铝合金熔体边缘流速达到最大值，流速从边缘到中心逐渐减小，搅拌漩涡在铝合金熔体中逐渐形成。随着延长搅拌的时间，金属熔体的流动速度减小，黏度增加，温度降低。搅拌时间为10秒与5s时熔体的动速度明显减小。
2.3温度场分布分析
在10 Hz的电流频率下， 把电流分别设定为10A、20A和30A，求出态温度场的瞬時分布。选取中心到边部路径上的点对金属熔体进行分析，在搅拌时间相同的条件下，随着电流的增大熔体整体的温度反而依次减小。电磁力随着电流的增大而增大，也增强了铝合金熔体的旋转运动，热量传输在熔体内部更快。电流为10A时熔体中心的温度为633.305℃，而电流为30A时熔体中心的温度反而降到617.161℃。因此，搅拌强度的提高能最大限度地增强熔体热交换的动力学过程，要获得好的搅拌效果必须控制好这些参数。
3 结束语
随着不断发展的强磁场技术，有关强磁场的数值计算模型需要建立。使人们对金属材料电磁加工过程可以通过计算模型有深刻地认识和理解，包括电磁搅拌下熔体的熔体内溶质的分布状况、凝固过程、凝固坯细晶组织的机理形成以及非金属物运动轨迹受电磁场影响情况，从而能够更加合理有效的利用电磁场搅拌。总之，在应用范围和计算精度上数字模型应该进一步扩展，使其在实际生产中能真正应用。
参考文献
[1] 张衬新.电磁搅拌作用下铝合金凝固组织的数值模拟[D].北京有色金属研究总院硕士，2012.
[2] 汤孟欧.环缝式电磁搅拌理论与工艺研究[D].北京有色金属研究总院博士，2011.
地震数值模拟实验报告 篇4
新结构传爆药柱多点起爆数值模拟及实验研究
根据冲击波汇聚技术原理及有效装药理论,设计出环锥形传爆药装药结构.用ANSYS/LY-DYNA软件分析起爆点个数对环锥形传爆药柱输出波形的影响,并用多点同步起爆网络起爆环锥形传爆药柱进行实验验证.结果表明,利用多点同步起爆网络起爆环锥形传爆药柱能有效提高传爆药柱的.输出威力,起爆点个数对爆轰波形有很大的影响.
作 者：张金勇 胡双启 曹雄 ZHANG Jin-yong HU Shuang-qi CAO Xiong 作者单位：中北大学环境与安全工程系,太原,030051刊 名：工业安全与环保 PKU英文刊名：INDUSTRIAL SAFETY AND ENVIRONMENTAL PROTECTION年，卷(期)：32(5)分类号：X5关键词：多点同步起爆网络 传爆药柱 起爆点 输出波形
清华大学数值分析实验报告 篇5
一、实验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
物流模拟实验报告 篇6
班 级：组（企业）号：姓 名：学
号 ：
电子商务1001班 no.xx 1002110105 目 录
1、概述.............................................................3 1.1实习目的..........................................................................................................................3 1.2实习手段..........................................................................................................................3 1.3实习进程安排.................................................................................................................3 1.4实习原理..........................................................................................................................3 1.5实习的平台.........................................................................................................................3
2、实习过程及实习内容...........................................................................................................3 2.1实习主要阶段性工作安排..........................................................................................3（1）实习准备.............................................................................................................3（2）第一阶段............................................................................................................3（3）第二阶段............................................................................................................4（4）第三阶段............................................................................................................4（5）第四阶段............................................................................................................4 2.2实习收获、感想、认识、评价等............................................................................4
3、实习总结.................................................................................................................................6
1、概述
1.1实习目的： 通过对奥拓物流软件进行模拟，掌握物流运作的基本流程，从而让我们了解物流的知识，更好的面对以后关于物流上面的工作和问题，小组成员能够通过这次实习能够熟悉掌握软件的使用，提高我们在物流管理中的信息技术能力的水平提高,学习现代物流思想，同时培养学生的实际操作能力和决策能力。1.2实习手段：
在15教的机房运用奥派物流软件进行实习1.3实习进程安排： 周一到周四下午 14：00-18：00 四个小时，师生答疑和上机共同进行。周五撰写实习报告，下周二交到指定教室。1.4实习原理： 1.5实习的平台：
奥派物流实践教学平台
2、实习过程及实习内容
2.1实习主要阶段性工作安排（1）实习准备： 在实习开始前，我们根据老师的要求，五个人一个小组，其中每个成员都分配好自己的角色，包括采购商、采购商仓储、供应商、供应商仓储、运输五个角色，分配好后输入相应的ip地址，进入奥派物流软件平台，个人注册完成，以及各个模块的注册完成，等待老师审核完成即可。（2）第一阶段：各个模块的注册完成，建立自己的集团以及自己的分属公司—张洁集团采购商公司。（3）第二阶段：
建立相应的客户关系。根据自己的角色，采购商应该和供应商和采购仓储建立相应的客户关系。（4）第三阶段：
根据采购的商品进行相应的汇款、运输等。（5）第四阶段： 交易完成，任务结束。2.2实习收获、感想、认识、评价等
在这次的物流角色的模拟中，自己的角色是采购商，根据自己的职责，做出相应的措施.在这次的实习过程中有顺利成功的地方也有不顺利失败的地方，成功的地方让我熟练的记住这个知识点，失败的有疑问的地方让我去寻找正确的答案，然后更好的丰富自己的知识。通过这一个周的实训加上这学期物流与供应链的知识的学习，让我把这两门很好的结合，让我更加进一步的认识到物流的相关知识，让我对物流的定义、基本功能、操作流程等有了很明确的概念和理解。首先，了解并学习了奥派物流软件平台的操作，从学生注册、模块注册，到各个公司的建立，采购商和供应商达成交易，到最后的商品出库等等，一系列的操作流程，让我从一个根本一无所知到后来正确的明白各个环节在物流过程中的作用和地位，这对我来都是收获颇丰和受益匪浅的。
现在就自己是一个采购商的角色，来具体的说说这次实验中自己的任务的完成情况。(1)因为自己的角色是采购商1002110107，所以建立属于自己组集团的分属公司—张洁集团采购商公司。自己采购商与供应商1002110112和采购仓储1002110109建立客户关系。（2）采购商采购商品，进行询价，询价成功供应商1002110112拟写合同，并且双方合同达成一致。（3）收到张洁集团供应商公司1002110112的缴费单信息，对自己采购的商品进行汇款。（4）供应商仓储1002110116进行发货，采购商将采购到的商品进行入库，转给张洁集团采篇二：第三方物流模拟实验报告
第 三 方 物 流 模 拟
实 验 报 告
实验项目：了解第三方物流实验指导教师：徐伟 班级： 物流0801学号： 0508106112姓名： 陈益群实验时间： 11月23号到11月27号第三方物流模拟实验报告 实验名称：第三方物流模拟
实验目的：随着第三方物流在社会上的推广，社会成立了越来越多的第三方物流公司。通过第三方物流模拟实验，我们可以掌握物流系统的各个环节的运作。同时，通过自己成立第三方物流公司，并在虚拟的市场上接受订单进行操作掌握第三方物流公司的运作流程。实验任务：(1)掌握第三方物流运作的基本流程(2)掌握物流系统中各个子系统的基本运作
(3)灵活运用所学的知识完成第三方物流的各项操作 实验仪器设备：计算机、物流管理模拟软件
实验原理：第三方物流tpl(lhird party logistics,又称物流代理)，是指由物流劳务的供方、需方之外的第三方去完成物流服务的运作方式。第三方提供物流交易双方的部分或全部物流功能的外部服务，把物流（商品实体的转移）以电子传递的形式通过网络来完成的模式。《3pl教学模拟平台》模拟了整个物流的过程，前台包括：管理中心（公司的一些资源管理、订单管理等）、调度中心（对订单进行调度配送）、运输中心（选择车辆运输配送的方案），仓储中心（仓库的进、出库的操作），物流市场（资源的购买）；后台包括班级、老师、学生管理、系统参数设置等
实验过程：这次的实习时间为期一周，扣除了周三和周五的下午。
周一的早上，早早的起床，吃完早饭，怀着期待的心情，我们来到了经管楼的五楼机房里，开时我们的实习.我们上机练习大致分为：
一关于公路`铁路`航空`等方面的报价
二是订单（长期订单、散单、外包）的受理，运输工具的选择和发货 三是关于仓库的出入库运单的处理 四是财务的上报和结算
其中包含了广告的发布，人员的安排，车辆的购买，仓库的购买和选址等。还有就是特殊情况的处理，如仓库货物被盗，车辆事故的处理。进入系统我们每个人都可以需申请了一个公司，注册资金是1000万。我申请的公司名称叫龙腾物流公司，账号是0508106112。我把总部方在福州，因为随着海西建设的号召，福建正在腾起，对物流的发展建设更是需要，而且容易形成融断,所以我判定在这里注册生意一定多。第一天，实验刚开始做还是很不上手，女生们总不敢轻易动，就怕做错。注册以后，我看到物流市场上有好多公共订单，刚开始我想，接的单子越多赚的钱也一定越多，而且我对公路运输情有独钟，所以只要看见是公路的单子不管三七二十一，一律拿下。接着我就按照订单上所需要的车辆进行买车辆，不过由于怕资金不够，我只敢买几辆车，其他都是租的。接着就开始买仓库，刚开始我买了两个设在福建的仓库，租了一个设在黑龙江的仓库。然后就是人员的配备。等一切弄好以后便开始营业了。就这样，我第一天就做了50多单生意
第二天上去一看我的财务报告，心里可高兴了，我的资金一下子多了100多万了，也就是说我第一天就赚了100多万。龙腾物流公司的排名一下子越升到了前六名。于是，就这样，接市场上的散订单，租车，配备人员，第二天很快过去了！
第三天一早，我依旧是进行财务上报，发现公司的排名下子被挤到后几名，自己的公司经过我昨天一天的营业，虽然没亏，可是居然分文不赚，也就是说昨天我赚的钱和亏的钱相平了。这下我傻眼了，怎么会这样呢呢？别人的资金都已经增了六七百万了。后来听老师说，不是每个定单都赚钱的，接太多的散单如果没有及时处理容易形成死单，这样就根本赚不了钱。而且报价也很重要。
吸取了这个教训以后，我开始对进行报价，可是报价是个相对比较困难的步骤，要经过：成本计算?报价，本来我以为报价低一点可以有更多的客户，可是我面对的是个机器啊，他没有“人性”，于是我报了老高的价。对公路，铁路，航空，仓库，搬运进行报价完后，我开始有了自己的长期合作者，而且在接单时，我都仔细查看单子的报价，以及他的价值。不看不知道一看吓一跳，举个例子从广州空运一批2吨多的服装到北京他的报价是3200元。你说这样的生意傻子才会接啊，想想以前我都是看都不看就接的，那亏的多万也就不为过了。有了长期合作者后我的订单都能及时完成，不禁感慨：“终于可以从自己的仓库出货而不用那么频繁的到处跑叫别人发货了
为了使资金增长，只好接多的单子（这次都是经过我查看确保赢利才接的），尽量每辆车都装满，有空的就把同路线的货拼装，做到车无虚发。在这期间，我还接了几个外包的订单，虽然具有一定的风险，但我想，越有风险越能赚钱，于是，经过价格修改到最后的签合同，我在外包上也赚了一大笔钱.这样到第四天的时候经过财务结算，有起色，资金多了600多万。生意越到后面越难做，很少接到好的单子。毕竟全班人都在抢单子，在接单子的时候我发现航空和铁路的单子很少有人接，也许是传统思维的影响。可我看了一下航空和铁路的报价，利润也可观，于是又“买”了n个搬运工，把报价调得老高，做起了航空和铁路的生意。有个好的客户和多元话经营，很快赢利就翻了翻。
到后面，有了长期合作者之后，我的仓库终于有了一定的库存，从出仓通知->备货->签发运输单->出库装货->出库登记->费用结算，这都需要一定的时间，而且当仓库没及时盘点，还会导致货物出库发生装备事故，当然，入库操作与货物出库操作基本相同 还有一个就是广告的效应，虽然在模拟系统里和现实有很大的区别，可是它却提示我们在当今的物流当中是不可缺少的，即使的的物流服务再好，价格再优惠，没有人知道那也等于零。很典型的例子，在一些车辆上我们经常可以看到某某物流公司，其实广告做到一定程度也是诚信的一种体现，而诚信是从商之本，可以想象广告对一个企业的影响。
等到第五天时，我的资金多了800多万，看一下自己的成绩测评，只有70多分，看一下评分的标准，原来要多买车辆，多买仓库就可以加分，租的仓库、租的车辆不加分。资金的排名只占评分的一部分，这些后来我都做到了，可是最让我头疼的是那个报价，之前，我报了老高的价为自己赢的了暴利，可是，报太高的价是会扣分的，于是矛盾出现了，如何平衡好赚钱与得分是关键！遗憾的是，我虽然在整个实验中对报价进行了四次的价格调整，然而系统依旧无法为我加分。
实验中出现的问题和不足：
1：系统无法让我们与其他同学合作经营注册的物流公司，否则容易造成系统瘫痪 2：网络太卡了，一直掉线 3：太晚买仓库导致自己的仓库没有任何的货物，货物集中在几个一开始就买好仓库的同学仓库里，造成那几位同学仓库货物过多，他们来不及给我们发货，容易形成死单 3：空运，海运，铁路等很难完成订单
4：报价很难决定，报高了可以为自己赚更多的钱，但会客户更少，且系统为法为自己的报价加分，可是如果报低了，可以加分却无法为自己的公司赢得相应的利润 5：对租的车辆和招聘的司机没有及时跟踪还返车辆及解雇司机，造成资金的损失 6：对车辆没有合理的安排，造成空车移动，增加物流成本
7：买的仓库不够大导致货物无法入库
8：合并单太多，以致不知道先找谁发货，容易形成死单
实习心得：经过这次实习，确实掌握了不少东西。以前学的课本知识，只是对物流的一个表面的了解。而且都是单一的，要么学运输，要么学仓储。而这次上机实习，是把这些进行了融合和连贯，将这些理论与实际的操作相结合,在实践中提高运用知识的能力,从而对物流操作流程有了更进一步的了解。实习是短暂的，可是他对我们的“思维方式”有很深的影响，这对我们以后走上工作岗位，我想是有很大的帮助
物流是支持经济生长的基础性产业，是经济生长到肯定历史阶段的一定产品，在我国更是极具生长潜力的新兴产业。随着物流服务的睁开，我国物流人力资源，尤其是物流高级管理人才紧张稀缺，远远不能满足今世物流生长的需求，物流人才缺乏是制约我国物流大规模生长的瓶颈，这就要求我们，新一代的大学生，扎扎实实学好仓储、配送、包装等每一块的专业知识，让自己成为复合型物流人才篇三：物流中心软件模拟实验报告 《物流中心设计与运作》
logistics center design and operation 实验报告
院系名称：商学院
实验班级： 11物管一班 姓 名： 陈玉玲
学 号： 20114211027 指导教师：王 力
地震数值模拟实验报告 篇7
一、地震波场数值模拟简介
地震波模拟是根据给定地下介质的结构模型和相应物理参数来模拟地震波的传播过程,从而研究地震波在地下介质中的传播规律。由于模拟过程中可直观、形象、动态地显示地震波动力学和运动学传播特征,非常容易调动学生的学习兴趣和求知欲望,可以收到事半功倍的教学效果。一般来说,地震波模拟可分为物理模拟和数值模拟两种方法。物理模拟是在实验室内将野外的地质构造和地质体按照一定的比例制作成物理模型,然后利用超声波或激光超声波等方法对野外地震勘探方法进行模拟;数值模拟就是利用有限差分、有限元等数值方法求解波动方程,从而获得已知模型的地震波传播。采用物理模拟方法存在费用高、选材困难等缺点,并且不适合于课堂理论教学;而数值模拟只需要一台较高计算速度的计算机就可全部解决问题,既简单方便,成本又低,且非常适合于课堂理论教学。
地震波场数值模拟方法主要有几何射线法和波动方程法。几何射线法也称为射线追踪法,其主要理论基础是,在高频近似条件下,地震波的主能量沿射线轨迹传播,即根据地震波的传播规律确定地震波在实际地层中传播的射线途径,并运用惠更斯原理和费马原理来重建射线路径,利用程函方程等计算射线的旅行时间。射线法的主要优点是概念明确、显示直观、运算方便、适应性强,缺点是旅行时的计算在一定程度上是近似的,特别是对复杂构造进行三维射线追踪时繁琐且误差较大。波动方程数值模拟方法是以弹性(粘弹性)理论及牛顿力学为基础的,求解双典型偏微分方程-波动方程为手段的一种数值模拟方法。这种方法不仅能保持地震波的运动学特征,而且还能保持地震波的动力学特征。根据地震波场数值模拟所采用的这两种不同方法的特点,我们在地震勘探原理课程讲授过程中采用波动方程弹性波数值模拟方法,这种方法的模拟结果可以使学生更好地了解地震波在传播过程中的运动学和动力学特征。目前用于波动方程数值模拟的方法主要有限差分法、有限元法、虚谱法等,在这些方法中有限差分法在计算精度、计算效率上都占有较大的优势,因此我们在数值模拟时采用了有限差分模拟方法[5]。
二、应用实例
地震波数值模拟在地震勘探原理课程教学中的应用比较广泛,在大部分的知识点讲解中都可采用数值模拟来对一些方法原理进行动态演示。在本文中主要通过两个知识点的应用来说明数值模拟方法在地震勘探原理课程教学中的效果。
(一)地震波的传播及界面处的反射、折射、透射和波形转换
在讲解地震波场的基本知识一节中,主要讲解地震波的传播特点,其中一些名词如波前、波后、球面波、平面波等,学生理解起来比较抽象,一些概念只能靠死记硬背。在讲解地震纵、横波的传播特点时,对于波在界面处发生的反射、透射、折射以及波形转换等,由于学生对波场的概念没有直观的认识,所以老师讲解起来也非常抽象和费力。在这种情况下,我们将数值模拟技术引入到课堂教学中,通过开发地震波场实时模拟软件,只需设置好模型和参数,软件将以动画的形式动态展示波场传播的全过程。
下面以两层介质为例进行说明:模型大小为600m×300m,网格大小为 △x =△z ?1m ,时间采样间隔为△t =0.1ms。模型上层纵波波速为3 000m/s,横波波速为1 732m/s,密度为1.8g/cm3;下层纵波波速为4 000m/s,横波波速为2 300m/s,密度为2.0g/cm3。图1为两层介质模型波场传播快照(t=20~100ms),从图中可见,当t=20ms时纵波和横波以震源为中心向下传播,此时由于传播时间较短,P波(纵波)和S波(横波)在图上还不能完全分清楚。当t=40ms时,从图中可以看出,由于P波波速较快传播在前,S波在后,P、S波都是以震源为圆心的同心半圆。根据波场快照很容易让学生理解波前和波后的概念。当t=40ms时,从图1(d)中可以看出,P波到达界面后发生反射和透射,同时产生转换S波。当t=80ms时,从图1(g)中可以看出,此时S波到达分界面后也发生反射和透射,同时产生转换P波。在图1(h)中对于不同类型的地震波做了相应的标识,从图中可以看出对于两层介质中弹性波传播时波场也比较复杂,有直达P和S波、透射P和S波、反射P和S波、PS转换波和SP转换波。由于上述波场较为复杂,在学生认识波场特征时容易混淆,鉴于这种情况,我们对于数值模拟技术做了相应的改进,在模拟过程中对P波和S波进行了分离,图2为两层模型t=90ms时分离后的P波和S波波场快照,从图2中可以看出分离后的波场快照波场更加清晰,学生在认识波场传播特性时也会更加容易理解和掌握。
(二)瑞雷波传播特性
瑞雷波方法是近年来发展较为迅速的一种工程地球物理方法,在讲解瑞雷波传播特点时主要是通过其传播速度、能量衰减和频散特点等三个方面展开的。由于瑞雷波属于面波,与纵、横波所属的体波相比其传播特点有较大的差别。在瑞雷波传播特性这一节课程的讲解过程中,我们就充分运用了数值模拟方法[6]。图3为均匀半空间模型下利用数值模拟方法获得的t=260ms时瑞雷波波场传播快照,其中模型参数为:纵波速度Vp =1000ms,横波速度Vs =577m/s,密度 ρ=2.0g/cm3。从图3中可以看出,瑞雷波沿地表传播,其传播速度和横波波速比较接近。从波场传播快照中计算可得瑞雷波的传播速度为Vr ≈ 531m/s,而模型中已知的横波速度为Vs =577m/s,瑞雷波速度和横波速度两者相比为0.92倍。由于所给模型为泊松体,从理论上分析可知,瑞雷波传播速度和横波传播速度比值为0.92,这样通过数值模拟,就很容易将瑞雷波和横波的关系讲解清楚,学生理解起来非常方便。
而对于瑞雷波传播深度的知识点,从理论分析可知,瑞雷波的传播深度为一个波长且波的能量呈指数衰减。对于瑞雷波的这个性质,由于和体波传播特性差别较大,学生们理解起来存在一定的困难。我们通过数值模拟,从波场传播快照上就可以对这个问题给予清楚说明。从图3的快照中可以测量出瑞雷波在模型空间中的传播深度约为50m,在50m以下基本上看不到瑞雷波的存在,通过计算模拟时瑞雷波传播的最大波长为44.2m,这充分说明瑞雷波的传播深度为一个波长。为了说明瑞雷波能量呈指数衰减的问题,我们从波场快照上提取了一道瑞雷波数据并进行指数拟合,得到该模型瑞雷波的能量衰减公式为: 1.370.055zE e?? ,其拟合曲线如图4所示。从图中可见,在深度为44.2m处能量衰减到原来的0.8%,在22.1m(即半个主波长)处能量衰减到原来的7.9%。因此,对于瑞雷波用于工程地质勘察中,通常取二分之一波长为其有效勘探深度,就可以得到较合理的解释了。
三、结论