谷根代数值分析--上机实习报告

华北电力大学数值分析实验报告实验题目 非线性方程求解专业、班级硕电力学号姓名课程编号实验类型验证性学时实验（上机）地点完成时间任课教师谷根代评分 一、实验目的及要求1 掌握非线性方程数值解法的基本概念、方法；2 理解二分法、牛顿法、割线法、斯蒂芬森法求解非线性方程的原理、方法、步骤，会用计算机进行编程求解；3 掌握运用Matlab进行数值计算的能力。二、研究、解答以下问题问题：分别用二分法、牛顿迭代法、割线法、斯蒂芬森迭代法求方程的根，观察不同初试值条件下的收敛性，并给出你的结论【解】：2.1 算法分析（1） 二分法二分法伪代码如下所示;输入：，定义函数，给定精度输出：待解方程的根。while if 由于本题中恒大于0，不满足二分法的条件，故本例无法通过二分法求出方程的根。（2）牛顿迭代法牛顿迭代法迭代法迭代公式为：a).对单根条件下，牛顿迭代法平方收敛，牛顿迭代法程序框图如下所示：b). 对重根条件下，此时迭代公式修改为此时，牛顿迭代法至少平方收敛。（3）割线法割线法的迭代公式为：割线法是超线性收敛，其程序流程图为（4）斯蒂芬森迭代法斯蒂芬森迭代法的计算格式为斯蒂芬森迭代法程序流程图与牛顿迭代法类似。2.2迭代结果分别取初试值为，运用牛顿迭代法、割线法（）、，斯蒂芬森法求得根如下(误差精度，数据结够为（format long）)牛顿迭代法割线法斯蒂芬森法未考虑重根考虑重根10.75246478873231.01478873239430.69997632410230.989561085620320.79542362451241.00003635756930.69995263720680.999984872531130.83069135397221.00000000022030.75243574328150.999999999968240.85971455036201.00000000000000.77847702762921.000000000000050.88364872892290.807419563780460.90342111644770.830533421415170.91977955805340.8514440083581640.99999581654950.9998927397768850.9999924537033表1 时，不同迭代方法迭代结果牛顿迭代法割线法斯蒂芬森法未考虑重根考虑重根11.85294117647051.11764705882351.9677760435023 0.918698062683021.72588259231141.00224878683861.94066438382440.918698062683031.61639394463971.00000084252241.94066438382440.918698062683041.52231054116801.00000000000011.94066438382440.918698062683051.44169908122671.940664383824461.37283109171661.584950464173771.31416248727131.5849504641737651.00000853976191.0003773276313681.00000494199540.9999924537033951.0000075015608表2 时，不同迭代方法迭代结果2.3 分析通过表1、表2不难分析出如下结论：（1）牛顿法、割线法、斯蒂芬森法本题下均收敛，可以求出误差范围内的根。（2）对于重根情况下，牛顿法的收敛速度会变慢变为线性收敛。为解决这一问题，可以对牛顿迭代公式进行修正，此时至少能够保证迭代是平方收敛，使得求解速度大大加快。（3）割线法的收敛速度小于牛顿法。（4）运用斯蒂芬森法可以对牛顿迭代法进行加速。三、本次实验的难点分析1、本次试验是数值分析课程的第一次试验，重温了matlab的基本语法与操作。由于本科时使用matlab主要是电力系统仿真即主要运用simlinuk模块，对于matlab中自定义函数的建立与运用还有欠缺。这是本次试验中耗时较多的部分。2、本次试验，原设计在牛顿迭代法程序中加入判断导数值是否是个小值（导数值为零）的过程，以防止分母出现较小的值。但是程序运行过程中导数值始终会接近于零而跳出，得不出解。为此不得不放弃这一思路。四、参考文献 1谷根代，杨晓忠等数值分析M北京：科学出版社，20112马东升，雷永军数值计算方法M北京：机械工业出版社，20083周路，陈渝等译数值方法（Matlab版）M北京：电子工业出版社，2010数值分析实验报告实验题目 插值与拟合专业、班级学号姓名课程编号实验类型验证性学时实验（上机）地点完成时间任课教师谷根代评分一、实验目的及要求1 掌握插值、拟合的基本概念、方法；2 理解牛顿插值法、三次样条插值原理、方法、步骤，理解最小二乘法的原理方法，会用计算机进行编程求解，懂得Runge现象及预防措施；3 掌握运用Matlab进行数值计算的能力。二、研究、解答以下问题问题：1 区间作等距划分：以为节点对函数进行插值逼近。（1）分别取，用牛顿插值对进行逼近，并在同一坐标系下作出函数的图形，进行比较，写出插值函数对的逼近程度与节点个数之间的关系，并分析原因；（2）使用三次样条插值对进行逼近，在同一坐标系下画出图形，观察样条插值函数对的逼近程度与节点个数的关系；（3）整体插值有和局限性？如何避免？2 已知一组数据如下，求其拟合曲线。01234567891023478101114161819106.42108.2109.5110109.93110.49110.59110.6110.76111111.2（1）求以上数据形如的拟合曲线及其平方误差；（2）求以上数据形如的拟合曲线及其平方误差；（3）通过观察（1）、（2）的结果，写出你对数据拟合的认识。2.1 算法分析2.1.1 插值（1）牛顿插值法Lagrange 插值公式结构紧凑，便于理论分析。利用插值基函数也容易到插值多项式的，Lagrange 插值公式的缺点是，当插值节点增加，或其位置变化时，全部插值基函数均要随之变化，从而整个插值公式的结构也发生变化，这在实际计算中是非常不利的。下面引入的Newton 插值公式可以克服这个缺点。Newton 插值多项式可以灵活地增加插值节点进行递推计算。该公式形式对称，结构紧凑，因而容易编写计算程序。Newton 插值多项式为：其中：称之为 f (x)在上的k 阶均差（或差商）。实际计算中，可以由差商表计算。，如下所示：表1 差商表牛顿迭代法程序流程图如下所示：（2）三次样条插值分段低次插值光滑性较差，对于像高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数。三次样条函数 记为 它是定义在区间a, b 上的函数， 满足以下两个条件：1). S(x) 在每一个小区间上是一个三次多项式函数 ；2). 在整个区间a,b上，其二阶导数存在且连续。 即在每个节点处的二阶导数连续。问题：给定函数在n+1个节点处的函数值为 。要求：求一个三次样条函数，使其满足：如何确定三次样条函数在每一个小区间上的三次多项式函数的系数呢？这里只简介确定系数的思想。对于三次样条插值呢，每一个区间上的三次函数的四个参数，在该区间上由两个端点的函数值只能够产生两个方程，仅此不足以唯一确定四个参数。注意到三次样条函数对整体光滑性要求，其二阶导数存在且连续，从全局的角度上考虑参数个数与方程个数的关系如下：参数：每个小段上4个，n个小段共计4n个。方程：1) 每个小段上由给定函数值得到2个，n个小段共计2n个；2) 光滑性要求每一个内部节点的一阶二阶导数连续，得出其左右导数相等，因此，每个节点产生2个方程，共计2(n-1) 个 。现在得到了4n-2个方程，还差两个。为此，常用的方法是对边界节点除函数值外附加要求，这就是所谓的边界条件。需要两个，正好左右两个端点各一个。常用如下三类边界条件：第一边界条件 给定两个边界节点的一阶导数值：m0,mn, 即： (x0)=m0, (xn)=mn。第二边界条件 给定两个边界节点的二阶导数值：M0,Mn, 即： (x0)=M0, (xn)=Mn。特别地，当 M0和Mn都为零时，称为自然边界条件。周期性边界条件 在两个边界的函数值，一阶导数值以及二阶导数值均相等：即 (x0)= (xn)； (x0)= (xn) 。2.1.2 拟合在科学实验或者统计研究中常常需要从一组测得的数据中去求得自变量和因变量之间的一个近似解析表达式使得误差最小, 这就是给定个点进行曲线拟合的问题。曲线拟合常用次多项式来进行拟合，即求解下面方程就可进行曲线拟合：对于指数形式的拟合曲线，可以令，则有，用上式依然可以解除2.2 数据处理与问题分析2.2.1 插值（1）运用所编程序，分别求出时，得出插值函数，并在同一坐标系下画出原函数与插值函数图像，分别如下所示：图1 n=1时原函数与插值函数图象图2 n=5时原函数与插值函数图象图3 n=10时原函数与插值函数图象图4 n=20时原函数与插值函数图象图5 n=25时原函数与插值函数图象(2) 运用所编程序，求出下三次样条插值函数，并在同一坐标系下作出原函数图象与插值函数图象。如下图所示：图6 n=1时原函数与插值函数图象图7 n=5时原函数与插值函数图象图8 n=10时原函数与插值函数图象2.2.2 拟合（1）形如的拟合，根据所编matlab程序，得到拟合方程为（2）形如的拟合，根据所编matlab程序，得到拟合方程为（3）将（1）、（2）与原数据散点图图象做在同一坐标系中，如下如所示：通过此图，可以发现，指数形式的拟合函数更加和原数据走势拟合，然而不管哪一种拟合形式，都不能保证所有的数据点都在拟合曲线上。所以，拟合可以表现出数据的真题趋势，而且，通过拟合可以合理的排除一些误差点。2.3 问题分析2.3.1 插值（1）由2.2.1 中时插值函数与原函数的图象分析可知，插值点越多，插值函数越逼近与原函数。然而，当把此题的区间扩大，改为时，情况却发生变化，时插值函数与原函数的图象如下所示：图9 n=3时原函数与插值函数图象图10 n=5时原函数与插值函数图象图11 n=10时原函数与插值函数图象图12 n=25时原函数与插值函数图象通过上述四图发现，随着插值节点数目的增多，区间中间的确越来越逼近原函数，然而区间两端发生了较大震荡，即发生了Runge现象。这是由插值余项的表达式所决定的。（2）为进一步说明三次样条插值下插值节点数目与插值函数的关系，同上，将区间扩大至，别求出时，得出插值函数，并在同一坐标系下画出原函数与插值函数图像，分别如下所示：图13 n=1时原函数与插值函数图象图14 n=5时原函数与插值函数图象图15 n=10时原函数与插值函数图象通过图68，图1315，可以发现：插值节点数目越多，三次样条插值函数越逼近与原函数。三次样条函数没有Runge现象，这是三次样条插值将区间分成了n份。（3）对于整体插值，如果节点数目较少，插值函数有不怎么精确，如图9，图10所示；如果节点数目增多，就有可能会发生Runge现象，如图11、12所示。解决方法是分段低次插值。2 拟合拟合部分问题分析见2.2.2三、本次实验的难点分析（1）由于自学，本次试验对牛顿插值法的差商公式理解有一点偏差。（3）本次试验中插值部分，对于牛顿插值，随着插值节点的增大，在区间内却不出现Runge现象，对于这点感到迷惑。为了说明问题，只能将插值区间选为。此时，随着插值节点数目的增多，出现了明显的Runge现象。四、参考文献 1谷根代，杨晓忠等数值分析M北京：科学出版社，20112马东升，雷永军数值计算方法M北京：机械工业出版社，20083周路，陈渝等译数值方法（Matlab版）M北京：电子工业出版社，2010数值分析实验报告实验题目 数值微分和数值积分专业、班级学号姓名课程编号实验类型验证性学时实验（上机）地点完成时间任课教师谷根代评分一、实验目的及要求1 掌握数值积分与数值微分的基本概念、方法；2 理解Newton-Cotes公式、复合求积公式、龙贝格积分公式、Gauss型积分公式等求解数值积分与数值微分的原理、方法、步骤，会用计算机进行编程求解；3 掌握运用Matlab进行数值计算的能力。二、研究、解答以下问题问题：设计区间分半求积算法、龙贝杰求积算法和自适应辛普森算法的程序，观察时，积分的结果，并给出相应评价。【解】： 1 区间半分积分法设将区间分成等分，既有个子区间，分点为，其中步长。对于子区间，再取其中点作新节点，求积分值若有，则取为积分值。流程图如下：输入命令格式如下：T2,T1,nn=Variable_step(-1,1,1e-5,n)，其中n为题中函数参数。迭代次数11010050010.540302-0.8390710.862319-0.88384920.489547-0.3486210.672401-0.38167830.480839-0.1018540.626418-0.28098840.478895-0.1304000.615014-0.25707750.478423-0.1377420.612169-0.25117560.478306-0.1395890.000711-0.00147170.478277-0.140038-0.007551-0.0080678-0.140152-0.0092420.0014079-0.140181-0.009656-0.00125410-0.009757-0.00173511-0.009772-0.00184812-0.0018762 龙贝格求积算法用复化梯形公式，取区间长度为，复化梯形公式的值与原积分值之间存在关系：该式称为欧拉-麦克劳林求和公式。上式还表明，用近似的截断误差为：龙贝格求积算法计算步骤如下：(1) 计算： (2) 对分区间并计算和：，其中为新分点的函数值之和。(3) 计算与：，(4) 利用外推公式：，()，()，直至求出。(5) 判断是否真，其中为给定的精度。若真，则，否则重复(2)、(3)、(4)的计算。对于的龙贝格积分结果，分别如下表所示：11.080605000020.5403020.36020200030.4895470.4726280.4801230040.4808390.4779370.4782910.478262050.4788950.4782470.4782680.4782670.478267=11-1.678143000002-0.839072-0.55938100003-0.348620-0.185136-0.1601870004-0.101855-0.019599-0.008564-0.006157005-0.130400-0.139915-0.147936-0.150148-0.15071306-0.137742-0.140189-0.140207-0.140084-0.140045-0.140035=1011.7246380.0000000.0000000.0000000.0000000.00000020.8623190.5748790.0000000.0000000.0000000.00000030.6724010.6090950.6113760.0000000.0000000.00000040.6264180.6110900.6112240.6112210.0000000.00000050.6150140.6112130.6112210.6112210.6112210.00000060.6121690.6112210.6112210.6112210.6112210.611221=1001-1.7676990.0000000.0000000.0000000.0000000.0000002-0.883849-0.5892330.0000000.0000000.0000000.0000003-0.381678-0.214287-0.1892910.0000000.0000000.0000004-0.280988-0.247425-0.249634-0.2505920.0000000.0000005-0.257077-0.249107-0.249219-0.249213-0.2492070.0000006-0.251175-0.249207-0.249214-0.249214-0.249214-0.249214=500 3 自适应求积方法自适应求积算法是按照被积函数在区间的变化形态来安排求积节点的。自适应Simpson积分法方法如下：输入：端点精度要求；输出：近似值,分层数，分点；置；while 1 n1=length(ep); n=length(p0); if n=1 break,end h=p0(2)-p0(1); s0=h/6*(f(p0(1),nn)+4*f(p0(1)+h/2),nn)+f(p0(1)+h,nn); s1=h/12*(f(p0(1),nn)+4*f(p0(1)+h/4),nn)+f(p0(1)+h/2),nn); s2=h/12*(f(p0(1)+h/2,nn)+4*f(p0(1)+3*h/4),nn)+f(p0(1)+h,nn); if abs(s0-s1-s2)=2 ep=ep(2:n1); end q=q+1; else m=m+1; p0=p0(1),p0(1)+h/2,p0(2:n); if n1=1 ep=ep(1)/2,ep(1)/2,ep(2:n1); end if q=0 p=p0; else p=p(1:q),p0; end end end当时，积分结果 ,产生的积分节点为12345678-1-0.5-0.2500.250.50.751下图表现了被积函数的曲线和自适应法产生的求积节点的关系当时，积分结果 ,产生的积分节点为1-1160.1252-0.9375170.18753-0.875180.254-0.8125190.3755-0.75200.43756-0.6875210.57-0.625220.1258-0.5625230.18759-0.5240.2510-0.4375250.37511-0.375260.437512-0.25270.513-0.1875280.562514-0.125290.62515030下图表现了被积函数的曲线和自适应法产生的求积节点的关系当时，积分结果 ,产生的积分节点为12345-1-0.500.51下图表现了被积函数的曲线和自适应法产生的求积节点的关系时，积分结果 ,产生的积分节点数目为985个，不在赘述。分析（1）综合（1）、（2）、（3）的结果发现在时三种积分方法基本相等，这表明在震荡较为平缓的情况下区间分半求积算法、龙贝杰求积算法和自适应辛普森算法均是适用的；（2）通过对龙贝杰求积算法与区间分半求积算法的迭代次数的比较，使用龙贝杰求积算法确实能够提高运算速度。（3）当逐渐增大，此时区间分半求积算法的求值结果最先出现偏差，故对震荡剧烈的被积函数区间分半求积算法适用性不强。（4）自适应求积算法是按照被积函数在区间的变化形态来安排求积节点的。对于震荡剧烈的被积函数该法较为适用。三、本次实验的难点分析本次试验难点主要有以下几点：（1）对于龙贝杰求积算法概念与原理的理解上有难度；（2）采用三种不同积分算法对较大的情况下的积分结果不同，理解原因有难度；（3）本题采用自适应积分法由于待积函数是偶函数且积分区间关于轴对称，理应积分节点也是关于左右对称的，但是对于的情况显然不对称，对于原因理解有难度。四、参考文献 1谷根代，杨晓忠等数值分析M北京：科学出版社，20112马东升，雷永军数值计算方法M北京：机械工业出版社，20083周路，陈渝等译数值方法（Matlab版）M北京：电子工业出版社，2010数值分析实验报告实验题目 常微分方程数值解法专业、班级学号姓名课程编号实验类型验证性学时实验（上机）地点完成时间任课教师谷根代评分 一、实验目的及要求1 掌握解常微分方程数值解法的基本概念、方法；2理解Euler法、隐式Euler法求解常微分方程的原理、方法、步骤，会用计算机进行编程求解；3 掌握运用Matlab进行数值计算的能力。二、研究、解答以下问题问题：分别取步长h=0.1, 0.05, 0.01,0.005,，用显式欧拉法和隐式欧拉法求解y=-50y, y(0)=100，由结果分析算法的稳定性。【解】： 1、显示欧拉法源程序：function y=xianshi(h)N=1000;a=;a(1)=100;for k=1:N a(k+1)=a(k)-50*a(k)*h; if (a(k)-a(k+1) xianshi(0.1)a = 100 -400 1600 xianshi(0.05)a = 100 -150 225 xianshi(0.01)a = Columns 1 through 13 100.0000 50.0000 25.0000 12.5000 6.2500 3.1250 1.5625 0.7813 0.3906 0.1953 0.0977 0.0488 0.0244 Columns 14 through 21 0.0122 0.0061 0.0031 0.0015 0.0008 0.0004 0.0002 0.0001 xianshi(0.005)a = Columns 1 through 13 100.0000 75.0000 56.2500 42.1875 31.6406 23.7305 17.7979 13.3484 10.0113 7.5085 5.6314 4.2235 3.1676 Columns 14 through 26 2.3757 1.7818 1.3363 1.0023 0.7517 0.5638 0.4228 0.3171 0.2378 0.1784 0.1338 0.1003 0.0753 Columns 27 through 39 0.0564 0.0423 0.0317 0.0238 0.0179 0.0134 0.0100 0.0075 0.0057 0.0042 0.0032 0.0024 0.0018 Columns 40 through 460.0013 0.0010 0.0008 0.0006 0.0004 0.0003 0.00022、隐式欧拉法源程序：function y=yinshi(h)N=1000;a=;a(1)=100;for k=1:N a(k+1)=a(k)/(1+50*h); if (a(k)-a(k+1) yinshi(0.1)a = 100.0000 16.6667 2.7778 0.4630 0.0772 0.0129 0.0021 0.0004 0.0001 0.0000 yinshi(0.05)a = 100.0000 28.5714 8.1633 2.3324 0.6664 0.1904 0.0544 0.0155 0.0044 0.0013 0.0004 0.0001 0.0000 yinshi(0.01)a = Columns 1 through 13 100.0000 66.6667 44.4444 29.6296 19.7531 13.1687 8.7791 5.8528 3.9018 2.6012 1.7342 1.1561 0.7707 Columns 14 through 26 0.5138 0.3425 0.2284 0.1522 0.1015 0.0677 0.0451 0.0301 0.0200 0.0134 0.0089 0.0059 0.0040 Columns 27 through 34 0.0026 0.0018 0.0012 0.0008 0.0005 0.0003 0.0002 0.0002 yinshi(0.005)a = Columns 1 through 13 100.0000 80.0000 64.0000 51.2000 40.9600 32.7680 26.2144 20.9715 16.7772 13.4218 10.7374 8.5899 6.8719 Columns 14 through 26 5.4976 4.3980 3.5184 2.8147 2.2518 1.8014 1.4412 1.1529 0.9223 0.7379 0.5903 0.4722 0.3778 Columns 27 through 39 0.3022 0.2418 0.1934 0.1547 0.1238 0.0990 0.0792 0.0634 0.0507 0.0406 0.0325 0.0260 0.0208 Columns 40 through 52 0.0166 0.0133 0.0106 0.0085 0.0068 0.0054 0.0044 0.0035 0.0028 0.0022 0.0018 0.0014 0.0011 Columns 53 through 57 0.0009 0.0007 0.0006 0.0005 0.0004三、本次实验的难点分析（1）（2）（3）四、参考文献1谷根代，杨晓忠等数值分析M北京：科学出版社，20112马东升，雷永军数值计算方法M北京：机械工业出版社，20083周路，陈渝等译数值方法（Matlab版）M北京：电子工业出版社，