曲线拟合的研究

曲线拟合的研究

(时振宇, 刘禹, 彭波)

1 综述(历史及应用)

插值在数学发展史上是个老问题，它和拉格朗日，牛顿，高斯等著名的数学家的名字联系在一起的，它最初来源于天体计算――由若干观测值计算任意时刻星球的位置(即插值点和插值)――的需要。现在插值仍在诸如机械加工等工程技术和数据处理等科学研究中有许多直接应用.

插值常用方法有拉格朗日多项式插值，分段线性插值，三次样条插值.

拉格朗日插值是高次多项式插值(n ＋1个节点上用不超过n 次的多项式), 插值曲线光滑，误差估计有表达式，但有振荡现象，收敛性不能保证，这种插值主要用于理论分析，实际意义不大. 分段线性和三次样条插值石低次多项式插值，简单实用，收敛性有保证，但不光滑，三次样条插值的整体光滑性已大有提高，应用广泛，唯误差估计较困难.

根据一组二维数据，即平面上的若干点，确定一个一元函数，即曲线，使这些点与曲线总体来说尽量接近，这就是曲线拟合.

线性最小二乘法是解决曲线拟合的最常用方法，基本思路是, 令:

)(...)()()(2211x r a x r a x r a x f m m +++=

其中rk(x)是事先选定的一组函数，a k 是待定系数，拟合标准是使n 个点(x i ,y i ) i=1,2,…n 与y ＝f (xi)的距离的平方和最小，称最小二乘准则.

本实验所用拟合方法使高次磨削法，原理在下面叙述.

2 问题分析和算法

3.1基本思路

我们首先考虑等步长情况，不等步长可以在此基础上稍做改进而得。如图,(x 1,y 1) (x 2,y 2) (x 3,y 3) (x 4,y 4) (x 5,y 5) 为所给原始数据点中的一部分，磨光过程中应对(x 2,y 2) (x 3,y 3) (x 4,y 4) 进行切削，图示第一次切削过程。然后还需对(p 3,q 3) (p 4,q 4)第二次磨削，如此重复。高次磨光后，最后一次切削所产生的(x 2,y 2)最右侧及(x 3,y 3)最左侧的折点均向(x 2,y 2) (x 3,y 3)中点(x 23,y 23)逼近，理想情况下两点在(x 23,y 23)处重合，则磨削点(x 2,y 2)和(x 3,y 3)工作完成。切削步长大小有要求，太小则无法对x 23附近的点进行磨削，太大则在x 23附近进行了多次磨削，设每次切削后步长变为原步长的1/n ，第一次切削步长为h/a ，m 为切削次数，则有:

h =+⋯+++))n

1n 1n 1(1h a lim(m 2; 即111=+n a ; 我们取a=2，n=2可符合要求，即每次的切削步长为h/2;

图3-1 磨削算法示意图

对点 (x j , y j ) 进行一次切割将产生两个点(x j ’, y j ’), (x j+1’,y j+1’);由前面分析可知:

⎪⎪⎩

⎪⎪⎨⎧--=-=-4)('4'1j j

j j j j y y y y h x x ⎪⎪⎩

⎪⎪⎨⎧--=+=+++4)('4'111j j

j j j j y y y y h x x 如果现在由n 个点，除去两个端点不需磨削，一次磨削后，将产生2*(n-2)+2个新的点，其坐标用行向量x ’，y ’表示；在求下一次磨削时，将x ’，y ’复制给行向量x ，y ，新的切削后的行纵坐标x ’，y ’可用上述方法得到。如果我们确定了循环次数，则拟合完成.

3.2 初始点的移动

为了使磨削后的点过原始点，需将原始点在磨削前移动一定距离，如图: x 代表原始点坐标, x ’代表移动后的坐标.

'8

1'43'8142''''32131222x x x x x x x x ++=+--

=

如此可得:

⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⋅⋅⋅⋅⋅⋅⋅⋅⋅

⋅⋅⋅⋅⋅⋅⋅⋅=1125.000125.075.0125.00125.0125.00125.075.0125.0000125.01A

T k x x x x X ],...,,,[321=

T X A X ⋅=

X A X T ⋅=-1

y 和x 的计算雷同, 故不再重复

图3-2 计算y 的示意图

3.3 不等步长的推广和单侧点的处理 当步长不等时，h 不是一个常数，可以在每次使用h 时计算出两点间的步长.如:

⎪⎪⎩

⎪⎪⎨⎧--=--=--4'4'1

1j j j j j j j y y yj y x x x x 当(x1,y1)及(xk , yk)不重合时，曲线不闭合，两个端点都只有单侧信息，我们只能对其一半进行磨削处理，剩下的一部分由直线连接；而当(x 1,y 1)和(x k , y k )重合后，磨削最后形成一个闭合的曲线，两端点亦需移动，矩阵B 如下。点的切削和前基本一致，唯一不同在于端点 (x 1,y 1)，(x k , y k )应联系(x 2,y 2) (x k-1,y k-1)进行磨削.

⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅=75.0125.000125.0125.075.0125.000125.0125.000125.075.0125.0125.000125.075.0B

3误差分析

f(xi+h/2)=(I+D*h/2+D^2*h^2/8)*f(xi)+o(h^2/2)

f(xi-/2)=(I-*h/2+D^2*h^2/8)*f(xi)+o(h^2/2)

S1(x) = f(xi+h/2)-(f(xi+h/2)-f(xi-h/2))/h*(-x+xi+h/2) = f(xi+h/2)-D*f(xi)(xi+h/2-x)

f(x)=(I+D*(x-xi)+D^2*(x-xi)^2/2)*f(xi)+o((x-xi)^2)

|f(x)-S1(x)|=D^2/2*((x-xi)^2-h^2/4)

当x=xi时，|f(x)-S1(x)|=D^2*f(xi)*h^2/8;

所以，若f(xi)的二次导数越小，点的连线越平和，误差越小4软件使用说明

4.1 界面功能介绍(curve.m)

①文件操作按钮