插值与曲线拟合实验报告

1、数值计算方法插值与拟合实验报告摘要：通过实验介绍插值方法中常见的拉格朗日插值，线性分段插值和牛顿前插公式，分析计算各种方法的插值余项。在曲线拟合方面使用两种不同类型的曲线来拟合同一组数据，并计算残差向量范数，比较不同曲线拟合的效果，在此例上给出优劣的判断。关键词：拉格朗日插值；线性分段插值；牛顿前插公式；曲线拟合引言在工程和科学计算中经常碰到只知道离散的数据测量点而需要匹配其变量之间的数学函数表达式的情况，这就需要插值和拟合的数值方法来解决这些问题。插值法是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，也是离散函数逼近的重要方法，利用它可通过函数在有限点处的取值状况

2、，估算出函数在其他点处的近似值。曲线拟合则是用连续曲线近似地刻画或比拟平面上离散点所表示的坐标之间的函数关系，在各个方面也有着愈加广泛的应用。1 算法介绍1.1 拉格朗日插值法1.1.1 算法理论对某个多项式函数，已知有给定的k+1个取值点：其中对应着自变量的位置，而对应着函数在这个位置的取值。假设任意两个不同的xj都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为：拉格朗日基本多项式 的特点是在 上取值为1，在其它的点 上取值为0。对于给定的 个点：，拉格朗日插值法的思路是找到一个在一点取值为1，而在其他点取值都是0的

3、多项式。这样，多项式在点取值为，而在其他点取值都是0。而多项式就可以满足在其它点取值为0的多项式容易找到，例如：它在点取值为：。由于已经假定两两互不相同，因此上面的取值不等于0。于是，将多项式除以这个取值，就得到一个满足“在取值为1，而在其他点取值都是0的多项式”：这就是拉格朗日基本多项式。拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，非常繁琐。此外，当插值点比较多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却

4、会和“实际上”的值之间有很大的偏差。这类现象也被称为龙格现象，解决的办法是分段用较低次数的插值多项式。1.2.2 算法描述对于已给定的点 和待估计的点的横坐标x，如上述理论，将其值代入计算出插值基函数的值，然后根据公式：计算出纵坐标的估计值，由此完成对该点的插值过程，其中k为该点插值的阶数。1.2 线性分段插值1.2.1 算法理论首先介绍线型插值，假设我们已知坐标 与 要得到 区间内某一位置x在直线上的值。根据图中所示，我们得到由于x值已知，所以可以从公式得到 y 的值已知求的过程与以上过程相同，只是与要进行交换。线性插值经常用于已知函数在两点的值要近似获得其它点数值的方法，这种近似方法的误差

5、定义为：其中p表示上面定义的线性插值多项式根据罗尔定理，我们可以证明：如果f有二阶连续导数，那么误差范围是 正如所看到的，函数上两点之间的近似随着所近似的函数的二阶导数的增大而逐渐变差。从直观上来看也是这样：函数的曲率越大，简单线性插值近似的误差也越大。而分段线型插值就是把需要插值的区间依已给定的点 分为k-1段，每段利用其端点进行线型插值。1.2.2 算法描述利用已给定的点 对插值区间分为段，将每段的端点与 作为数据点利用公式在所构成的区间进行线性插值。1.3 牛顿插值法1.3.1 算法理论注意在下面的讲述中着重介绍差商的概念，设：问题是如何根据插值条件来计算待定系数 由 知，。由 因而，其

6、中， 称为函数 在点的一阶差商。由 知：因而其中， 称为函数在 点的二阶差商。实际上，它是一阶差商的差商。一般地，如果已知一阶差商，那么就可以计算二阶差商类似于上述过程不断地推导下去，可得其中， 称为函数f (x)在相应点处的n 阶差商。这些结果说明，要保证上面构造多项式的方法具有期望的优点，就必须要有一个好的计算系数 的方法。按上述方式构造插值多项式的方法叫做牛顿插值法。根据插值多项式的惟一性知，其截断误差与拉格朗日插值法相同，即：但也可以表示成差商形式。这是因为以 为节点的多项式从而 于是 的截断误差可表为顺便指出，因为牛顿插值多项式具有性质：所以，类似于逐次线性插值法，也可以把上述和式中

