单轮对纵向动力学数值分析
单轮对纵向动力学数值分析
严晓明，罗世辉，马卫华
（西南交通大学牵引动力国家重点实验室，四川成都６１００３１）
摘要：轮对的纵向振动会影响机车车辆动力学性能，而且是轮轨非正常磨耗的一个重要因素。但是机车车辆动力学研究中，对轮对的纵向动力学特点的研究往往被忽略。文章建立了一个包括ｘ方向的运动、轮对的摇头、点头扰动和轮对的横移的４自由度单轮对计算模型，并对该模型进行数值仿真，研究其纵向振动现象。最后讨论了系统参数对纵向动力学行为的影响，认为一系纵向刚度、轴重和黏着系数对纵向振动影响很大。
关键词：轮对；纵向振动；数值仿真；动力学；参数；影响中图分类号：Ｕ２６０．３３１＋．１
文献标识码：Ａ
文章编号：１６７２－１１８７（２００５）０６－００２２－０３
Ｎｕｍｅｒｉｃａｌａｎａｌｙｓｉｓｏｆｌｏｎｇｉｔｕｄｉｎａｌｄｙｎａｍｉｃｓｏｆｓｉｎｇｌｅｗｈｅｅｌｓｅｔ
ＹＡＮＸｉａｏ－ｍｉｎｇ，ＬＵＯＳｈｉ－ｈｕｉ，ＭＡＷｅｉ－ｈｕａ
（ＴｒａｃｔｉｏｎＰｏｗｅｒＳｔａｔｅＫｅｙＬａｂｏｒａｔｏｒｙ，ＳｏｕｔｈｗｅｓｔＪｉａｏｔｏｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｅｎｇｄｕ６１００３１，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：Ｔｈｅｌｏｎｇｉｔｕｄｉｎａｌｖｉｂｒａｔｉｏｎｏｆｗｈｅｅｌｓｅｔｅｆｆｅｃｔｓｔｈｅｄｙｎａｍｉｃｓｐｅｒｆｏｒｍａｎｃｅｏｆｒａｉｌｗａｙｖｅｈｉｃｌｅ，ａｎｄｉｔｍａｙｂｅｃａｕｓｅａｂｎｏｒｍａｌｗｈｅｅｌ／ｒａｉｌｃｏｎｔａｃｔｆａｔｉｇｕｅｐｒｏｂｌｅｍ．Ｂｕｔｔｈｅｌｏｎｇｉｔｕｄｉｎａｌｄｙｎａｍｉｃｓｂｅｈａｖｉｏｒｏｆｗｈｅｅｌｓｅｔｉｓｏｆｔｅｎｎｅｇｌｅｃｔｅｄｗｈｉｌｅｃａｒｒｙｉｎｇｏｕｔｒａｉｌｗａｙｖｅｈｉｃｌｅｄｙｎａｍｉｃｓａｎａｌｙｓｉｓ．Ｔｈｅｄｙｎａｍｉｃｓｃａｌｃｕｌａｔｉｏｎｍｏｄｅｌｆｏｒｓｉｎｇｌｅｗｈｅｅｌｓｅｔｉｎｃｌｕｄｉｎｇｌｏｎｇｉｔｕｄｉｎａｌｄｉｓｐｌａｃｅｍｅｎｔ，ｙａｗ，ｐｉｔｃｈａｎｄｌａｔｅｒａｌｄｉｓｐｌａｃｅｍｅｎｔｉｓｓｅｔｕｐｉｎｔｈｉｓｐａｐｅｒ．Ｓｉｍｕｌａｔｉｏｎｏｆｔｈｉｓｍｏｄｅｌｉｓｄｏｎｅ，ａｎｄｔｈｅｌｏｎｇｉｔｕｄｉｎａｌｖｉｂｒａｔｉｏｎｐｅｒｆｏｒｍａｎｃｅｉｓｉｎｖｅｓｔｉｇａｔｅｄ．Ｔｈｅｅｆｆｅｃｔｏｆｔｈｅｓｙｓｔｅｍｐａｒａｍｅｔｅｒｓｉｓｄｉｓｃｕｓｓｅｄ，ｔｈｅｌｏｎｇｉｔｕｄｉｎａｌｖｉｂｒａｔｉｏｎｉｓｉｎｆｌｕｅｎｃｅｄｈｉｇｈｌｙｂｙｐｒｉｍａｒｙｌｏｎｇｉｔｕｄｉｎａｌｓｔｉｆｆｎｅｓｓ，ａｘｌｅ－ｌｏａｄａｎｄｃｏｅｆｆｉｃｉｅｎｔｏｆｆｒｉｃｔｉｏｎ．
Ｋｅｙｗｏｒｄｓ：
ｗｈｅｅｌｓｅｔ；ｌｏｎｇｉｔｕｄｉｎａｌｖｉｂｒａｔｉｏｎ；ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ；ｄｙｎａｍｉｃｓ；ｐａｒａｍｅｔｅｒｓ；ｉｎｆｌｕｅｎｃｅ电力机车与城轨车辆
ＥｌｅｃｔｒｉｃＬｏｃｏｍｏｔｉｖｅｓ＆ＭａｓｓＴｒａｎｓｉｔ
Ｖｅｈｉｃｌｅｓ第２８卷第６期２００５年１１月２０日Ｖｏｌ．２８Ｎｏ．６Ｎｏｖ．２０ｔｈ，２００５
收稿日期：２００５－０６－１３
作者简介：严晓明（１９８０－），男，在读硕士研究生，研究方向为机车、城轨车辆动力学理论仿真及应用。
０引言
在轨道车辆动力学研究中，总是事先假定车辆以
固定的名义速度前进，忽略纵向上速度变化的影响，但实际上车辆的前进速度应该是名义的滚动速度上叠加一个随机分量。在文献［１］中，
借助于一个整车模型，观察到了轮对相对于构架实际上始终存在纵向振动，而且在一定的条件下，这种振动可能发展为强烈的共振，如果牵引装置纵向连接刚度很大，它还会引起车体的点头和纵向颤抖，大大恶化了整车的垂向平稳性。另一方面，轮对相对构架的强烈振动，将在轮轨接触点产生强烈的粘滑作用，使轮轨接触点承受很大的疲劳载荷，成为车轮非正常磨耗（即剥离）以及某些波长轨道波磨的一个重要因素［１］。
长期以来，为了实现更高的运行速度和更好的平稳性，铁道车辆的动力学分析一直非常关注横向动力
学问题。机车车辆纵向动力学问题的研究，主要集中在列车纵向动力学和机车纵向动力学，机车纵向动力学研究主要集中在传动系统的扭转刚度能否满足大的起动黏着利用，防止动轮起动时和运行时滑动和空转［１］，对轮对纵向振动问题的研究一直不太为人注意。
本文从一个最简单最基本的４自由度单轮对模型出发，通过对该模型的计算机仿真，得到轮对时域范围的纵向振动情况，并讨论了一些主要的系统参数对纵向振动的影响。
１动力学模型
１．１模型自由度
系统数值仿真的计算动力学模型如图１所示，坐
标系的取法为：以车辆的前进方向为ｘ轴，ｙ轴平行于轨面指向右方，ｚ轴垂直轨道平面向下。建模时，假设轮对在刚性、平直的轨道上运行，轮对具有一系纵向、
－２２－
严晓明等?单轮对纵向动力学数值分析?２００５年第６期
横向刚度和阻尼。轮对在平直钢轨上运动时，由于轮对的踏面有斜度，轮对在横移时，伴随着侧滚，同时其质心产生摆动；轮对摇头时，质心作上下运动。由于轮轨接触几何关系已知，侧滚和质心沉浮的运动规律均可求得，因此它们不是独立的自由度。所以系统具有４个自由度，包括ｘ向的运动ｘｗ，轮对的摇头ψｗ，点头扰动βｗ和轮对的横移ｙｗ。
图１单轮对模型
１．２轮轨接触模型
轮轨接触采用ＪＭ１磨耗型机车踏面和６０ｋｇ／ｍ钢
轨的配合，轮轨接触几何关系是轮对横移量的非线性函数，这些几何参数如图２所示。
（ａ）滚动圆半径与轮对横移量（ｂ）轮对侧滚角与轮对横移量
图２轮轨接触几何关系
２模型的运动微分方程及求解
对轮对做受力分析（见图３），可以得到系统的力学
方程：
图３轮对受力图
ｍｘ
!ｗ!＝ＦＣｘＬ＋ＦＣｘＲ＋ＦＳＵＳＰｘ
（１）ｍｙ
!ｗ＝ＦＣｙＬ＋ＦＣｙＲ＋ＦＮｙＬ＋ＦＮｙＲ＋ＦＳＵＳＰｙ（２）Ｉｙβ!ｙｗ＝ＦＣｘＬｒＬ＋ＦＣｘＲｒＲ＋ＭＣｙＬ＋ＭＣｙＲ＋ｘｗＬ（ＦＣｚＬ＋ＦＮｚＬ）－ｘｗＲ（ＦＣｚＲ＋ＦＮｚＲ）（３）
Ｉｚψ!ｗ＋Ｉｙβ"ｗΦ"ｗ＝（ＦＣｘＬ－ＦＣｘＲ）ａ＋ＭＣｚＬ＋ＭＣｚＲ＋ａψｗ（ＦＣｙＬ＋
ＦＮｙＬ－ＦＣｙＲ－ＦＮｙＲ）＋ＭＳＵＳＰ
（４）
式中：!ｘｗｉ—左右轮轨接触点超前量，（ｉ表示左右轮，ｉ＝Ｌ，
Ｒ）；ＦＣ，ＦＮ—蠕滑力和法向力；ＦＳＵＳＰ，ＭＳＵＳＰ—悬挂力和悬挂力矩；Φｗ—轮对侧滚角；ｍ—轮对质量；ａ—两轨接触点距离之半；Ｉｊ—ｊ方向的转动惯量（ｊ＝ｘ，ｙ，ｚ）。
超前量ｘｗｉ＝ψｗｔａｎ（δｉ）ｒｉ，其中δｉ为轮轨接触角；ｒｉ为滚动半径［２］。
关于轮轨蠕滑力的计算，首先用Ｋａｌｋｅｒ线性理论计算，然后采用Ｓｈｅｎ－Ｈｅｄｒｉｃｋ－Ｅｌｋｉｎｓ理论进行非线性修正［３］。蠕滑率计算公式如下：
ζｘｉ＝１ｖｉ［ｖａ－ｒｉ
（ｖａｒ０
＋β"ｗ）ａψ"ｗζｙｉ＝１ｖｉ
［ｙ"ｗ＋ｒｉΦ"ｗ－ｖψｗ］／ｃｏｓ（δｉ±Φｗ）ζＳＰｉ＝１ｖｉ［!!!（ｖａｒ０＋β"ｗ）ｓｉｎ（δｉ±Φｗ）＋ψ"ｗｃｏｓ（δｉ±Φｗ）］式中：ζｘｉ，ζｙｉ，ζＳＰｉ—纵向、横向、自旋蠕滑率；ｒ０—名义滚动半径；ｖａ—轮对前进速度；ｖｉ—左右轮的名义速度，ｖｉ＝１２（ｖａ＋ｒｉｒ０
ｖａｃｏｓψｗ）。左右接触斑上的法向力根据约束方程推导出，采用文献［４］所用方法。最终系统的微分方程整理成：
［Ｍ］｛ｘ
%｝＋［Ｃ］｛ｘ"｝＋［Ｋ］｛ｘ｝＝｛Ｆ｝式中：［Ｍ］为质量矩阵；［Ｃ］为阻尼矩阵；［Ｋ］为刚度矩阵；｛ｘ｝为坐标向量；｛Ｆ｝为包括激励的广义力矩阵。
方程求解采用四阶龙格—库塔法，计算步长为０．１ｍｓ，计算过程中加入了德国高干扰轨道谱激励，包括水平不平顺和方向不平顺。
计算中的系统参数具体值为：轮对质量２４８０ｋｇ，轴重２１ｔ，名义滚动半径０．５２５ｍ，轮轨接触点距离之半０．７４６５ｍ，轮对摇头转动惯量４９８ｋｇ?ｍ２，轮对点头转动惯量１０６ｋｇ?ｍ２，一系纵向刚度１．２ＭＮ／ｍ，一系横向刚度０．３０６ＭＮ／ｍ，一系纵向、横向阻尼６ＫＮ?ｓ／ｍ。
３
纵向振动的数值仿真结果
３．１
纵向振动现象
仿真过程中，加入了水平和方向不平顺，轮对在前
进时必然伴随着摇头和侧滚运动，产生了轮对质心在纵向上的振动。图４（ａ）为初始速度为８０ｋｍ／ｈ时的纵向加速度数值仿真结果。从图中可以明显观察到纵向的振动现象，加速度在－０．７ｍ／ｓ２到０．６ｍ／ｓ２范围内产生了振动。轮对的实际前进速度是８０ｋｍ／ｈ的初始速
度叠加一个随机速度分量。纵向的振动还表现在黏着
－２３－
电力机车与城轨车辆?２００５年第６期
系数的波动，图４（ｂ）是黏着系数随时间变化的曲线。从图中可以看出，黏着系数波动严重，经常超过了最大许用黏着系数０．１５，也就是说左轮经常出现滑动现象。一旦左轮或右轮出现了滑动，必然导致轮对在纵向上振动的加强，可以想象，如果左右轮都出现滑动的话，整个轮对的振动加速度就会明显增大。如果机车车辆在运营中，轮对经常发生这样强烈的纵向振动，轮轨接触点处就产生强烈的粘滑作用，轮轨接触点承受的疲劳载荷加大，恶化了轮轨接触情况，很可能导致轮对出现非正常磨耗、轨道出现某些特定波长的波磨。
（ａ）纵向加速度时间响应曲线
（ｂ）左轮黏着系数时间变化曲线
图４８０ｋｍ／ｈ下纵向振动仿真结果
３．２参数对纵向振动的影响
由于轮对的纵向振动总是存在的，有时候加速度振幅会比较大，加大了轮轨接触的动态载荷，产生了一些不利于轮轨接触的影响。所以有必要讨论系统参数对纵向动力学特性的影响，为提出抑制轮对纵向振动方法提供判断依据。在讨论参数影响时，采用纵向振动加速度的最大值和纵向振动加速度的３σ值作为评价振动的指标（σ值是振动的标准差，σ值的大小能够反映振动的剧烈程度）。
３．２．１一系定位刚度的影响
一系轴箱定位刚度对纵向振动的影响很大，本文在保持其他计算条件不变的情况下，对一系定位刚度在初值１．２×１０６Ｎ／ｍ附近进行仿真计算，计算速度都为８０ｋｍ／ｈ，计算结果见图５。图中可以看出纵向定位刚度的影响是很大的，有些范围内即使是刚度值只改变０．１ＭＮ／ｍ，纵向振动加速度的３σ值都可能发生很大的变化。而且纵向加速度的大小不是严格的随纵向定位刚度增大或减小而线性变化，没有明显的变化趋势。所以实际在确定系统的参数时，首先在保证整车的稳定性、牵引和制动的前提下，确定出一系纵向刚度。为了减小纵向振动，可以在初选的一系纵向刚度附近模拟出纵向振动的曲线，根据所得到的曲线，避开会引起较大振动的一系纵向定位刚度值，对一系定位刚度值进行优化选取；而且建议一系纵向刚度的选取尽量避开对纵向振动影响变化较大、较为敏感的区域。
图５一系纵向定位刚度的影响曲线
３．２．２簧上质量的影响
不同的簧上质量将引起轮轨接触处不同的法向力，从而导致接触斑的形状不同，进而影响到纵向蠕滑率和蠕滑力。在保持相同的计算条件下，分别对１６，２１，２５，３３ｔ不同的簧上质量进行仿真计算。计算时，轮对的质量和转动惯量这些条件都不改变，只是改变簧上质量，计算结果如表１所示。从表中可以看出，随着轴重的增加，轮对纵向振动的幅值有所减小，但是３σ这个统计值却增大了，反应出振动的强度增大了。
表１不同簧上质量时纵向振动情况
簧上质量／ｔ１６２１２５３３
纵向振动加速度／ｍ?ｓ－２－０．７１５９４－０．６７６０８－０．６７５２６－０．６６８６７振动加速度３σ值／ｍ?ｓ－２０．４６２３０．４７６００．４７９１０．４７９８
３．２．３许用黏着系数的影响
在同样的仿真条件下，分别取许用黏着系数为０．１５，０．２５，０．４的情况进行仿真计算，计算结果如表２所示。从表中可以看出，随着许用最大黏着系数的增大，振动的幅度增大了，而且振动也更加地剧烈。从仿真结果还可以看出，由于黏着系数的提高，黏着系数不容易超过最大黏着系数，轮对更不容易出现空转。但是由于振动的加剧，使得轮轨接触区的疲劳动载荷加强了，
加剧了踏面与轨面的磨耗。（下转第２７页）－２４－
表２不同许用黏着系数时纵向振动情况
许用黏着系数大小０．１５０．２５０．４纵向振动加速度／ｍ?ｓ－２－０．６７６０８－０．７６５３５－２．８１６８１纵向加速度的３σ值／ｍ?ｓ－２０．４７６００．５３９９０．７１０２
４结论
（１）纵向振动总是存在的，不仅表现在纵向加速度上的振动，纵向蠕滑率和黏着系数也是剧烈变化的。振动严重时会使黏着系数经常超过最大黏着系数，导致轮对空转，并在纵向上滑动。
（２）轴箱纵向定位刚度对纵向振动影响很大，增大轴重和最大黏着系数，会导致纵向振动的加剧。
（３）强烈的纵向振动会在轮轨接触点处产生大的粘滑作用，使得轮轨接触点承受的疲劳载荷加大，恶化了轮轨接触情况。
参考文献：
［１］罗世辉，金鼎昌，陈清．轮对纵向动力学——
—一种尚待证实的动力学现象［Ａ］．Ｉｎ：Ｐｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅ６ｔｈＩｎｔｅｒｎａｔｉｏｎａｌＣｏｎｆｅｒｅｎｃｅｏｎＦｒｏｎｔｉｅｒｓｏｆＤｅｓｉｇｎａｎｄＭａｎｕｆａｔｕｒｉｎｇ［Ｃ］．ｘｉ＇ａｎ：［ｓ．ｎ．］，２００４：４２－５１．［２］ＷｉｃｋｅｎｓＡＨ．ＦｕｎｄａｍｅｎｔａｌｓｏｆＲａｉｌＶｅｈｉｃｌｅＤｙｎａｍｉｃｓ［Ｍ］．Ａｍｅｒｏｎｇｅｎ：Ｓｗｅｔｓ＆ＺｅｉｔｌｉｎｇｅｒＰｕｂｌｉｓｈｅｒｓ，２００３：４４．
［３］翟婉明．车辆—轨道耦合动力学［Ｍ］．北京：中国铁道出版社，２００１：６９－７８．
［４］ＭｏｈａｎＡ，ＡｈｍａｄｉａｎＥ．ＮｏｎｌｉｎｅａｒＩｎｖｅｓｔｉｇａｔｉｏｎｏｆｔｈｅＥｆｆｅｃｔｏｆＰｒｉｍａｒｙＳｕｓｐｅｎｓｉｏｎｏｎｔｈｅＨｕｎｔｉｎｇＳｔａｂｉｌｉｔｙｏｆＡＲａｉｌＷｈｅｅｌｓｅｔ［Ａ］．Ｉｎ：２００４ＰｒｏｃｅｅｄｉｎｇｓｏｆＡＳＭＥ／ＩＥＥＥＪｏｉｎｔＲａｉｌＣｏｎｆｅｒｅｎｃｅ［Ｃ］．Ｂａｌｔｉｍｏｒｅ：［ｓ．ｎ．］，２００４：５３－６１．
洪常勤等?钢丝螺套在机车维修及改造中的应用工艺?２００５年第６期
值的改变又首当其冲。
图３钢丝螺套型面结构
通常情况下，钢丝螺套的使用者均认为钢丝螺套的型面尺寸Ｂ为一固定值，这无疑给使用者带来了负面影响。而实际上，使用者可根据自己的需要来选定Ｂ值，但应注意所选定的Ｂ值一定要在相关标准规定的公差带内，同时还应充分考虑其互换性。也就是说在攻丝后，所加工出的ＳＴ或ＳＭ内螺纹的中径如果出现了超差或有特殊的要求时，用改变Ｂ值的方法进行螺纹的修正（使其满足标准内螺纹的要求或其它特殊要求），是一种最为方便、快速、有效、经济的方法。当然，其修正的范围是以满足Ｂ值在相关标准
规定的公差带内为前提的。
４．２Ｂ值和型面半角误差的要求
Ｂ（上限）≥Ｂ＋ｆ△α／２（４）Ｂ（下限）≤Ｂ（５）式中：Ｂ为实际尺寸（Ｂ本身数值）；ｆ△α／２为半角误差的Ｂ值补偿值。
另外，同一型面的左右两侧Ｂ值之差不大于０．０１ｍｍ。
４．３Ｂ值修正的应用举例
以在２０号钢上加工ＳＴ２７×２－５Ｈ螺孔为例，加工ＳＴ２７×２螺纹孔时其中径为２８．６８，该值已超出了ＳＴ２７×２螺纹孔标准中径（２８．２９９＋０．１５４
０
）的公差范围，当
装入Ｂ＝１．２９９－０．００７
－０．０１３
（标准规定Ｂ＝１．２９９０
－０．０１３
）的钢丝螺套后，标准内螺纹（Ｍ２７×２－５Ｈ）检验为不合格。现将钢丝螺套的Ｂ值进行修正，使Ｂ值处于其公差带的上偏差范围内（此时Ａｒ，Ｋ１，Ｋ２也会有相应的改变）。经装入后检查，该Ｂ值的钢丝螺套完全能通过标准内螺纹的通止规检验要求。
５结束语
上述研究与应用是通过了多次的研究、分析、试验所总结出来的，已成为了一种维修及改造的成熟经验。其在ＤＪＪ１型机车电机悬挂螺栓改造中得到了确认，并为今后的电力机车在特定条件下钢丝螺套的扩大应用提供了依据。
参考文献：
［１］石文祥．钢丝螺套的基本结构及安装使用［ＥＢ／ＯＬ］．ｈｔｔｐ：／／ｗｗｗ．ｍｗ３５．ｃｏｍ，２００５－０６－０１．
［２］陈小贤，陈广有，刘剑霄．普通型钢丝螺套技术条件［Ｚ］．北京：中国航空工业总公司第三○一研究所，
１９９７．
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!（上接第２４页）
－２７－
数值分析实验报告1
实验一误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 其中ε（1.1）和（1.221,,,a a 的输出b ”和“poly ε。 （1（2 （3）写成展 关于α solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 实验过程： 程序： a=poly(1:20); rr=roots(a); forn=2:21 n form=1:9 ess=10^(-6-m);
ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n
很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号：06450210 姓名：万轩 实验二插值法
数值分析实验报告176453
实验报告 插值法 数学实验室 数值逼近 算法设计 级 ____________________________ 号 ____________________________ 名 _____________________________ 实验项目名称 实验室 所属课程名称 实验类型 实验日期
实验概述： 【实验目的及要求】 本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容，掌握三种插 多项式插值，三次样条插值，拉格朗日插值，并比较三种插值方法的 优劣。 本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码，并 去实现。 【实验原理】 《数值分析》第二章“插值法”的相关内容，包括：牛顿多项式插值，三次样条插值， 拉格朗日 插值的相应算法和相关性质。 【实验环境】（使用的软硬件） 软件： MATLAB 2012a 硬件： 电脑型号：联想 Lenovo 昭阳E46A 笔记本电脑 操作系统： Win dows 8专业版 处理器：In tel ( R Core ( TM i3 CPU M 350 @2.27GHz 2.27GHz 实验内容： 【实验方案设计】 第一步，将书上关于三种插值方法的内容转化成程序语言，用 MATLA B 现; 第二步，分别用牛顿多项式插值，三次样条插值，拉格朗日插值求解不同的问题。 【实验过程】（实验步骤、记录、数据、分析） 实验的主要步骤是：首先分析问题，根据分析设计 MATLA 程序，利用程序算出 问题答案，分析所得答案结果，再得出最后结论。 实验一： 已知函数在下列各点的值为 试用4次牛顿插值多项式 P 4（ x ）及三次样条函数 S （ x ）（自然边界条件）对数据进行插值。 用图给出{（ X i , y i ）, X i =0.2+0.08i , i=0 , 1, 11, 10 } , P 4 （ x ）及 S （ x ）。 值方法：牛顿 在MATLAB 件中
数值分析实验报告
数值分析实验报告 姓名：周茹 学号： 912113850115 专业：数学与应用数学 指导老师：李建良
线性方程组的数值实验 一、课题名字：求解双对角线性方程组 二、问题描述 考虑一种特殊的对角线元素不为零的双对角线性方程组（以n=7为例） ?????????? ?????? ? ???? ?d a d a d a d a d a d a d 766 55 44 3 32 211??????????????????????x x x x x x x 7654321=?????????? ? ???????????b b b b b b b 7654321 写出一般的n （奇数）阶方程组程序（不要用消元法，因为不用它可以十分方便的解出这个方程组） 。 三、摘要 本文提出解三对角矩阵的一种十分简便的方法——追赶法，该算法适用于任意三对角方程组的求解。 四、引言 对于一般给定的d Ax =，我们可以用高斯消去法求解。但是高斯消去法过程复杂繁琐。对于特殊的三对角矩阵，如果A 是不可约的弱对角占优矩阵，可以将A 分解为UL ，再运用追赶法求解。
五、计算公式（数学模型） 对于形如????? ?? ????? ??? ?---b a c b a c b a c b n n n n n 111 2 2 2 11... ... ...的三对角矩阵UL A =，容易验证U 、L 具有如下形式： ??????? ????? ??? ?=u a u a u a u n n U ...... 3 3 22 1 ， ?? ????? ? ?? ??????=1 (1) 1132 1l l l L 比较UL A =两边元素，可以得到 ? ?? ??-== = l a b u u c l b u i i i i i i 111 i=2, 3, ... ,n 考虑三对角线系数矩阵的线性方程组 f Ax = 这里()T n x x x x ... 2 1 = ,()T n f f f f ... 2 1 = 令y Lx =，则有 f Uy = 于是有 ()?????-== --u y a f y u f y i i i i i 1 1 11 1 * i=2， 3, ... ,n 再根据y Lx =可得到
岩土工程计算原理和方法
岩土工程数值计算原理与方法 随着计算机的计算速度和存储能力的飞速发展以及计算方法的日益完善，数值模拟方法已经成为研究未知领域的强有力的工具。在岩土工程计算与分析中数值计算原理与方法也发展很快。特别是有限元的发展，促进了岩土工程研究、工程预测、优化设计和计算机辅助设计等的发展。但在工程实际中使用数值计算原理与方法却存在一些问题：例如有些人因缺乏对有限元和工程性质的深入了解，而有限元的迅速发展给他们造成一种假象，认为它是万能的，可以处理几乎所有的岩土工程问题；同时他们又被有限元计算结果的精度所迷惑，不了解这些精确结果后面所隐藏的不确定性，也不了解这些数值方法所采用本构模型的局限性以及相应参数的不确定性；因这些不确定性导致数值计算原理与方法的预测结果与实际情况和实际经验相差很大，又由于部分人计算偏于保守，使得岩土工程师难以接受现代数值计算原理与方法。 1. 岩土工程数值计算原理与方法也具有两面性。 有些人偏向于用其进行岩土工程的分析计算的原因在于： (1)数值计算原理与方法能够做任何传统的分析方法所能做到的分析与计算，而且做得更多、更好。 (2)数值计算原理与方法能够给出复杂数学模型的解。因而能够从机理上预测工程性质，而不是统计和经验性的描述，这是一大优点；而简化或经验分析方法有时只能描述其表面或形式上(统计)的关系，缺乏物理机制的描述和探讨。 (3)该方法既能处理简单问题，也能处理复杂问题。 数值计算原理与方法难以被其他人接受的原因在于： (1)使用复杂，难以被很好的掌握。 (2)数值计算原理与方法本身的不确定性(指与精确的解析方法相比所产生的不确定性，特别是在岩土动力非线性问题中这种不确定性会很大)导致预测结果与工程实际不符。 (3)数值计算原理与方法所使用的物理模型或本构模型有局限性，难以反映实际情况，导致预测结果与工程实际不符。 (4)采用复杂模型要求较多的参数，而这些参数难以用简单试验获得。 (5)既然数值计算原理与方法和传统的分析方法都具有很大的不确定性，还不如采用传统的分析方法，因为传统的方法简单、实用。 (6)精确的数值分析结果会误导使用者迷信这些结果的精确性，而没有认识到其后面隐
数值计算实验报告
(此文档为word格式，下载后您可任意编辑修改！) 2012级6班###(学号)计算机数值方法 实验报告成绩册 姓名:宋元台 学号： 成绩：
数值计算方法与算法实验报告 学期： 2014 至 2015 第 1 学期 2014年 12月1日课程名称: 数值计算方法与算法专业:信息与计算科学班级 12级5班 实验编号： 1实验项目Neton插值多项式指导教师：孙峪怀 姓名：宋元台学号：实验成绩： 一、实验目的及要求 实验目的： 掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。 实验要求： 1. 给出Newton插值算法 2. 用C语言实现算法 二、实验内容 三、实验步骤(该部分不够填写.请填写附页)
1.算法分析： 下面用伪码描述Newton插值多项式的算法： Step1 输入插值节点数n，插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)* Newton=newton+temp*f(i); Step5 输出f（x）的近似数值newton(x)=newton. 2.用C语言实现算法的程序代码 #include
数值分析实验报告模板
数值分析实验报告模板 篇一：数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二：数值分析实验报告 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远，Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收
敛，但精度不够。熟悉Matlab语言编程，学习编程要点。体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x);
数值分析实验报告1
实验一 误差分析 实验（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b =
的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots + 上述简单的Matlab 程序便得到（）的全部根，程序中的“ess ”即是（）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。 如果扰动项的系数ε很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何 （2）将方程（）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象 出现 （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将 方程（）写成展开的形式， ) 3.1(0 ),(1920=+-= x x x p αα 同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系为什么你发现了什么现象，哪些根关于α的变化更敏感 思考题一：（上述实验的改进） 在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。
数值分析2016上机实验报告
序言 数值分析是计算数学的范畴，有时也称它为计算数学、计算方法、数值方法等，其研究对象是各种数学问题的数值方法的设计、分析及其有关的数学理论和具体实现的一门学科，它是一个数学分支。是科学与工程计算（科学计算）的理论支持。许多科学与工程实际问题（核武器的研制、导弹的发射、气象预报）的解决都离不开科学计算。目前，试验、理论、计算已成为人类进行科学活动的三大方法。 数值分析是计算数学的一个主要部分,计算数学是数学科学的一个分支,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现。现在面向数值分析问题的计算机软件有：C,C++,MATLAB,Python,Fortran等。 MATLAB是matrix laboratory的英文缩写，它是由美国Mathwork公司于1967年推出的适合用于不同规格计算机和各种操纵系统的数学软件包，现已发展成为一种功能强大的计算机语言，特别适合用于科学和工程计算。目前，MATLAB应用非常广泛，主要用于算法开发、数据可视化、数值计算和数据分析等，除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。 本实验报告使用了MATLAB软件。对不动点迭代，函数逼近（lagrange插值，三次样条插值，最小二乘拟合），追赶法求解矩阵的解，4RungeKutta方法求解，欧拉法及改进欧拉法等算法做了简单的计算模拟实践。并比较了各种算法的优劣性，得到了对数值分析这们学科良好的理解，对以后的科研数值分析能力有了极大的提高。
目录 序言 (1) 问题一非线性方程数值解法 (3) 1.1 计算题目 (3) 1.2 迭代法分析 (3) 1.3计算结果分析及结论 (4) 问题二追赶法解三对角矩阵 (5) 2.1 问题 (5) 2.2 问题分析（追赶法） (6) 2.3 计算结果 (7) 问题三函数拟合 (7) 3.1 计算题目 (7) 3.2 题目分析 (7) 3.3 结果比较 (12) 问题四欧拉法解微分方程 (14) 4.1 计算题目 (14) 4.2.1 方程的准确解 (14) 4.2.2 Euler方法求解 (14) 4.2.3改进欧拉方法 (16) 问题五四阶龙格-库塔计算常微分方程初值问题 (17) 5.1 计算题目 (17) 5.2 四阶龙格-库塔方法分析 (18) 5.3 程序流程图 (18) 5.4 标准四阶Runge-Kutta法Matlab实现 (19) 5.5 计算结果及比较 (20) 问题六舍入误差观察 (22) 6.1 计算题目 (22) 6.2 计算结果 (22) 6.3 结论 (23) 7 总结 (24) 附录
数值分析实验报告
学生实验报告实验课程名称 开课实验室 学院年级专业班 学生姓名学号 开课时间至学年学期
if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位 （2）建立MATLAB界面 利用MA TLAB的GUI建立如下界面求解线性方程组： 详见程序。 五、计算实例、数据、结果、分析 下面我们对以上的结果进行测试，求解：
? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和，得到如下结果： 更改以上数据进行测试，求解如下方程组： 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果：
岩土工程数值分析学习笔记(DOC)
岩土工程数值分析读书笔记 摘要：阅读笔记分为两部分：理论学习和plaxis模拟相关问题。 理论部分 0岩土工程数值分析简介 岩土工程问题解析分析是以弹塑性力学理论和结构力学作为理论依据，适用于解决连续介质、各向同性材料、未知量少、边界条件简单的工程问题，存在很大的局限性。 岩土工程问题数值分析是借助于计算机的计算能力，适用于解决材料复杂、边界条件复杂、任意荷载、任意几何形状，适用范围广。 岩土工程数值分析发展过程： 20世纪40年代，使用差分法解决了土工中的渗流及固结问题，如土坝渗流及浸润线的求法、土坝及地基的固结等。 20世纪60年代，使用有限元法成解决了土石坝的静力问题的求解。 20世纪70年代，使用有限元法解决了土石坝及高楼(包括地基)的抗震分析。 20世纪80年代，边界元法异军突起，解决了半无限域的边界问题；地基的静力及动力问题都使用边界元法得到了有效地解决。 岩土工程数值分析的方法有两类，一类方法是将土视为连续介质，随后又将其离散化，如有限单元法、有限差分法、边界单元法、有限元线法、无单元法以及各种方法的耦合。另一类计算方法是考虑岩土材料本身的不连续性，如裂缝及不同材料间界面的界面模型和界面单元的使用，离散元法，不连续变形分析，流形元法，颗粒流等数值计算方法。 1数值分析过程中存在的问题及解决措施 问题：(1)对岩土工程数值分析方法缺乏系统的知识和深入的理解，出现问题时不知道在什么情况下属于理论问题或数学模型问题；在什么情况下是属于计算方法问题或本构模型问题；在什么情况下是参数的确定问题或计算本身的问题等。 (2)各种本构模型固有的局限性。具有多相性土的物理力学性质太复杂，难以准确地用数学模型和本构模型描述。例如邓肯一张模型不能反映剪胀性，不能反映压缩与剪切的交叉影响； (3)现有的试验手段和设备不能提供适当、合理和精确的参数。靠少数样本点所获得的参数难以准确地描述整个空间场地的物理力学性能；土的参数因土样扰动难以高质量的获取，其精度很差。 (4)数学模型还会给人造成一种错觉，让人觉得其计算结果也一定会更好、更可靠。这样可能使人们忽略了精确的数学公式也照样会有出错的可能性。只有当输入参数的质量和精度很高，并能与数学模型的精度相匹配时，才有可能得到较为准确的计算结果。 措施：(1)加强对土的本构模型的教学与培训，了解和掌握各种土的本构模型的优点和局限性以及模型参数的离散性。 (2)在使用数值分析方法的同时，不断地积累使用经验，包括他人的经验。
(完整版)哈工大-数值分析上机实验报告
实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。
Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0;
%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);
数值分析实验报告资料
机电工程学院 机械工程 陈星星 6720150109 《数值分析》课程设计实验报告 实验一 函数插值方法 一、问题提出 对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。 数据如下： （1 求五次Lagrange 多项式5L ()x ，计算(0.596)f ,(0.99)f 的值。（提示：结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈） 实验步骤： 第一步：先在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]; >>p = lagran(x,y) 回车得到： P = 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.6257 得到f（0.596）=0.6257 同理得到f（0.99）=1.0542
数值分析实验报告
实验报告 实验项目名称 估计水塔的水流量 实验项目类型 验证 演示 综合 设计 其他 指导教师 成绩 一、实验目的 （1）学会对实际问题的分析方法 （2）学会利用所学的知识解决实际问题 （3）设计出相应的算法，编制相应的应用程序 二、实验内容 某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m ，塔的直径为17.4m ，水塔是由水泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降低至最低水位，约8.2m 时，水泵自动启动加水。当水位升至最高水位，约10.8m 时，水泵停止工作。 下表给出了某一天的测量记录，测量了28个时刻的数据，但由于水泵正向水塔供水，由3个时刻无法测量到水位（表中为—）。 试建立数学模型，计算居民的用水速度和日总用水量。 三、实验原理、方法（算法）、步骤 时刻 0 0.921 1.843 2.949 3.871 4.978 5.900 水位 9.677 9.479 9.308 9.125 8.982 8.814. 8.686 时刻 7.006 7.928 8.967 9.981 10.925 10.945 12.032 水位 8.525 8.388 8.220 — — 10.820 10.500 时刻 12.954 13.875 14.982 15.903 1 6.826 1 7.931 19.037 水位 10.210 9.936 9.653 9.409 9.180 8.921 8.662 时刻 1 9.959 20.839 22.015 22.958 23.880 24.986 25.908 水位 8.433 8.220 — 10.820 10.591 10.354 10.180
关于岩土工程的数值计算方法的综述
题目：关于岩土工程的数值计算方法的综述学院：资源与土木工程学院 专业：岩土工程 学号： 姓名：
关于岩土工程的数值计算方法的综述 我通过学习和查阅相关资料文献了解到，近年来，数值计算模拟分析在岩土工程中越来越受欢迎，随着城市的建设，地下工程所处的环境越来越复杂，影响的因素也是越来越多，所以依靠传统的解析计算难以实现，计算机的数值模拟恰恰解决的了岩土的计算的问题，它可以模拟各种复杂情况下岩土问题。就岩土工程而言，由于岩土介质涉及本构关系、力学参数、自身构造以及边界条件等的复杂多变性，在未采用计算机数值方法以前，对于复杂、重要的岩土工程，如果用传统的弹性力学或弹塑性力学的解析法难以求解时，只好采用物理模拟或其他方法从宏观上把握工程的受力和变形特征。随着计算机数值分析方法的出现和发展，情况发生了巨大的变化。计算机数值方法已经能够较好的模拟非均匀质体、各向异性介质面临的复杂边界条件问题，也能处理岩土工程中不连续性界面、渗流问题、岩土损伤断裂问题以及复杂的岩土工程结构分析问题，对于涉及时间因素的动力问题、蠕变问题，特别是耦合问题，数值模拟计算方法极大的加强了解决岩土工程的能力。 数值计算方法其主要有有限单元法、有限差分法、边界元法、离散元法和流形元法等。 有限单元法：有限单元法发展非常迅速，至今已经成为求解复杂工程问题的有力工具，并在岩土工程领域广泛的采用，主要的分析软件ANSYS。 有限单元法的最基本的元素是单元和节点，基本计算步骤的第一步为离散化，问题域的连续体被离散为单元与节点的组合，连续体内部分的应力及位移通过节点传递，每个单元可以具有不同的物理特征，这样，便可以得到在物理意义上与原来的连续体相近似的模型。第二步为单元分析，一般以位移法为基本方法，建立单元的刚度矩阵。第三步由单元的刚度矩阵集合成总体刚度矩阵，并由此建立系统的整体方程组。第四步进入计算模型的边界条件，求解方程组，求得节点位移。第五步求出各单元的应变、应力及主应力。 有限差分法：有限差分法在岩土工程中是应用非常广泛的方法，在数值计算模拟上有很大的贡献，主要的应用软件为FLAC3D。 基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方
数值分析实验报告1
实验一 误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（1.1）中19x 的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve =
车辆动力学仿真
车辆动力学仿真 课程编码：202060 课程英文译名：Dynamics Simulation of Vehicle System 课程类别：专业课 开课对象：车辆工程专业开课学期：第7学期 学分：2.5学分；总学时： 40学时；理论课学时：32学时；上机学时： 8学时 先修课程：理论力学、材料力学、机械原理、机械设计、机械振动 教材：车辆动力学模拟及其方法，威鲁麦特（德），北京理工大学出版社, 1998.5 ，第1版 参考书：【1】汽车系统动力学，张洪欣，同济大学出版社, 1996 ，第1版【2】汽车系统动力学及仿真，雷雨成，国防工业出版社, 1997 ，第1版一、课程的性质、目的和任务 《车辆系统动力学仿真》是车辆工程专业理论性较强的专业课。本课程的目的是，使学生初步学会汽车动力学分析方法，能够解决工程实际问题，以便增强其研究和解决车辆动力学问题的能力。本课程的任务，是以数学力学模型为基础，结合虚拟样机仿真技术，讲授汽车的垂直动力学、横向动力学、纵向动力学，为继续学习和掌握汽车新科技创造条件。 二、课程的基本要求 对汽车动力学有一定的了解，掌握有关的基本概念、基本理论和基本方法及其应用，掌握汽车多体动力学仿真的方法。具体要求为： 1.对汽车动力学仿真的基本概念和基本分析方法有明确的认识； 2.掌握单自由度系统的振动系统，自由振动、强迫振动的微分方程的建立方法； 3.掌握多自由度系统的振动系统的微分方程，初步掌握多自由度系统振动的模态分析方法； 4.了解随机振动的一些基本概念，掌握路面不平度功率谱密度的概念及其计算方法； 5.掌握汽车垂直动力学模型的建立方法，以及路面激励对汽车振动的影响； 6.掌握汽车弹簧、减震器、橡胶金属部件、轮胎等部件垂向动力学的特性； 7.掌握汽车纵向动力学微分方程，掌握滚动阻力、爬坡阻力、加速阻力的计算方法； 8.掌握驱动附着率、制动附着率对行驶极限的影响； 9.掌握汽车横向动力学的微分方程建立方法，及其横向动力学微分方程的特性； 10.掌握汽车操作稳定性的概念及其影响汽车操作稳定性的因素； 11.掌握轮胎的真实特性，初步掌握轮胎动力学的初步概念。
数值分析实验报告
实验五 解线性方程组的直接方法 实验5.1 （主元的选取与算法的稳定性） 问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。 实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。 实验要求： （1）取矩阵?? ? ?? ?? ?????????=????????????????=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。取n=10计算矩阵的 条件数。让程序自动选取主元，结果如何？ （2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。 （3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。 思考题一：（Vadermonde 矩阵）设 ?? ??????????????????????=? ? ? ?????????????=∑∑∑∑====n i i n n i i n i i n i i n n n n n n n x x x x b x x x x x x x x x x x x A 0020 10022222121102001111 ，， 其中，n k k x k ,,1,0,1.01 =+=， （1）对n=2,5,8，计算A 的条件数；随n 增大，矩阵性态如何变化？ （2）对n=5，解方程组Ax=b ；设A 的最后一个元素有扰动10-4，再求解Ax=b （3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。 （4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？ 相关MATLAB 函数提示： zeros(m,n) 生成m 行，n 列的零矩阵 ones(m,n) 生成m 行，n 列的元素全为1的矩阵 eye(n) 生成n 阶单位矩阵 rand(m,n) 生成m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x 的元素构成的对角矩阵 tril(A) 提取矩阵A 的下三角部分生成下三角矩阵
数值分析实验报告2
一、实验名称 复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式及自适应辛普森积分。 二、实验目的及要求 1. 掌握复合梯形求积计算积分、复合辛普森求积计算积分、龙贝格求积计算积分和自适应辛普森积分的基本思路和步骤. 2. 培养Matlab 编程与上机调试能力. 三、实验环境 计算机，MATLAB 软件 四、实验内容 1.用不同数值方法计算积分9 4 ln 1 0-=? xdx x 。 （1）取不同的步长h 。分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确指比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善。 （2）用龙贝格求积计算完成问题（1）。 （3）用自适应辛普森积分，使其精度达到10-4。 五、算法描述及实验步骤 1.复合梯形公式 将区间[a,b]划分为n 等份，分点x k =a+ah,h=(b-a)/h,k=0,1,...,n ，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用梯形公式（），得 )]()([2 )(b f a f a b dx x f b a +-≈ ? （） )]()(2)([2)]()([21 1 110b f x f b f h x f x f h T n k k k n k k n ++=+=∑∑-=+-= （） ),(),(12 )(' '2b a f h a b f R n ∈-- =ηη （） 其中Tn 称为复合梯形公式，Rn 为复合梯形公式的余项。 2.复合辛普森求积公式 将区间[a,b]划分为n 等份，在每个子区间[x k ,x k +1](k=0,1,...,n-1)上采用辛普森公式（），得 )]()2 (4)([6b f b a f a f a b S +++-= （）
岩土工程数值分析读书报告
岩土工程数值分析读书报告 一.岩土与数值分析 在很多岩土工程的实际问题中，例如档土墙、板桩、基础梁和板等工程，由于岩土的非均质、非线性的性状以及几何形状的任意性、不连续性等因素，在多数情况下不能获得解析解。最近二十多年来，随着电子计算机的迅速兴起，在岩土工程中，数值分析受到了极大的重视，各种数值方法在岩土工程中都得到了广泛地应用，而岩土工程中的各种复杂问题的解决又深化和丰富了数值分析的内容。 目前．在岩土工程的数值分析中，用的最为普遍的是有限元法和差分法，其他方法如边界元法正在兴起。变分法与加权余量法既可以独立地作为数值方法运用于土工实际问题的求解，又可作为推导前几种数值方法的手段。当数值分析中的差分法首先盛行于工程科学时，土工中的渗流及固结问题在四十年代后期也开始采用差分法成功地解决了某些实际问题，如土坝渗流及浸润线的求法、土坝及地基的固结等。五十年代及六十年代初，弹性地基上的梁与板以及板桩也用差分法来求解。六十年代，土石坝的静力问题用有限元法来求解。由于有限元解法的灵活性，使差分法在土工中的应用暂时趋丁停滞。进入七十年代之后，土石坝及高楼（包括地基）成功地使用有限无法解决了抗震分析。七十午代后期及八十年代，边界元法异军突起。这方法特别适宜于半无限域课题，这些是土力学及地基工程学科经常遇到的边界情况。近十年来，地基的静力及动力问题，例如桩基及强夯（即
动力固结）等，都使用边界元法得到了有效地解决。 岩土工程数值分析的方法有两类，一类方法是将土视为连续介质，随后又将其离散化，如有限单元法、有限差分法、边界单元法、有限元线法、无单元法以及各种方法的耦合。另一类计算方法是考虑岩土材料本身的不连续性，如裂缝及不同材料间界面的界面模型和界面单元的使用，离散元法（DEM），不连续变形分析（DDA），流形元法 （MEM），颗粒流（PFC）等数值计算方法迅速发展。 二.土的本构关系 材料的本构关系（constitutive relationship）是反映材料的力学性状的数学表达式，表示形式一般为应力-应变-时间的关系，也称为本构定律（constitutive law）、本构方程（constitutive equation）,还可称为本构关系数学模型（mathematical model）简称为本构模型。 （一）土的弹性模型 在线弹性模型中，假定材料符合弹性力学规律，应力-应变关系式为： { } = [D]{ } 这里刚度矩阵称为弹性矩阵，由广义虎克定律 L，"~，1■- 1 J —V 式中包含了弹性模量和泊松比柑两个常数。。它们可以用另外两 个弹性常数，剪切模量G和体积模量K来代替。它们之间的关系为
数值分析实验报告3
实验报告 实验项目名称数值积分与数值微分实验室数学实验室 所属课程名称数值逼近 实验类型算法设计 实验日期 班级 学号 姓名 成绩
实验概述： 【实验目的及要求】 本次实验的目的是熟练《数值分析》第四章“数值积分与数值微分”的相关内容，掌握复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式。 本次试验要求编写复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的程序编码，并在MATLAB软件中去实现。 【实验原理】 《数值分析》第四章“数值积分与数值微分”的相关内容，包括：复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的相应算法和相关性质。 【实验环境】（使用的软硬件） 软件： MATLAB 2012a 硬件： 电脑型号：联想 Lenovo 昭阳E46A笔记本电脑 操作系统：Windows 8 专业版 处理器：Intel（R）Core（TM）i3 CPU M 350 @2.27GHz 2.27GHz 实验内容： 【实验方案设计】 第一步，将书上关于复合梯形求积公式、复合辛普森求积公式、龙贝格求积公式以及高斯-勒让德公式的内容转化成程序语言，用MATLAB实现；第二步，分别用以上求积公式的程序编码求解不同的问题。 【实验过程】（实验步骤、记录、数据、分析） 实验的主要步骤是：首先分析问题，根据分析设计MATLAB程序，利用程序算出问题答案，分析所得答案结果，再得出最后结论。 实验：用不同数值方法计算积分 (1) 取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？ (2) 用龙贝格求积计算完成问题（1）。 (3)用勒让德多项式确定零点，再代入计算高斯公式，使其精度达到10-4 （1）在MATLAB的Editor中建立一个M-文件，输入程序代码，实现复合梯形求积公式的程序代码如下：