下载文档原格式
实验报告课程名称：计算方法院系：数学科学系专业班级：数应1001学号：1031110139学生姓名：姚海保指导教师：沈林开课时间：2012至2013学年第一学期一、学生撰写要求按照实验课程培养方案的要求，每门实验课程中的每一个实验项目完成后，每位参加实验的学生均须在实验教师规定的时间内独立完成一份实验报告，不得抄袭，不得缺交。
学生撰写实验报告时应严格按照本实验报告规定的内容和要求填写。
字迹工整，文字简练，数据齐全，图表规范，计算正确，分析充分、具体、定量。
二、教师评阅与装订要求1.实验报告批改要深入细致，批改过程中要发现和纠正学生实验报告中的问题，给出评语和实验报告成绩，签名并注明批改日期。
实验报告批改完成后，应采用适当的形式将学生实验报告中存在的问题及时反馈给学生。
2.实验报告成绩用百分制评定，并给出成绩评定的依据或评分标准（附于实验报告成绩登记表后）。
对迟交实验报告的学生要酌情扣分，对缺交和抄袭实验报告的学生应及时批评教育，并对该次实验报告的分数以零分处理。
对单独设课的实验课程，如学生抄袭或缺交实验报告达该课程全学期实验报告总次数三分之一以上，不得同意其参加本课程的考核。
3.各实验项目的实验报告成绩登记在实验报告成绩登记表中。
本学期实验项目全部完成后，给定实验报告综合成绩。
4.实验报告综合成绩应按课程教学大纲规定比例（一般为10-15%）计入实验课总评成绩；实验总评成绩原则上应包括考勤、实验报告、考核（操作、理论）等多方面成绩；5.实验教师每学期负责对拟存档的学生实验报告按课程、学生收齐并装订，按如下顺序装订成册：实验报告封面、实验报告成绩登记表、实验报告成绩评定依据、实验报告（按教学进度表规定的实验项目顺序排序）。
装订时统一靠左侧按“两钉三等分”原则装订。
024********-1-0.8-0.6-0.4-0.20.20.40.60.813、画出2222)sin(yxyxz++=所表示的三维曲面。
现浇钢筋混凝土空心板自然频率与振型的差分法研究的开题报告一、选题背景现浇钢筋混凝土空心板是一种常用于建筑工程中的楼板结构形式，具有良好的承载能力和隔音性能。
然而，随着建筑物的高度和地震等自然灾害频率的增加，现浇钢筋混凝土空心板的抗震性能成为了一个备受关注的问题。
因此，研究现浇钢筋混凝土空心板的自然频率和振型对于提升该结构的抗震性能具有重要意义。
二、研究内容本研究旨在通过差分法研究现浇钢筋混凝土空心板的自然频率和振型。
首先，根据空心板的几何尺寸、构造形式和钢筋混凝土的材料力学性质，建立数学模型，并通过有限元软件（如ANSYS）进行建模和数值仿真。
然后，利用差分法计算空心板的自然频率和振型。
最后，通过对仿真结果进行分析和比较，得出了空心板在不同荷载情况下的振动特性。
三、研究意义本研究的结果可以帮助工程师更好地理解现浇钢筋混凝土空心板的振动特性，进而设计出更加安全可靠的结构。
此外，本研究的方法也可以为其他钢筋混凝土结构的振动分析提供参考。
四、研究方法本研究采用差分法计算空心板的自然频率和振型。
具体来说，差分法是一种将微分方程转化为代数方程的方法。
通过将二阶微分方程进行差分，可以得到空心板的特征方程和振型方程。
然后，利用MATLAB等工具进行数值解算和结果可视化。
五、研究计划本研究计划分为以下几个阶段：1. 研究文献综述，了解现浇钢筋混凝土空心板的特点和研究动态；2. 建立空心板的有限元模型，并进行数值仿真；3. 利用差分法计算空心板的自然频率和振型；4. 分析仿真结果，得出空心板在不同荷载情况下的振动特性；5. 编写论文并进行答辩。
六、预期成果本研究预期可以得出现浇钢筋混凝土空心板在不同荷载情况下的自然频率和振型，并分析振动特性，为该结构的抗震设计提供依据和参考；同时，本研究的方法和分析结果也可适用于其他钢筋混凝土结构的振动分析。
结构力学程序大作业报告一．题目要求要求编写连续梁内力计算程序。
二．程序总框图三．程序变量声明NE——单元总数N——节点转角未知量总数BL(NE)——单元杆长数组2EI(NE)——单元抗弯刚度数组，后存单元线刚度P(N)——等效节点荷载数组，后存节点转角KE(N,2)——整体刚度矩阵数组JD(NE,2)——单元定位向量数组FM(NE,2)——单元固端弯矩数组FJ(NE,2)——单元杆端弯矩数组四．程序源代码见附件，或者所附的Matlab文件。
五．程序算例例题一：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[1,2;2,3;3,4;4,5;5,6;6,0]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：1 232 33 44 55 66 0单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.0000-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,0,-8,0,10,0]节点转角未知量总数： 6直接节点力矩： 0 0 -8 0 10 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载：10.6670 -1.6670 -11.0000 15.3330 -8.3330 15.0000单元节点位移：11.3842 -1.4345 -8.9803 14.0534 -10.1918 10.0480单元杆端弯矩0 14.9246-14.9246 -0.6976-7.3024 12.3755-12.3755 18.167945-8.1679 7.9520 -7.9520 23.0240做弯矩图：例题一结构弯矩图例题二：计算结果： 输入单元总数：6输入单元杆长：[4,6,6,8,4,6] 输入单元抗弯刚度：[1,1.5,1,2,1,1.5] 输入单元定位向量:[0,1;1,2;2,3;3,4;4,5;5,6]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18] 单元数： 67单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：0 11 22 33 44 55 6单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.0000-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,-8,0,10,0,0]节点转角未知量总数： 6直接节点力矩： 0 -8 0 10 0 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载：-1.6670 -11.0000 15.3330 -8.3330 15.0000 -18.0000单元节点位移：1.6865 -10.0801 14.8707 -12.1830 17.1951 -26.5976单元杆端弯矩6-9.8237 12.3535-12.3535 -0.2368-7.7632 12.5538-12.5538 16.5854-6.5854 14.1037-14.1037 0弯矩图：例题二结构弯矩图例题三：计算结果：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]7输入单元定位向量:[0,1;1,2;2,3;3,4;4,5;5,0]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：0 11 22 33 44 55 0单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.0000-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,-8,0,10,0]节点转角未知量总数： 5直接节点力矩： 0 -8 0 10 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载：-1.6670 -11.0000 15.3330 -8.3330 15.00008单元节点位移：1.6537 -9.9489 14.2640 -10.2480 10.0620 单元杆端弯矩-9.8401 12.3207-12.3207 -0.1221-7.8779 12.1930-12.1930 18.2170-8.2170 7.9380-7.9380 23.0310弯矩图：例题三机构弯矩图例题四：计算结果：输入单元总数：69输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[1,2;2,3;3,4;4,5;5,6;6,7]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：1 22 33 44 55 66 7单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.0000-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,0,-8,0,10,0,0]节点转角未知量总数： 7直接节点力矩： 0 0 -8 0 10 0 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.250010折算固端弯矩后的节点荷载：10.6670 -1.6670 -11.0000 15.3330 -8.3330 15.0000 -18.0000 单元节点位移：11.3653 -1.3965 -9.1131 14.6603 -12.1263 17.1789 -26.5895 单元杆端弯矩0 14.9531-14.9531 -0.8114-7.1886 12.7358-12.7358 16.5368-6.5368 14.1158-14.1158 0弯矩图：第四题结构弯矩图例题五：(工况情况1 and 工况情况3)11计算结果：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[1,2;2,3;3,4;4,5;5,6;6,7]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：1 22 33 44 55 66 7单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.0000-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,0,0,0,0,0,0]节点转角未知量总数： 712直接节点力矩： 0 0 0 0 0 0 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载：10.6670 -1.6670 -3.0000 15.3330 -18.3330 15.0000 -18.0000单元节点位移：12.0951 -2.8562 -4.0044 15.3061 -17.6848 18.7671 -27.3835单元杆端弯矩0 13.8584-13.8584 3.5675-3.5675 14.8693-14.8693 11.3013-11.3013 12.9247-12.9247 0弯矩图：例题六：（荷载工况4）13计算结果：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[0,1;1,2;2,3;3,4;4,5;5,0]输入单元固端弯矩：[-10.667,10.667;0,0;0,0;-21.333,21.333;0,0;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：0 11 22 33 44 55 0单元固端弯矩：-10.6670 10.66700 00 0-21.3330 21.33300 0-18.0000 18.0000确认输入数据正确（Y/N）：Y输入直接节点力矩：[0,-8,0,10,0]节点转角未知量总数： 514直接节点力矩： 0 -8 0 10 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载：-10.6670 -8.0000 21.3330 -11.3330 18.0000单元节点位移：-3.4807 -7.4113 18.2775 -13.3183 12.3296单元杆端弯矩-12.4073 7.1863-7.1863 -9.15161.1516 9.7146-9.7146 17.1535-7.1535 5.6704-5.6704 24.1648弯矩图：例题六结构弯矩图例题七：（载荷工况2）15计算结果：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[0,1;1,2;2,3;3,4;4,5;5,6]输入单元固端弯矩：[0,0;-9,9;-6,6;0,0;-3,3;0,0]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：0 11 22 33 44 55 6单元固端弯矩：0 0-9 9-6 60 0-3 316170 0确认输入数据正确（Y/N ）：Y输入直接节点力矩：[0,-8,0,10,0,0]节点转角未知量总数： 6直接节点力矩： 0 -8 0 10 0 0节点线刚度矩阵： 0.2500 0.2500 0.1667 0.2500 0.2500 0.2500 折算固端弯矩后的节点荷载： 9 -11 -6 13 -3 0单元节点位移：6.3944 -7.5778 -4.70278.7277 -4.2079 2.1040单元杆端弯矩3.1972 6.3944-6.3944 4.6194-12.6194 0.3389-0.3389 6.37633.6237 3.1559-3.1559 0.0000弯矩图：例题七结构弯矩图例题八：（载荷工况5）计算结果：输入单元总数：6输入单元杆长：[4,6,6,8,4,6]输入单元抗弯刚度：[1,1.5,1,2,1,1.5]输入单元定位向量:[1,2;2,3;3,4;4,5;5,6;6,0]输入单元固端弯矩：[-10.667,10.667;-9,9;-6,6;-21.333,21.333;-3,3;-18,18]单元数： 6单元杆长： 4 6 6 8 4 6单元抗弯刚度： 1.0000 1.5000 1.0000 2.0000 1.0000 1.5000 单元定位向量：1 22 33 44 55 66 0单元固端弯矩：-10.6670 10.6670-9.0000 9.0000-6.0000 6.0000-21.3330 21.3330-3.0000 3.00001819-18.0000 18.0000确认输入数据正确（Y/N ）：Y输入直接节点力矩：[0,0,0,0,0,0]节点转角未知量总数： 6直接节点力矩： 0 0 0 0 0 0节点线刚度矩阵：0.2500 0.2500 0.1667 0.2500 0.2500 0.2500折算固端弯矩后的节点荷载：10.6670 -1.6670 -3.0000 15.3330 -18.3330 15.0000单元节点位移：12.1146 -2.8952 -3.8676 14.6811 -15.6926 11.4231单元杆端弯矩0 13.8291-13.8291 3.6847-3.6847 14.4982-14.4982 12.9810-12.9810 6.5769-6.5769 23.7116弯矩图：例题八结构弯矩图附件：程序源代码disp('输入单元总数：'); % NE表示单元总数，为1*n矩阵（设NE=n）NE=input('');disp('输入单元杆长：'); % BL表示单元杆长矩阵，为1*n矩阵BL=input('');disp('输入单元抗弯刚度：'); % EI表示单元抗弯刚度矩阵，为n*2矩阵EI=input('');disp('输入单元定位向量:'); % JD表示单元定位向量，为n*2矩阵JD=input('');disp('输入单元固端弯矩：'); % FM表示单元固端弯矩，为n*2矩阵FM=input('');disp('单元数：');disp(NE);disp('单元杆长：');disp(BL);disp('单元抗弯刚度：');disp(EI);disp('单元定位向量：');disp(JD);disp('单元固端弯矩：');disp(FM);disp('确认输入数据正确（Y/N）：'); % 确认所输入数据是否有误A1=input('','s');if A1=='N'disp('请重启程序，重新输入');pause(2);close;endN=JD(NE,1); % 确定N的维数if JD(NE,2)~=0N=JD(NE,2);Enddisp('输入直接节点力矩：')P=input(''); % P表示节点载荷，为1*N阶矩阵（N为未知变量数，n为单元数，有N≤n）20disp('节点转角未知量总数：');disp(N);disp('直接节点力矩：');disp(P);for I=1:NE % 得到单元线刚度矩阵EI EI(I)=EI(I)/BL(I);enddisp('节点线刚度矩阵：');disp(EI);for I=1:NE % 将固端弯矩折算到节点载荷中I1=JD(I,1);I2=JD(I,2);if I1~=0P(I1)=P(I1)-FM(I,1);endif I2~=0P(I2)=P(I2)-FM(I,2);endenddisp('折算固端弯矩后的节点荷载：');disp(P);for I=1:N % 在形成整体刚度矩阵前，将其初值赋为0 for J=1:2KE(I,J)=0.0;endendfor I=1:NE % 求解整体刚度矩阵，为n*2阶矩阵I1=JD(I,1);I2=JD(I,2);if I1~=0KE(I1,1)=KE(I1,1)+4*EI(I);if I2~=0KE(I1,2)=2*EI(I);endendif I2~=0KE(I2,1)=KE(I2,1)+4*EI(I);end21endfor I=1:N-1 % 高斯消元的过程C=KE(I,2)/KE(I,1);KE(I+1,1)=KE(I+1,1)-C*KE(I,2);P(I+1)=P(I+1)-C*P(I);endP(N)=P(N)/KE(N,1);for I=N-1:-1:1P(I)=P(I)-KE(I,2)*P(I+1);P(I)=P(I)/KE(I,1);Enddisp('单元节点位移：');disp(P); % 输出计算得出的节点位移for I=1:NE % 计算杆端弯矩I1=JD(I,1);I2=JD(I,2);DZ(1)=0;DZ(2)=0;if I1~=0DZ(1)=P(I1);endif I2~=0DZ(2)=P(I2);endFJ(I,1)=4*EI(I)*DZ(1)+2*EI(I)*DZ(2)+FM(I,1); FJ(I,2)=2*EI(I)*DZ(1)+4*EI(I)*DZ(2)+FM(I,2); enddisp('单元杆端弯矩');disp(FJ);22。
差分编程算法
差分编程算法是一种用于解决问题的常见算法思想，它通过计算相邻元素之间的差异来获取有用的信息。
差分算法常被用于处理序列数据，如时间序列、图像处理等领域。
下面我将通过一个具体的例子来说明差分编程算法的应用。
假设有一组表示每天气温的序列数据，我们想要找出连续五天中气温增加最快的那一天。
我们可以使用差分编程算法来解决这个问题。
我们需要计算相邻两天的气温差异。
假设气温序列数据为[20, 18, 22, 25, 28, 26, 30, 32]，我们可以得到差分序列[-2, 4, 3, 3, -2, 4, 2]。
然后，我们可以找出差分序列中的最大值所对应的索引，即找出气温增加最快的那一天。
在这个例子中，差分序列中的最大值为4，对应的索引为5。
因此，第六天的气温增加最快。
通过差分编程算法，我们可以快速找到气温增加最快的那一天。
这个算法的时间复杂度为O(n)，其中n为序列数据的长度。
差分编程算法在处理序列数据的增减趋势分析、峰值检测等问题上具有广泛的应用。
总结一下，差分编程算法是一种用于解决问题的常见算法思想，它通过计算相邻元素之间的差异来获取有用的信息。
在处理序列数据的增减趋势分析、峰值检测等问题上具有广泛的应用。
希望通过这
个例子能够让你对差分编程算法有一个初步的了解。
中心差分法的基本理论与程序设计
中心差分法（Central Difference Method）是一种常用的数值计算方法，用于近似求解函数的导数。
该方法基于导数的定义，通过两个函数值的差分来近似求解导数的值。
本文将介绍中心差分法的基本理论和程序设计过程。
一、中心差分法的基本理论
f'(x)≈(f(x+h)-f(x-h))/(2h)
其中，h是步长，决定了计算函数值的邻近点的间距。
二、中心差分法的程序设计
下面我们将介绍中心差分法的程序设计过程，包括步骤和代码实现。
步骤：
1.定义函数f(x)，表示需要求解导数的函数。
2.定义步长h，决定计算函数值的邻近点的间距。
3.根据中心差分法的公式，计算导数的近似值：
f'(x)≈(f(x+h)-f(x-h))/(2h)
4.返回导数的近似值。
代码实现：
```
def central_difference(f, x, h):
"""
计算函数f在点x处的导数值
:param f: 函数f(x)
:param x: 需要求导数的点
:param h: 步长
:return: 导数的近似值
"""
df = (f(x + h) - f(x - h)) / (2 * h)
return df
```
以上就是中心差分法的基本理论和程序设计过程。
中心差分法是一种简单有效的数值计算方法，可以用于求解导数的近似值。
在实际应用中，可以根据需要选择合适的步长，并通过增加点的数量来提高计算结果的精度。
分支结构程序设计实验报告实验目的本次实验旨在通过编写分支结构程序，加深对分支结构的理解，并提高编写程序的能力。
实验内容本次实验分为两个部分，分别为编写一个判断分数等级的程序和一个银行利率计算程序。
判断分数等级编写一个程序，根据输入的分数输出对应的等级。
等级规则如下：- 90分以上：优秀- 80-89分：良好- 70-79分：中等- 60-69分：及格- 60分以下：不及格银行利率计算编写一个程序，根据输入的存款金额和存款期限计算出存款到期后的本息合计。
假设银行利率为年利率，根据存款期限的不同，利率如下：- 存款期限小于等于1年：年利率为3%- 存款期限大于1年小于等于3年：年利率为4% - 存款期限大于3年：年利率为5%实验步骤判断分数等级程序思路1. 接收用户输入的分数。
2. 使用if-elif-else语句判断分数所处的等级范围。
3. 根据判断结果，输出对应的等级。
银行利率计算程序思路1. 接收用户输入的存款金额和存款期限。
2. 使用if-elif-else语句判断存款期限范围。
3. 根据判断结果，计算存款到期后的本息合计。
4. 输出存款到期后的本息合计。
实验代码判断分数等级程序代码pythonscore = float(input("请输入分数："))if score >= 90:print("优秀")elif score >= 80:print("良好")elif score >= 70:print("中等")elif score >= 60:print("及格")else:print("不及格")银行利率计算程序代码pythonamount = float(input("请输入存款金额：")) period = float(input("请输入存款期限："))if period <= 1:interest_rate = 0.03elif period <= 3:interest_rate = 0.04else:interest_rate = 0.05total_amount = amount * (1 + interest_rate)print("存款到期后的本息合计为：", total_amount)实验结果与分析判断分数等级程序当输入分数为85时，程序输出良好，符合预期。
实验报告课程名称：计算方法院系：数学科学系专业班级：数应1001学号：1031110139学生姓名：姚海保指导教师：沈林开课时间：2012至2013学年第一学期一、学生撰写要求按照实验课程培养方案的要求，每门实验课程中的每一个实验项目完成后，每位参加实验的学生均须在实验教师规定的时间内独立完成一份实验报告，不得抄袭，不得缺交。
学生撰写实验报告时应严格按照本实验报告规定的内容和要求填写。
字迹工整，文字简练，数据齐全，图表规范，计算正确，分析充分、具体、定量。
二、教师评阅与装订要求1.实验报告批改要深入细致，批改过程中要发现和纠正学生实验报告中的问题，给出评语和实验报告成绩，签名并注明批改日期。
实验报告批改完成后，应采用适当的形式将学生实验报告中存在的问题及时反馈给学生。
2.实验报告成绩用百分制评定，并给出成绩评定的依据或评分标准（附于实验报告成绩登记表后）。
对迟交实验报告的学生要酌情扣分，对缺交和抄袭实验报告的学生应及时批评教育，并对该次实验报告的分数以零分处理。
对单独设课的实验课程，如学生抄袭或缺交实验报告达该课程全学期实验报告总次数三分之一以上，不得同意其参加本课程的考核。
3.各实验项目的实验报告成绩登记在实验报告成绩登记表中。
本学期实验项目全部完成后，给定实验报告综合成绩。
4.实验报告综合成绩应按课程教学大纲规定比例（一般为10-15%）计入实验课总评成绩；实验总评成绩原则上应包括考勤、实验报告、考核（操作、理论）等多方面成绩；5.实验教师每学期负责对拟存档的学生实验报告按课程、学生收齐并装订，按如下顺序装订成册：实验报告封面、实验报告成绩登记表、实验报告成绩评定依据、实验报告（按教学进度表规定的实验项目顺序排序）。
装订时统一靠左侧按“两钉三等分”原则装订。
4、复数矩阵的生成及运算A=[1,3;2,4]-[5,8;6,9]*iB=[1+5i,2+6i;3+8*i,4+9*i]C=A*BA = 1.0000 - 5.0000i 3.0000 - 8.0000i2.0000 - 6.0000i 4.0000 - 9.0000iB =1.0000 + 5.0000i 2.0000 + 6.0000i3.0000 + 8.0000i4.0000 + 9.0000iC =1.0e+002 *0.9900 1.1600 - 0.0900i1.1600 + 0.0900i 1.3700。
中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明1.1 中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
1.2 中心差分法原理中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为：t u u u i i ∆-=-+•211 (a) 2112t u u u u i i i ∆+-=-+•• (b) 而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……)由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c) 将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程：02211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m （d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c t m u t m k u t c t m (e)由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
1.3 初始条件转化假设给定的初始条件为),0(),0(00••==u u u u （g ）由式（g ）确定1-u 。
在零时刻速度和加速度的中心差分公式为：t u u u ∆-=-•2110 （h ） ` 210102t u u u u ∆+-=-•• （i ）将式（i ）消去1u 得：020012•••-∆+∆-=u t u t u u （j ）而零时刻的加速度值0••u 可以用t ＝0时的运动方程0000=++•••ku u c u m确定 即 )(1000ku u c m u --=••• （k ） 这样就可以根据初始条件00,•u u 和初始荷载0P ，就可以根据上式确定1-u 的值。
偏微分中心差分格式实验报告实验目的：1.掌握偏微分的中心差分格式；2.理解中心差分格式的精度和稳定性。
实验原理：中心差分是一种常用的数值求解偏微分方程的格式，其基本思想是用函数在两个点的导数的平均值来近似函数在这两个点中间的导数值。
具体来说，对于一维的偏微分方程，中心差分格式可以表述为：f'(x)=(f(x+h)-f(x-h))/(2h)其中f'(x)表示x点处的导数，h表示步长。
实验步骤：1.编写一个计算函数在任意给定点x处的导数值的中心差分程序；2.给定一个函数f(x)，例如f(x)=x^2，计算在一定范围内的该函数在每个点处的导数值；3.比较计算的导数值与理论值的差异，并分析中心差分格式的精度；4.对给定步长h，逐渐减小h，计算导数值，并观察数值的变化，分析中心差分格式的稳定性。
实验结果与分析：以函数f(x)=x^2为例，给定步长h=0.1，计算在范围[-1,1]内的函数f(x)在每个点处的导数值。
实验结果如下表所示：x，f'(x)，理论值，误差-1.0，-1.999，-2，0.001 -0.9，-1.899，-1.8，0.099 -0.8，-1.698，-1.6，0.098 -0.7，-1.397，-1.4，0.003 -0.6，-0.996，-1，0.004 -0.5，-0.495，-0.5，0.005 -0.4，0.204，0，0.204-0.3，0.615，0.6，0.015 -0.2，1.216，1.2，0.016 -0.1，1.797，1.8，0.003 0.0，1.996，2，0.0040.1，2.193，2.2，0.007 0.2，2.792，2.8，0.008 0.3，3.293，3.3，0.007 0.4，3.594，3.6，0.006 0.5，3.896，3.9，0.004 0.6，4.437，4.4，0.037 0.7，4.998，5，0.0020.9，6.795，6.8，0.0051.0，7.993，8，0.007从实验结果可以看出，随着x的增大，计算的导数值与理论值之间的误差也在增大，但整体上相对较小。
中心差分法计算程序编程姓名：张泽伟 学号： 电话：一、中心差分法程序原理说明 1.1 中心差分法思路中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
1.2 中心差分法原理中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为：t u u u i i ∆-=-+•211 (a) 2112t u u u u i i i ∆+-=-+•• (b)而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……)由体系运动方程为：)()()(=++•••t ku t u c t u m i (c)将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t时刻的运动方程：2211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m（d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c t m u t m k u t c t m (e)由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
1.3 初始条件转化假设给定的初始条件为),0(),0(00••==u u u u （g ） 由式（g ）确定1-u 。
在零时刻速度和加速度的中心差分公式为：t u u u ∆-=-•2110 （h ） `210102t u u u u ∆+-=-•• （i ）将式（i ）消去1u 得：020012•••-∆+∆-=u t u t u u （j ）而零时刻的加速度值0••u 可以用t ＝0时的运动方程 0000=++•••ku u c u m确定即 )(1000ku u c m u --=••• （k ）这样就可以根据初始条件00,•u u 和初始荷载P ，就可以根据上式确定1-u 的值。
二阶常微分方程的中心差分求解学校:中国石油大学（华东）理学院 姓名：张道德一、 实验目的1、 构造二阶常微分边值问题:22,(),(),d uLu qu f a x bdx u a u b αβ⎧=-+=<<⎪⎨⎪==⎩其中,q f 为[,]a b 上的连续函数，0;,q αβ≥为给定常数的中心差分格式；2、 根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。
二、 中心差分格式的构造将区间[,]a b 分成N 等分，分点为： 0,1,2,,i x a ih i N =+=()/h b a N =-。
于是我们得到区间的一个网络剖分。
称为网格的节点称为步长。
得中心差分格式为：11202,1,2,,1,,.i i i h i i i i N u u u L u q u f i N h u u αβ+--+⎧=-+==-⎪⎨⎪==⎩其中式中(),()i i i i q q x f f x ==。
三、 差分格式求解根据中心差分格式可以构造出：1112222222233322212211210012101201001200N N N u f q h h u f q h h h u f q h hh q u f h h ---⎡⎤⎡⎤⎡⎤+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。
四、 举例求解我们选取的二阶常微分方程边值问题为：222242,01(0)1,(1),x d u Lu x u e x dx u u e ⎧=-+=-<<⎪⎨⎪==⎩其精确解为：2x u e=。
则我们具体求解出的解如下：1、 当10N =时，数值解与精确解为： (1) 表1、10N =时，数值解与精确解统计表x 的值 0.10.20.30.40.5 u 的数值解 1.011069 1.042744 1.096904 1.176896 1.28789 u 的精确解 1.01005 1.040811 1.094174 1.173511 1.284025 两者之差 0.001019 0.001934 0.002729 0.003385 0.003864 x 的值 0.60.70.80.9u 的数值解 1.437443 1.636363 1.900001 2.250209 u 的精确解 1.433329 1.632316 1.896481 2.247908 两者之差0.0041140.0040460.003520.002301将两者绘于同一图中如下：(2)结论：可以看出数值解与精确解之间的误差很小， 在 210-这样一个数量级上。
中心差分法的基本理论与程序设计1程序设计的目的与意义该程序通过用C语言（部分C++语言）编写了有限元中用于求解动力学问题的中心差分法，巩固和掌握了中心差分法的基本概念，提高了实际动手能力，并通过实际编程实现了中心差分法在求解某些动力学问题中的运用，加深了对该方法的理解和掌握。
2程序功能及特点该程序采用C语言（部分C++语言）实现了用于求解动力学问题的中心差分法，可以求解得到运动方程的解答，包括位移，速度和加速度。
计算简便且在算法稳定的条件下，精度较高。
3中心差分法的基本理论在动力学问题中，系统的有限元求解方程（运动方程）如下所示：阵和结点载荷向量，并分别由各自的单元矩阵和向量集成。
与静力学分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。
常微分方程的求解方法可以分为两类，即直接积分法和振型叠加法。
中心差分法属于直接积分法，其对运动方程不进行方程形式的变换而直接进行逐步数值积分。
通常的直接积分是基于两个概念，中心差分法的基本思路是用有限差分代替位移对时间的求导，将运动方程中的速度和加速度用位移的某种组合表示，然后将常微分方程组的求解问题转换为代数方程组的求解问题，并假设在每个小的时间间隔内满足运动方程，则可以求得每个时间间隔的递推公式，进而求得整个时程的反应。
在中心差分法中，加速度和速度可以用位移表示，即：(2t =∆(2t=∆而得到。
为此将加速度和速度的表达式代入上式中，即可得到中心差分法的递推公式：个离散时间点的解的递推公式，这种数值积分方法又称为逐步积分法。
需要指出和速度的表达式可知：2a +中心差分法避免了矩阵求逆的运算，是显式算法，且其为条件稳定算法，利型载荷引起的波传播问题的求解，而对于结构动力学问题则不太合适。
4 中心差分法的有限元计算格式利用中心差分法逐步求解运动方程的算法步骤如下所示： 1. 初始计算(1)(2)(3)ω=；(4)(5)(6)2.(1)(2)(3)5 程序设计5.1 程序流程图1 程序流程图各子程序主要功能为：ArrayLU ：LU Inverse ：求矩阵的转置矩阵； ArrayMVector ：矩阵和向量的乘法； LUSolve5.2 输入数据及变量说明5.2.1 输入数据该程序的原始输入数据应包括三个部分： (1)(2)(3) 确定时间步长，其中为了保证该算法的稳定性，需要满足ω=5.2.2 变量说明该程序的各个变量含义如下： (1) num ，timeStep ，dtnum ——矩阵维度； timeStep ——时间步数； dt ——时间步长；(2) M ，C ，K ，X ，V ，A ，P ，MM ，PT ，c0，c1，c2，c3M ——质量矩阵； C ——阻尼矩阵； K ——刚度矩阵； X ——位移矩阵； V ——速度矩阵；A ——加速度矩阵； P ——载荷向量； MM ——有效质量矩阵； PT —— c0,c1,c2,c3——积分常数；6算例6.1问题描述应用本程序计算一个三自由度系统，它的运动方程是：⎣当时，利用公式，可以计算得到6.2理论计算6.2.1中心差分法（理论解）(1)则起步条件为⎢⎥⎣⎦对于每一时间步长，需先计算有效载荷：由上式得到的每一时间步长的位移结果如表1所示：表10.00 0.00 0.00 0.03 0.13 0.36 0.79 1.46 2.373.42 0.00 0.03 0.19 0.58 1.26 2.24 3.434.695.846.770.401.482.974.525.826.717.227.517.858.45(2)按照相同的步骤，所得结果如下：再计算下去，位移将继续增大，这是不稳定的典型表现。
(1)中心差分法求解反应谱程序：clcload('GM.txt'); %输入原始数据n=length(GM);t=0.01;Tn=0.02:0.02:4;m=length(Tn);MaxAcc=zeros(1,m);MaxVel=zeros(1,m);MaxDis=zeros(1,m);Damp=0.05; %确定阻尼比mm=1;for T=0.02:0.02:4 %中心差分法求解结构响应时程Dis=zeros(1,n+1);Vel=zeros(1,n);Acc=zeros(1,n);w=2*pi/T;c=2*w*Damp;k=w*w;Acc(1,1)=-1*GM(1,2);Vel(1,1)=0;Dis(1,2)=0;Dis(1,1)=(t^2)/2*Acc(1,1);for i=2:1:nKe=1/(t^2)+c/(2*t);P=-1*GM(i-1,2)-(1/(t*t)-c/(2*t))*Dis(1,i-1)-(k-2/(t*t))*Dis(1,i);Dis(1,i+1)=P/Ke;endfor ii=2:1:n-1Vel(ii)=(Dis(ii+2)-Dis(ii))/(2*t);Acc(ii)=(Dis(ii+2)-2*Dis(ii+1)+Dis(ii))/(t*t)+GM(ii);endMaxAcc(1,mm)=max(abs(Acc)); %求解确定周期下结构最大响应MaxVel(1,mm)=max(abs(Vel));MaxDis(1,mm)=max(abs(Dis));mm=mm+1;endfigure(1) %绘制反应谱plot(Tn,MaxAcc(1,:))title('加速度反应谱')xlabel('自振周期(s)')ylabel('最大绝对加速度Sa(g)')gridfigure(2)plot(Tn,MaxVel(1,:))title('速度反应谱')xlabel('自振周期(s)')ylabel('最大相对速度(g*m/s)') gridfigure(3)plot(Tn,MaxDis(1,:))title('位移反应谱)xlabel('自振周期(s)')ylabel('最大相对位移(g*m*m)') grid。
第1篇一、实验目的本次实验旨在通过数值分析的方法，对工程实际问题进行建模、求解和分析。
通过学习数值方法的基本原理和算法，提高解决实际工程问题的能力。
二、实验内容1. 线性方程组的求解2. 矩阵特征值与特征向量的计算3. 函数插值与曲线拟合4. 数值微分与积分三、实验步骤1. 线性方程组的求解（1）编写程序实现高斯消元法、克劳斯消元法和列主元素法（2）设计输入界面，用户输入增广矩阵的行和列，填写系数及常数项（3）分别运用三种方法求解线性方程组，比较求解结果的正确性、数值稳定性和计算效率2. 矩阵特征值与特征向量的计算（1）编写程序实现幂法、QR算法和逆幂法（2）设计输入界面，用户输入矩阵的行和列，填写矩阵元素（3）分别运用三种方法计算矩阵的特征值与特征向量，比较求解结果的准确性和计算效率3. 函数插值与曲线拟合（1）编写程序实现拉格朗日插值、牛顿插值和样条插值（2）设计输入界面，用户输入函数的自变量和函数值，选择插值方法（3）分别运用三种方法进行函数插值，比较插值结果的准确性和光滑性4. 数值微分与积分（1）编写程序实现有限差分法、龙格-库塔法和辛普森法（2）设计输入界面，用户输入函数的导数或积分的上下限，选择数值方法（3）分别运用三种方法进行数值微分和积分，比较求解结果的准确性和计算效率四、实验结果与分析1. 线性方程组的求解通过实验，我们发现列主元素法在求解线性方程组时具有较好的数值稳定性，计算效率也较高。
而高斯消元法和克劳斯消元法在处理大型稀疏矩阵时存在一定的困难。
2. 矩阵特征值与特征向量的计算实验结果表明，QR算法和逆幂法在计算矩阵特征值与特征向量时具有较高的准确性和计算效率。
幂法在处理大型稀疏矩阵时表现出较好的性能。
3. 函数插值与曲线拟合在函数插值和曲线拟合实验中，样条插值方法具有较好的准确性和光滑性。
拉格朗日插值和牛顿插值方法在处理简单函数时表现良好，但在处理复杂函数时可能存在精度问题。
中心差分格式1、考虑问题考虑二阶常微分方程边值问题：f qu dxu d Lu =+-=22 （1） βα==)(,)(b u a u其中q ，f 为[a,b]上的连续函数，βα,为常数。
2、网格剖分与差分格式将区间[a,b]分成N 等分,分点为N i ih a x i ,,1,0,⋅⋅⋅=+=,h=(b-a)/N,于是我们得到区间I=[a,b]的网格剖分，i x 为网格节点，h 为步长。
差分格式为：.,,1,,2,120211βα==-⋅⋅⋅==++--=-+N i i i i i i i h u u N i f u q hu u u u L3、截断误差将方程（1）在节点离散化，由泰勒公式展开得)()(12)()()(2)(344222211h dx x u d h dx x u d h x u x u x u ii i i i O +⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡=+--+ 所以截断误差为)()(12)(3442h dx x u d h u R ii O +⎥⎦⎤⎢⎣⎡-= 4、数值例子x x q e x u x sin 1)()(+== 其中[]1,0∈x5、求解 由f qu dxu d Lu =+-=22，且已知 x x q e x u xsin 1)()(+== 可得x e x f x sin )(=将向量式的差分格式用矩阵形式表示出来，得到矩阵形式为⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+--+--+-2122212112112h q h q h q N⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-121N u u u =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡++-βα122212N f h f h f h 系数矩阵A=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡+--+--+-2122212112112h q h q h q N ，我们求出矩阵A 极其逆便可求得u （x ）的数值解。
6、参考文献《偏微分方程数值解法》李荣华 高等教育出版社《科学计算中的有限差分法》《MATLAB 程序设计教程》刘卫国 中国水利水电出版社附件1、程序流程图开始计算各初始条件构造系数矩阵d计算矩阵d的逆D构造格式右端矩阵B计算数值解u=D*B输出u的数值解结束2 、程序代码%zah.m%u(x)=exp(x), q(x)=1+sinx, f(x)=exp(x)*sinx 0<x<1a=input('input a:');%输入区间条件ab=input('input b:');%输入区间条件bN=input('input N:');%输入等分数Nh=(b-a)/N;%步长for i=1:N-1x(i)=a+i*h;%网格节点q(i)=1+sin(x(i));%连续函数q网格节点上各值f(i)=exp(x(i))*sin(x(i));%连续函数f网格节点上各值endc1=linspace(2+q(1)*h*h,2+q(N-1)*h*h,N-1); %系数矩阵主对角线上各值c2=linspace(-1,-1,N-2); %系数矩阵次对角线上各值d1=diag(c1,0);d2=diag(c2,-1);d3=diag(c2,1);d=d1+d2+d3 %构造系数矩阵D=inv(d); %系数矩阵的逆b1=linspace(h*h*f(1),h*h*f(N-1),N-1);for i=1:N-1if i==1b2(i)=exp(a);elseif i==N-1b2(i)=exp(b);elseb2(i)=0;endendendB1=b1+b2;%构造差分格式等式右边的向量B=B1' %差分格式等式右边向量的转置u=D*B %输出u（x）数值解输出结果：>> clear>> zahinput a:0input b:1input N:10d =2.0110 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0119 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0127 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0136 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0144 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0153 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0161 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0170 -1.0000 0 0 0 0 0 0 0 -1.0000 2.0178B =1.00110.00340.00560.00790.01020.01250.01470.01702.7375u =1.11171.23451.36851.51431.67271.84512.03312.23932.4664。
中心差分法计算单自由度体系的动力反应北京工业大学结构工程组员:胡建华 S201204111马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。
时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。
由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。
在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。
1中心差分法原理中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为： tu u u i i ∆-=-+•211 (a)2112tu u u u i i i ∆+-=-+•• (b) 而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c)将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程：02211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m（d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c tm u t m k u t c t m (e) 由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
中心差分法计算单自由度体系的动力反应北京工业大学结构工程组员:胡建华 S201204111马 恒 S201204112 陈相家S201204083 张力嘉S201204022 0前言时域逐步积分法是数值分析方法，它只假设结构本构关系在一个微小的时间步距内是线性的。
时域逐步积分法是结构动力问题中一个得到广泛研究的课题，并在结构动力反应计算中得到广泛应用。
由于引进的假设条件不同，可以有各种不同的方法，比如中心差分法，线性加速度法，平均常加速度法，Wilson-θ法等，其中中心差分法精度最高。
在本文中，通过编制中心差分法计算单自由度体系的动力反应的程序来了解其应用及稳定性。
1中心差分法原理中心差分法的基本思路：是将运动方程中的速度向量和加速度向量用位移的某种组合来表示，将微分方程组的求解问题转化为代数方程组的求解问题，并在时间区间内求得每个微小时间区间的递推公式，进而求得整个时程的反应。
中心差分法只在相隔t ∆一些离散的时间区间内满足运动方程，其基于有限差分代替位移对时间的求导（即速度和加速度），如果采用等时间步长，t t i ∆=∆，则速度与加速度的中心差分近似为： tu u u i i ∆-=-+•211 (a)2112tu u u u i i i ∆+-=-+•• (b) 而离散时间点的运动为)(),(),(i i i i i i t u u t u u t u u ••••••=== ( =i 0,1,2,3,……) 由体系运动方程为：0)()()(=++•••t ku t u c t u m i (c)将速度和加速度的差分近似公式（a ）和式（b ）代入式（c ）可以得到i t 时刻的运动方程：02211211=+∆-+∆+--+-+i i i i i i ku t u u c t u u u m（d ）在（d ）式中，假设i u 和1-i u 是已知的，即在i t 及i t 以前时刻的运动已知，则可以把已知项移到方程的右边，整理得到：12212)2()2()2(-+∆-∆-∆--=∆+∆i i i u t c tm u t m k u t c t m (e) 由式（e ）就可以根据i t 及i t 以前时刻的运动，求得1+i t 时刻的运动，如果需要可以用式（a ）和式（b ）求得体系的速度和加速度。
假设给定的初始条件为),0(),0(00••==u u u u （g ）由式（g ）确定1-u 。
在零时刻速度和加速度的中心差分公式为：tu u u ∆-=-•2110 （h ） `210102tu u u u ∆+-=-•• （i ） 将式（i ）消去1u 得：020012•••-∆+∆-=u t u t u u （j ） 而零时刻的加速度值0••u 可以用t ＝0时的运动方程 0000=++•••ku u c u m 确定即 )(1000ku u c mu --=••• （k ）这样就可以根据初始条件00,•u u 和初始荷载0P ，就可以根据上式确定1-u 的值。
2中心差分法程序编制原理：（1） 基本数据准备和初始条件计算)(1000ku u c mu --=•••020012•••-∆+∆-=u t u t u u （2） 计算等效刚度和中心差分计算公式中的相关系数t c t m k ∆+∆=2222t mk a ∆-=t ctm b ∆-∆=22（3） 根据i t 及i t 以前时刻的运动，计算1+i t 时刻的运动 1---=i i bu au P k P u i =+1 如果需要，可计算t u u u i i ∆-=-+•2112112tu u u u i i i ∆+-=-+•• （4）下一步计算用i+1代替i ，对于线弹性结构体系，重复第3步，对于非线性结构体系，重复第2步和第3步。
以上为中心差分法逐步计算公式，其具有2阶精度，即误差)(02t ∆∝ε；并且为有条件稳定，稳定条件为： πnT t ≤∆上式中，n T 为结构的自振周期，对于多自由度结构体系则为结构的最小自振周期。
3算例对于一个单层框架结构，假设楼板刚度无限大，且结构质量集中于楼层，其质量M=2000kg 、刚度K ＝50KN/m 、阻尼系数C ＝3KNs/m ，假设结构处于线弹性状态，用中心差分法计算结构的自由振动反应。
采用MATLAB 语言编程，并以单自由度体系为例进行计算,设初位移u0＝0和初速度v0=0，取不同的步长分别计算，以验证中心差分法的稳定条件πnT t ≤∆。
先计算t ∆，由稳定条件nnT t ωπ2=≤∆，而5200050000===M K n ωrad/s,则4.0522===≤∆nnT t ωπ所以本次计算取t ∆＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分别进行计算4 MATLAB 程序清单[u,v,ac]=[M,C,K,u0,v0,time,dt]% 本程序采用中心差分法计算结构的动力响应% 本程序是既可以计算单自由度体系又可以计算多自由度体系，且均假设结构体系处于线弹性状态；% ---------%%%%%输入参数%%%%%%%------------% M------------质量矩阵% C------------阻尼矩阵% K------------刚度矩阵% u0-----------初始位移% v0-----------初始速度% time---------模拟时间% dt-----------时间步长% ---%%%%%%输出值%%%%%%%%------% u--------------位移% v--------------速度% ac-------------加速度% -------%%%%%%%%中心差分法主要公式及原理%%%%%%%%%%----------- % MX''+CX'+KX=0% M*(X(t+dt)-2*X(t)+X(t-dt))/(dt^2)+C*(X(t+dt)-X(t-dt))/(2*dt)+K*X(t)=0% (M/dt^2+C/2*dt)*(X(t+dt))=-(K-2*M/dt^2)*X(t+dt)-(M/dt^2-C/2*dt)*X(t-dt)%----------------- 等效刚度Ke等效荷载Pe和相关系数a,b------------------------- % Ke=M/dt^2+C/2*dt% a=K-2*M/dt^2% b=M/dt^2-C/2*dt% Pe=-a*X(t)-b*X(t-dt)% X(t+dt)=Pe/Ke% X(t)'=(X(t+dt)-X(t-dt))/(2*dt)% X(t)''=(X(t+dt)-2*X(t)+X(t-dt))/(dt^2)% ------------------初始条件---------------------% X0''=(-C*X0'-K*X0)% X(-1)=X0-X0'*dt+X0''*(dt^2)/2% --------@Copyright by zhouhuaping(S201004232)-----clear allM=input('输入质量矩阵M ：');C=input('输入阻尼矩阵C：');K=input('输入刚度矩阵K：');u0=input('输入初始位移u0：');v0=input('输入初始速度v0：');time=input('输入模拟时间time：');dt=input('输入时间步长dt ：');[m,m]=size(K);n=time/dt; %计算步数u=zeros(m,floor(n)+1); %设定存储位移矩阵v=zeros(m,floor(n)+1); %设定存储速度矩阵ac=zeros(m,floor(n)+1); %设定存储加速度矩阵P=zeros(m,floor(n)+1); %设定存储荷载矩阵u(:,2)=u0; %给定初位移v(:,2)=v0; %给定初速度Ke=M/(dt^2)+((C)/(2*dt)); %等效刚度Ke 及系数a 、b a=K-2*M/dt^2; b=M/dt^2-C/(2*dt); for i=3:1:floor(n)+1; t=(i-2)*dt;ac(:,2)=M\(-K*u(:,2)-C*v(:,2)); %计算初加速度u(:,1)=u(:,2)-v(:,2)*dt+(ac(:,2)*(dt^2))/2; %计算(0-dt)时刻位移 Pe= -a*u(:,i-1)-b*u(:,i-2); %计算等效荷载Peu(:,i)=Ke\Pe; %计算位移 v(:,i)=(u(:,i)-u(:,i-2))/(2*dt); %计算速度 ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2))/(dt^2); %计算加速度 end%--------%%%%%%%%%%绘制位移、速度、加速度时程曲线%%%%%%%%%%%-------- t=0:dt:time;subplot(2,2,1),plot(t,u(m,:),'k-'),grid,xlabel('时间(s)'),ylabel('位移(m)'),title('顶层位移的时程曲线');subplot(2,2,2),plot(t,v(m,:),'r-'),grid,xlabel('时间(s)'),ylabel('速度(m/s)'),title('顶层速度的时程曲线');subplot(2,2,3),plot(t,ac(m,:),'b-'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),title('顶层加速度的时程曲线'); %----------end运行centraldifferent.M 文件 输入参数：K=50000; M=2000; C=3000; u0=0; v0=0;time ＝20s ；dt ＝？5中心差分法计算结果稳定性分析由以上时程图可以得到当t ∆＝0.1, 0.3, 0.4时逐步计算结果给出的结构运动趋向收敛的，即计算结果是稳定的；当t ∆＝0.41,0.42, 0.45时逐步计算结果给出的结构运动趋向发散的，即结果是不稳定的，且随着步长t ∆的增加，计算结果发散得越来越快。