7、的第二项看成是估计 的一种实用误差估计式。与差商概念密切联系的另一个概念是差分，它是指在等距节点上函数值的差。所谓等距节点，是指对给定的常数h（称为步长），节点称为 处的一阶向前差分；称为 处的一阶向后差分。一阶差分的差分称为二阶差分，即 称为 处的二阶向前差分。一般地，m 阶向前和向后差分可定义如下：差商与差分之间存在关系：其中 是介于 之间的某个数。将此式代入， 即可得到牛顿前插公式。1.3.2 算法描述利用已给定的点 和待估计的点的横坐标，代入牛顿前插公式，先计算差分，然后算出该点纵坐标的估计值。1.4 曲线拟合1.4.1 算法理论令待求的未知量为，它们可由个直接测量通过下列函数关系求得

8、：若为真值，由上述已知函数求出真值，若其测量值为，则对应的误差为 最小二乘法可定量表示为： 将拟合函数取线性函数是一种简单的数据拟合方法，利用数据点确定线性拟合函数 称为对数据的线性拟合。对于线性拟合问题，需要求函数 的最小值点，该问题的几何背景是寻求一条直线，使该直线与数据表所确定的平面散点的纵向距离的平方和最小。由函数对两个变量求导得： 其余等于零，得正规方程组 也可将其矩阵形式写出来即：解得的值，将其代入 即可得到拟合线性函数。为了确定数据拟合问题，选用幂函数作为函数类，则 这就是多项式拟合函数.为了确定拟合函数的系数，需要求解正规方程组 也可以用矩阵形式表示为 解得最优解即可，将其代入

9、即可得到拟合多项式。若用指数函数对数据进行拟合，则使用线性化方法，对等式两边求对数得：令得 对和进行线性拟合求出B， C,由此可进一步求出A。1.2.2 算法描述抛物线拟合数据,令然后对和进行线性拟合，求出参数a， b；指数曲线拟合数据，先对纵坐标数据取对数，然后如上述理论对lny 和 x 进行线性拟合，进一步求出参数a，b。3结果分析3.1 拉格朗日插值法和分段线性插值对下表中给定离的散点序列进行五次插值和分段线性插值xi0.300.420.500.580.660.72yi1.044031.084621.118031.156031.198171.23223求 时的函数近似值，得到五次插值结果

10、为1., 1., 1,分段线性插值结果为1., 1.14178, 。现在分析两种方法的插值余项：拉格朗日插值法的截断误差为在本例中 对 均只能确定是0.30到0.72之间的某个值，而n为5；对线性分段插值，其误差项与拉格两日插值类似，其n的取值为1， 但 的取值范围比较小，均在插值点最近的两个数据点之间，对，在0.30, 0.42区间之内，为0.50, 0.58, 为0.66, 0.72。由于随着插值次数的增高，误差项之中所除的阶乘部分快速增大而使误差相对缩小，但是次数过高时由于高阶导数上下界无法确定，而且多项式次数升高也可能会产生龙格现象使误差增大。插值曲线如下图，其中第一幅图为五次插值曲线

11、，另一幅为分段线性插值曲线。由图中可以看出，相比之下，分段插值曲线较为曲折，有些点处不可导。3.2 利用牛顿前插公式进行插值对下表中给定离的散点序列利用牛顿前插公式进行插值xi0.1250.2500.3750.5000.6250.750f(xi)0.7960.7730.7440.7040.6560.602先得到差分表：X00.0.0.0.0.0.X1-0.-0.-0.-0.-0.X2-0.-0.-0.-0.X3-0.0.0.X40.-0.X5-0.再分别求 x = 0.158 和 x = 0.636 的三次插值的值，求得的结果为0., 0.。牛顿插值与拉格朗日插值法插值余项相同，即：这里计算的

12、是三次插值，所以n取3，但由于当 x = 0.636 时离前四个数据点较远，会使 增大，所以估计x = 0.636 是插值余项的值会比较大。3.3 曲线拟合试分别用抛物线 和指数曲线拟合下列数据：xi12.53.54yi3.81.5026.033.0画出数据点和两条拟合曲线，其中实线为抛物线拟合结果，虚线为指数曲线拟合：求得抛物线拟合的残差向量的2范数为10.524， 指数曲线拟合为15.012， 可见抛物线对这组数据拟合效果较好。4 结论拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，

13、非常繁琐。这时可以用牛顿插值法来代替，而牛顿插值法的插值余项与拉格朗日插值法相同。此外，当插值点比较多的时候，高次插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差。这时采用的办法则是采用分段插值以降低多项式次数。参考文献： 1John H.Mathews, Kurtis D.Fink. 数值方法M. 北京：电子工业出版社，2010.2孙志忠, 吴宏伟, 袁慰平, 闻震初. 计算方法与实习M. 第5版. 南京:东南大学出版社, 2011.附录（python语言代码）五次插值和分段插值程序：from num

14、py import *from matplotlib.pyplot import *Xi = array(0.30, 0.42, 0.50, 0.58, 0.66, 0.72)Yi = array(1.04403, 1.08462, 1.11803,1.15603, 1.19817, 1.23223)X = array(0.32, 0.55, 0.68)def lagrint(X, Xi, Yi):n = len(Xi)L = array(0.0 * n)Y = array(0.0 * len(X)for i in range(len(X):dxi = Xi - Xifor j in rang

15、e(0, n):tmp = Xij - Xinum = prod(dxik for k in range(n) if k != j)den = prod(tmpk for k in range(n) if k != j)Lj = num / denYi = sum(Yi * L)return Ydef seglinearint(X, Xi, Yi):n = len(X)Y = array(0.0 * n)for i in range(n):j = 0;while Xi Xij: j = j + 1x0, y0 = Xij-1, Yij-1x1, y1 = Xij, YijYi = y0 + (

16、Xi - x0) * (y1 - y0) / (x1 - x0)return Yif _name_ = _main_:print(LAGRANGE INTERPOLATION:)print(lagrint(X, Xi, Yi)print(PIECEWISE INTERPOLATION:)print(seglinearint(X, Xi, Yi)X = linspace(0.30, 0.72, 1000)Y1 = lagrint(X, Xi, Yi)Y2 = seglinearint(X, Xi, Yi)plot(X, Y1, color = red)scatter(Xi, Yi)title(F

17、igure of Lagrange Interpolation)figure()plot(X, Y2, color = blue)scatter(Xi, Yi)title(Figure of Piecewise Method)show()牛顿插值公式程序：from numpy import *from sys import *Xi = array(0.125, 0.250, 0.375, 0.500, 0.625, 0.750)Yi = array(0.796, 0.773, 0.744, 0.704, 0.656, 0.602)X = array(0.158, 0.636)def diffe

18、rence(index, degree, Yi):if degree = 0:return Yiindexelse:diff1 = difference(index+1, degree-1, Yi)diff2 = difference(index, degree-1, Yi)return (diff1 - diff2)def printDifTable(Yi):n = len(Yi)print(DIFFERENCE TABLE:)for i in range(n):stdout.write(X%d % i)for j in range(n - i):diffe = difference(j,

19、i, Yi)stdout.write(%10f % diffe)stdout.write(n)stdout.write(n)def newtonint(X, Xi, Yi, degree):h = Xi1 - Xi0Y = array(0.0 * len(X)dif = array(0.0 * (degree + 1)L = array(1.0 * (degree + 1)for i in range(degree):difi = difference(0, i, Yi)for i in range(len(X):t = (Xi - Xi0) / htfact = 1.0for j in ra

20、nge(degree):Lj = tfact / prod(range(1, j+1)tfact = tfact * (t - j)Yi = sum(L * dif)return Yif _name_ = _main_:printDifTable(Yi)print(INTERPOLATION RESULT:)print(newtonint(X, Xi, Yi, degree = 3)曲线拟合程序：from numpy import *from matplotlib.pyplot import *Xi = array(1.0, 2.5, 3.5, 4.0)Yi = array(3.8, 1.5, 26.0, 33.0)norm = lambda X: sqrt(sum(X * 2)def fitAplusBx2(Xi, Yi):Xi = Xi * 2xmean = mean(Xi)ymean = mean(Yi)sumx2 = sum(Xi - xmean) * 2)