曲线拟合的数值方法 |
您所在的位置：网站首页 › 指数方程的解法步骤 › 曲线拟合的数值方法 |
曲线拟合的数值方法
|
《数值计算方法》系列总目录 第一章 误差序列实验 第二章 非线性方程f(x)=0求根的数值方法 第三章 CAD模型旋转和AX=B的数值方法 第四章 插值与多项式逼近的数值计算方法 第五章 曲线拟合的数值方法 第六章 数值微分计算方法 第七章 数值积分计算方法 第八章 数值优化方法 第五章 一、算法原理1、最小二乘拟合曲线2、曲线拟合原理3、样条插值4、贝塞尔曲线5、三角多项式原理 二、实验内容及核心算法代码1、最小二乘拟合曲线原理实现2、曲线拟合原理实现3、样条函数插值原理实现4、三角多项式原理实现5、贝塞尔曲线原理实现 三、结果分析1、最小二乘拟合曲线实例分析2、曲线拟合实例分析3、样条函数插值实例分析4、三角多项式实例分析5、贝塞尔曲线实例分析 四、结论 一、算法原理 1、最小二乘拟合曲线 在科学和工程试验中,经常产生一组数据(
x
1
,
y
1
x_1,y_1
x1,y1),…,(
x
N
,
y
N
x_N,y_N
xN,yN),这里的横坐标{
x
k
x_k
xk}是明确的的。数值方法的目标之一是确定一个将这些变量联系起来的函数
y
=
f
(
x
)
y = f(x)
y=f(x)，即通过算法形成一种数学刻画的公式来寻求其中的规律。通常会选择一类可用的函数并确定它们的系数，然而这类函数不是随便选取的，必须根据实际的实验数据进行判断，选取最符合的函数模型，因此选择函数的可能性是多种多样的。一般会根据物理情况采用一个基本数学模型来确定函数的形式。这一节主要分析线性函数,其形式为
y
=
f
(
x
)
=
A
x
+
B
y = f(x) = Ax + B
y=f(x)=Ax+B 根据插值多项式逼近原理，如果所有的数值{
x
k
x_k
xk},{
y
k
y_k
yk}已知有多位有效数字精度,则能成功地使用多项式插值;否则,不能使用多项式插值。因为插值多项式是根据已有的精确的实验数据进行多项式拟合的，其拟合结果一定会经过插值点，如果插值点数据本身存在较大误差，则逼近结果就会出现更大误差。然而,许多试验数据可能并没有如此高的精度,而且通常在试验中还存在试验误差,所以尽管记录的{
x
k
x_k
xk},{
y
k
y_k
yk}有多位有效数字,但还是与其真实值存在范围外的误差，即真实值
f
(
x
k
)
f(x_k)
f(xk)满足
f
(
x
k
)
=
y
k
+
e
k
f(x_k)=y_k+e_k
f(xk)=yk+ek，其中，
e
k
e_k
ek表示测量误差，因此需要寻找一种经过测试点附近(不总是穿过)的最佳线性逼近表达式。下面首先对误差(又称偏差或残差)概念进行一个介绍：
e
k
=
f
(
x
k
)
−
y
k
,
f
o
r
1
≤
k
≤
N
{e_k} = f({x_k}) - {y_k},{\rm{ }}for{\rm{ }}1 \le k \le N
ek=f(xk)−yk,for1≤k≤N 其中，
y
k
y_k
yk测量实验值，
f
(
x
k
)
f(x_k)
f(xk)表示真实值，误差还有其他表示如下形式：
E
∞
(
f
)
=
max
1
≤
k
≤
N
{
∣
f
(
x
k
)
−
y
k
∣
}
{E_\infty }(f) = \mathop {\max }\limits_{1 \le k \le N} \{ |f({x_k}) - {y_k}|\}
E∞(f)=1≤k≤Nmax{∣f(xk)−yk∣}
E
1
(
f
)
=
1
N
∑
k
=
1
N
∣
f
(
x
k
)
−
y
k
∣
{E_1}(f) = \frac{1}{N}\sum\limits_{k = 1}^N {\left| {f({x_k}) - {y_k}} \right|}
E1(f)=N1k=1∑N∣f(xk)−yk∣
E
2
(
f
)
=
(
1
N
∑
k
=
1
N
∣
f
(
x
k
)
−
y
k
∣
2
)
1
/
2
{E_2}(f) = {\left( {\frac{1}{N}\sum\limits_{k = 1}^N {{{\left| {f({x_k}) - {y_k}} \right|}^2}} } \right)^{1/2}}
E2(f)=(N1k=1∑N∣f(xk)−yk∣2)1/2 那么如何确定系数
A
，
B
A，B
A，B使得曲线拟合达到我们想要的效果呢？其中，最小二乘法是在进行曲线拟合时常用的一种求解拟合曲线表达式的方法，无论是拟合何种曲线，其系数求解基本思想都是基于最小二乘法的。分析可知，设 是一个N个点的集合，其中横坐标是确定的。那么，理想的最小二乘拟合曲线 的拟合结果应满足是均方根误差最小的曲线。那么，所选定的
A
，
B
A，B
A，B构成的拟合曲线应当使
E
2
(
f
)
E_2(f)
E2(f)的值当且仅当所选取的
A
，
B
A，B
A，B值处最小。换言之，也就是说数据点到曲线的垂直距离平方和最小，即误差的平方和，其最小值为
N
(
E
2
(
f
)
2
)
=
∑
k
=
1
N
∣
A
x
k
+
B
−
y
k
∣
2
N({E_2}{(f)^2}) = \sum\limits_{k = 1}^N {{{\left| {A{x_k} + B - {y_k}} \right|}^2}}
N(E2(f)2)=k=1∑N∣Axk+B−yk∣2。那么，基于实验结果的所求出选取
A
，
B
A，B
A，B的值的表达式到底是怎样的呢？因此，有如下推导过程： 设
{
(
x
k
,
y
k
)
}
k
=
1
N
\{ ({x_k},{y_k})\} _{k = 1}^N
{(xk,yk)}k=1N有N个点，其中横坐标是确定的，当数据点到曲线的垂直距离平方和最小时，
A
，
B
A，B
A，B取值点应在误差表达式的极值处，即表达式对于
A
，
B
A，B
A，B两变量的偏导数为0，即
∂
E
(
A
,
B
)
∂
A
=
∑
k
=
1
N
2
(
A
x
k
+
B
−
y
k
)
x
k
=
0
∂
E
(
A
,
B
)
∂
B
=
∑
k
=
1
N
2
(
A
x
k
+
B
−
y
k
)
=
0
\frac{{\partial E(A,B)}}{{\partial A}} = \sum\limits_{k = 1}^N {2\left( {A{x_k} + B - {y_k}} \right){x_k}} = 0\\ \frac{{\partial E(A,B)}}{{\partial B}} = \sum\limits_{k = 1}^N {2\left( {A{x_k} + B - {y_k}} \right)} = 0
∂A∂E(A,B)=k=1∑N2(Axk+B−yk)xk=0∂B∂E(A,B)=k=1∑N2(Axk+B−yk)=0 通过移项变换即得到关于系数
A
，
B
A，B
A，B的表达式如下
(
∑
k
=
1
N
x
k
2
)
A
+
(
∑
k
=
1
N
x
k
)
B
=
∑
k
=
1
N
x
k
y
k
(
∑
k
=
1
N
x
k
)
A
+
N
B
=
∑
k
=
1
N
y
k
\left( {\sum\limits_{k = 1}^N {x_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)B = \sum\limits_{k = 1}^N {{x_k}{y_k}} \\ \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)A + NB = \sum\limits_{k = 1}^N {{y_k}}
(k=1∑Nxk2)A+(k=1∑Nxk)B=k=1∑Nxkyk(k=1∑Nxk)A+NB=k=1∑Nyk 并称式（5）表示的方程组为正规方程。得到的曲线称为最小二乘拟合曲线。 除此之外，对于形如
f
(
x
)
=
A
x
M
f(x) = A{x^M}
f(x)=AxM幂函数拟合，其中M为已知常数，根据线性曲线的拟合推导过程，写出幂函数数据点到曲线的垂直距离平方和表达式，即
E
(
A
)
=
∑
k
=
1
N
(
A
x
k
M
−
y
k
)
2
E(A) = \sum\limits_{k = 1}^N {{{\left( {Ax_k^M - {y_k}} \right)}^2}}
E(A)=k=1∑N(AxkM−yk)2 求出A为其极小值点时候的值，即表达式对于A的偏导数为0，即
E
′
(
A
)
=
2
∑
k
=
1
N
(
A
x
k
M
−
y
k
)
x
k
M
=
2
∑
k
=
1
N
(
A
x
k
2
M
−
x
k
M
y
k
)
=
0
E'(A) = 2\sum\limits_{k = 1}^N {\left( {Ax_k^M - {y_k}} \right)} x_k^M = 2\sum\limits_{k = 1}^N {\left( {Ax_k^{2M} - x_k^M{y_k}} \right) = 0}
E′(A)=2k=1∑N(AxkM−yk)xkM=2k=1∑N(Axk2M−xkMyk)=0 最终可以得出幂函数拟合曲线的系数A的表达式为 对于不是线性关系的曲线拟合，可以对其进行线性化操作再利用最小二乘法求出其系数的表达式即可。也就是说，对于一组不具有线性关系的数据点集合，可以对其进行诸如取对数、指数的方法将其变换为具有线性关系的表达式，即将数据点
x
k
,
y
k
{xk,yk}
xk,yk变换成
X
Y
XY
XY平面上具有线性关系的点集
X
k
,
Y
k
{Xk,Yk}
Xk,Yk,再利用最小二乘法求出变换后的曲线表达式，最后再变换为原变量的表达式。这个过程称为数据线性化。例如，对于非线性函数[y = C{e^{Ax}}]，可以对其进行取对数变化，即
ln
(
y
)
=
A
x
+
ln
(
C
)
\ln (y) = Ax + \ln (C)
ln(y)=Ax+ln(C) 并引入变量， ，得到变换后的线性关系式 ，根据最小二乘法的结果，系数
A
，
B
A，B
A，B表达式为
(
∑
k
=
1
N
X
k
2
)
A
+
(
∑
k
=
1
N
X
k
)
B
=
∑
k
=
1
N
X
k
Y
k
(
∑
k
=
1
N
X
k
)
A
+
N
B
=
∑
k
=
1
N
Y
k
\left( {\sum\limits_{k = 1}^N {X_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {X_k^{}} } \right)B = \sum\limits_{k = 1}^N {{X_k}{Y_k}} \\ \left( {\sum\limits_{k = 1}^N {X_k^{}} } \right)A + NB = \sum\limits_{k = 1}^N {{Y_k}}
(k=1∑NXk2)A+(k=1∑NXk)B=k=1∑NXkYk(k=1∑NXk)A+NB=k=1∑NYk
C
=
e
B
C = {e^B}
C=eB 则可以得到拟合曲线的表达式。然而，也可以直接写出数据点到曲线的垂直距离平方和表达式，并求出其极值点，即为拟合曲线的系数表达式，其推导过程如下：
E
(
A
,
C
)
=
∑
k
=
1
N
(
C
e
A
x
k
−
y
k
)
2
l
e
t
{
∂
E
∂
A
=
2
∑
k
=
1
N
(
C
e
A
x
k
−
y
k
)
(
C
x
k
e
A
x
k
)
=
0
∂
E
∂
C
=
2
∑
k
=
1
N
(
C
e
A
x
k
−
y
k
)
(
e
A
x
k
)
s
o
{
C
∑
k
=
1
N
x
k
e
2
A
x
k
−
∑
k
=
1
N
y
k
x
k
e
A
x
k
=
0
C
∑
k
=
1
N
e
A
x
k
−
∑
k
=
1
N
y
k
e
A
x
k
=
0
E(A,C) = \sum\limits_{k = 1}^N {{{\left( {C{e^{A{x_k}}} - {y_k}} \right)}^2}} \\ {\rm{let}}\left\{ \begin{array}{l} \frac{{\partial E}}{{\partial A}} = 2\sum\limits_{k = 1}^N {\left( {C{e^{A{x_k}}} - {y_k}} \right)} \left( {C{x_k}{e^{A{x_k}}}} \right){\rm{ = 0}}\\ \frac{{\partial E}}{{\partial C}} = 2\sum\limits_{k = 1}^N {\left( {C{e^{A{x_k}}} - {y_k}} \right)} \left( {{e^{A{x_k}}}} \right) \end{array} \right.\\ {\rm{so}}\\ \left\{ \begin{array}{l} C\sum\limits_{k = 1}^N {{x_k}{e^{2A{x_k}}} - } \sum\limits_{k = 1}^N {{y_k}{x_k}{e^{A{x_k}}}} = 0\\ C\sum\limits_{k = 1}^N {{e^{A{x_k}}} - } \sum\limits_{k = 1}^N {{y_k}{e^{A{x_k}}}} = 0 \end{array} \right.
E(A,C)=k=1∑N(CeAxk−yk)2let⎩⎪⎪⎨⎪⎪⎧∂A∂E=2k=1∑N(CeAxk−yk)(CxkeAxk)=0∂C∂E=2k=1∑N(CeAxk−yk)(eAxk)so⎩⎪⎪⎨⎪⎪⎧Ck=1∑Nxke2Axk−k=1∑NykxkeAxk=0Ck=1∑NeAxk−k=1∑NykeAxk=0 其中未知数
A
，
C
A，C
A，C是非线性的，可用牛顿法求解，可以看出，直接利用最小二乘法求解更加复杂。除此之外，在求解其他函数的拟合曲线的时候，也可以进行类似的数据线性变换。具体如下图所示。 对N+1个点 的多项式插值经常不令人满意。一个N阶多项式可能有N-1个相对极大值和极小值,同时曲线可能会摆动,以经过这些点。另一个方法是将图形分段,每段为一个低阶多项式
S
(
x
)
S(x)
S(x),并在相邻点
(
x
k
,
y
k
)
(x_k,y_k)
(xk,yk)和
(
x
k
+
1
,
y
k
+
1
)
(x_{k+1},y_{k+1})
(xk+1,yk+1)之间进行插值.这样可以保证插值多项式没有高阶多项式插值的摆动性，而且这种插值方法一般是N段的分段函数。这样的分段拟合的函数称之为样条函数。其中一个重要的性质就是插值多项式经过结点。那么，接下来要思考的问题就是每一小段的插值多项式应该选择什么样的多项式。 首先考虑一阶多项式，即一阶拉格朗日多项式。然而不难想到这种拟合方式在大多数情况下下是不适用的，因为拟合后的曲线为一段一段的折线段，显然不符合拟合的要求。然后再来考虑二阶多项式，二阶多项式能满足平滑的特点，但是会发现其在偶数结点处其曲率变化很大，这很有可能导致有非期望的弯曲或畸变。而且二次样条的导数在其结点处不连续。那么如果使用三次样条，则可以使其一阶二阶导数都连续，然而这个条件仅仅需要其多项式系数满足一定要求即可。不再继续考虑更高阶的样条函数是因为三次样条已经能满足需求，不必要牺牲产生震荡失真的代价去考虑高次样条，而且高次样条的系数增多，使求解变得极其困难。那么三次样条的多项式系数该如何确定呢？ 首先，我们需要给出三次样条需要满足的条件使其是一个光滑连续的函数，即一阶二阶导数在结点处连续。  Mark：紧压样条在端点有斜率。紧压样条可想像为用外力使柔软而有弹性的木杆经过数据点,并在端点处使其具有固定斜率。这样的样条对于画经过多个点的光滑曲线的绘图员相当有用。natural样条。存在惟一的三次样条曲线,其一阶导数的边界条件是 =0和 =0，于是求解的方程组为  Mark：natural样条是柔软有弹性的木杆经过所有数据点形成的曲线,但让端点的斜率自由地在某一位置保持平衡,使得曲线的摇摆最小。它在对有多位有效数字精度的试验数据进行曲线拟合时很有用。外推样条。存在惟一的三次样条曲线,其中通过对点
x
1
x_1
x1和
x
2
x_2
x2进行外推得到 ,同时通过对点
x
N
−
1
x_{N-1}
xN−1和
x
N
−
2
x_{N-2}
xN−2进行外推得到 。于是求解的方程组为  Mark：外推样条等价于端点处的三次多项式曲线是相邻三次多项式曲线的扩展形成的样条，也就是说,样条曲线在区间
[
x
0
,
x
2
]
[x_0,x_2]
[x0,x2]内形成单个三次多项式曲线,同时在区间
[
x
N
−
2
,
x
N
]
[x_{N-2},x_N]
[xN−2,xN]内形成另一个三次多项式曲线。抛物线终结样条。存在唯一的三次样条，其中在区间
[
x
0
,
x
1
]
[x_0,x_1]
[x0,x1]内 ,而在
[
x
N
−
1
,
x
N
]
[x_{N-1},x_N]
[xN−1,xN]内， .于是求解的方程组为  Mark：在区间
[
x
0
,
x
1
]
[x_0,x_1]
[x0,x1]内
S
"
(
x
)
=
0
S"(x)=0
S"(x)=0使得在区间
[
x
0
,
x
1
]
[x_0,x_1]
[x0,x1]内三次多项式曲线退化为二次多项式曲线,同时在区间
[
x
N
−
1
,
x
N
]
[x_{N-1},x_N]
[xN−1,xN]内也发生同样的情况。端点曲率调整样条。存在唯一的三次样条曲线，其中二阶导数的边界条件
S
’
’
(
a
)
S’’(a)
S’’(a)和
S
’
’
(
b
)
S’’(b)
S’’(b)使确定的。于是求解的方程组为  Mark：直接对
S
’
’
(
a
)
S’’(a)
S’’(a)和
S
’
’
(
b
)
S’’(b)
S’’(b)赋值，可调整每个端点的曲率。
综上所述，无论是哪种策略，其都给出了在端点处的一阶或二阶导数取值以及相应的线性方程组，以此来得出除端点外的所有结点二阶导数值，即 m k m_k mk。然后再通过五种策略的根据其他点二阶导数值来确定端点的二阶导数值，从而确定所有系数，最终求出三次样条插值函数。 4、贝塞尔曲线20世纪70年代,雷诺汽车公司的Pierre Bezier和雪铁龙汽车公司的Paul de Castaliau各自独立推导出了CAD/CAM应用中的贝塞尔(Bezier)曲线。这些参数多项式是一类逼近样条。贝塞尔曲线是计算机图形学(CAD/CAM,计算机辅助几何设计)中曲线和曲面表示的基本方法。最初推导的贝塞尔曲线是用递归方法隐式定义的。将它们用伯恩斯坦多项式(Bemnstein polynomial)显式地定义,有助于推导贝塞尔曲线的性质。 其中，N阶伯恩斯坦多项式的定义为 B i , N ( t ) = ( N i ) t i ( 1 − t ) N − i i = 0 , 1 , 2 , . . . , N ( N i ) = N ! t ! ( N − i ) ! {B_{i,N}}(t) = \left( \begin{array}{l} N\\ i \end{array} \right){t^i}{(1 - t)^{N - i}}{\rm{ }}i = 0,1,2,...,N{\rm{ }}\left( \begin{array}{l} N\\ i \end{array} \right) = \frac{{N!}}{{t!(N - i)!}} Bi,N(t)=(Ni)ti(1−t)N−ii=0,1,2,...,N(Ni)=t!(N−i)!N! 一般，N阶伯恩斯坦多项式由N+1项。其中，伯恩斯坦多项式具有几个重要的性质。首先，伯恩斯坦多项式具有递归关系，表现为 B 0 , 0 ( t ) = 1 , B i , N ( t ) = 0 , i < 0 o r i > N B i , N ( t ) = ( 1 − t ) B i , N − 1 ( t ) + t B i − 1 , N − 1 ( t ) = 0 , i = 1 , 2 , 3 , . . . , N − 1 {B_{0,0}}(t) = 1,{B_{i,N}}(t) = 0,{\rm{ }}i < 0{\rm{ }}or{\rm{ }}i > N\\ {B_{i,N}}(t) = \left( {1 - t} \right){B_{i,N - 1}}(t) + t{B_{i - 1,N - 1}}(t) = 0,{\rm{ }}i = 1,2,3,...,N - 1 B0,0(t)=1,Bi,N(t)=0,iNBi,N(t)=(1−t)Bi,N−1(t)+tBi−1,N−1(t)=0,i=1,2,3,...,N−1 因此，任意一个伯恩斯坦多项式可由其递归关系生成，并为后来的贝塞尔函数求解提供手段。其次，伯恩斯坦多项式还具有非负性，即在区间 [ 0 , 1 ] [0,1] [0,1]上是非负函数。因为伯恩斯坦多项式可以看作n次重复实验的伯努利分布的概率公式，其中t可以视为概率p，即只能在区间 [ 0 , 1 ] [0,1] [0,1]内取值，因此，其值也就是概率一定是处在 [ 0 , 1 ] [0,1] [0,1]内的非负数。第三个性质，伯恩斯坦多项式还具有规范性，即 ∑ i = 0 N B i , N ( t ) = 1 \sum\limits_{i = 0}^N {{B_{i,N}}(t)} = 1 i=0∑NBi,N(t)=1 其规范性可以由二项式定理证明。第四个性质，也就是伯恩斯坦多项式的导数性质，即 d d t B i , N ( t ) = N ( B i − 1 , N − 1 ( t ) − B i , N − 1 ( t ) ) \frac{d}{{dt}}{B_{i,N}}(t) = N\left( {{B_{i - 1,N - 1}}(t) - {B_{i,N - 1}}(t)} \right) dtdBi,N(t)=N(Bi−1,N−1(t)−Bi,N−1(t)) 可以对伯恩斯坦多项式的定义是直接求导，即可证明其导数性质。最后一个性质，也是其可以推导出贝塞尔曲线的最重要的一个性质，就是其可以作为基，即N阶伯恩斯坦多项式组成阶数小于等于N的所有多项式的一个基空间。性质5说明,任何阶数小于等于N的多项式都可以惟一地表示为伯恩斯坦多项式的线性组合。 给出贝塞尔曲线的定义，给定一个控制点集 ，一条N阶贝塞尔曲线定义为N阶伯恩斯坦多项式的加权和。即N阶贝塞尔曲线可以表示成 P ( t ) = ∑ i = 0 N P i B i , N ( t ) , i = 0 , 1 , . . . , N , t ∈ [ 0 , 1 ] P(t) = \sum\limits_{i = 0}^N {{P_i}{B_{i,N}}(t)} ,{\rm{ }}i = 0,1,...,N,t \in [0,1] P(t)=i=0∑NPiBi,N(t),i=0,1,...,N,t∈[0,1] 其中 为N阶伯恩斯坦多项式。控制点是表示平面中x和y坐标的有序对。可将控制点作为向量,而对应的伯恩斯坦多项式作为标量处理,这样公式(28)可参数化表示为 P ( t ) = ( x ( t ) , y ( t ) ) P(t)=(x(t),y(t)) P(t)=(x(t),y(t)).其中 x ( t ) = ∑ i = 0 N x i B i , N ( t ) , y ( t ) = ∑ i = 0 N y i B i , N ( t ) , x(t) = \sum\limits_{i = 0}^N {{x_i}{B_{i,N}}(t)} ,{\rm{ }}y(t) = \sum\limits_{i = 0}^N {{y_i}{B_{i,N}}(t)} ,{\rm{ }} x(t)=i=0∑NxiBi,N(t),y(t)=i=0∑NyiBi,N(t), 其中 0 < t < 1 0{P_i}{B_{i,N}}(0)} = {P_0},P(1) = \sum\limits_{i = 0}^N {{P_i}{B_{i,N}}(1)} = {P_N} |
除此之外，还可以使用线性最小二乘法来进行拟合。即设有N个数据点
{
(
x
k
,
y
k
)
}
\{(x_k,y_k)\}
{(xk,yk)}，并给定了M个线性独立的函数
f
i
(
x
)
{f_i(x)}
fi(x)。利用这些线性独立的函数的线性组合表示函数
f
(
x
)
f(x)
f(x)，即
f
(
x
)
=
∑
j
=
1
M
c
j
f
j
(
x
)
f(x) = \sum\limits_{j = 1}^M {{c_j}{f_j}(x)}
f(x)=j=1∑Mcjfj(x) 同样，利用最小二乘法求解系数的思想，使其数据点到曲线的垂直距离平方和最小，即其表达式对于每一个参数
c
i
c_i
ci的偏导数为0，则可以得到其正规方程：
∑
j
=
1
M
(
∑
k
=
1
N
f
i
(
x
k
)
f
j
(
x
k
)
)
c
j
=
∑
k
=
1
N
f
i
(
x
k
)
y
k
,
i
=
1
,
2
,
.
.
.
,
M
\sum\limits_{j = 1}^M {\left( {\sum\limits_{k = 1}^N {{f_i}({x_k}){f_j}({x_k})} } \right)} {c_j} = \sum\limits_{k = 1}^N {{f_i}({x_k})} {y_k},i = 1,2,...,M
j=1∑M(k=1∑Nfi(xk)fj(xk))cj=k=1∑Nfi(xk)yk,i=1,2,...,M 不难看出，等式左边可以表达为矩阵乘积的结果，因此构造如下矩阵： 
则等式右边可以写成如下形式 
而方程左边未知参数c的系数可以写为矩阵F与其逆矩阵的乘积，因此，方程组可以表达为如下线性方程组的形式：
F
′
F
C
=
F
′
Y
F'FC = F'Y
F′FC=F′Y 对于不太大的M，可以利用第三章求解线性方程组的方法求解系数c，从而得出函数的表达式。因此，当所取线性独立函数
f
i
(
x
)
{f_i(x)}
fi(x)为幂函数时，则函数
f
(
x
)
f(x)
f(x)可以表示为M阶多项式的形式，然而事实上，这种多项式拟合通常具有多项式摆动的特性，因此会造成拟合结果的失真，而且阶数越高越容易发生，因此通常情况下只会取到二阶多项式进行拟合，除非已知使用的多项式是真实的。因此，进行最小二乘抛物线拟合的函数表达式为
y
=
f
(
x
)
=
C
+
B
x
+
A
x
2
y = f(x) = C + Bx + A{x^2}
y=f(x)=C+Bx+Ax2 同样，利用最小二乘法求解系数的思想，使其数据点到曲线的垂直距离平方和最小，即其表达式对于每一个参数的偏导数为0，则可以得到其正规方程：
(
∑
k
=
1
N
x
k
4
)
A
+
(
∑
k
=
1
N
x
k
3
)
B
+
(
∑
k
=
1
N
x
k
2
)
C
=
∑
k
=
1
N
y
k
x
k
2
(
∑
k
=
1
N
x
k
3
)
A
+
(
∑
k
=
1
N
x
k
2
)
B
+
(
∑
k
=
1
N
x
k
)
C
=
∑
k
=
1
N
y
k
x
k
(
∑
k
=
1
N
x
k
2
)
A
+
(
∑
k
=
1
N
x
k
)
B
+
N
C
=
∑
k
=
1
N
y
k
\left( {\sum\limits_{k = 1}^N {x_k^4} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^3} } \right)B + \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)C = \sum\limits_{k = 1}^N {{y_k}x_k^2} \\ \left( {\sum\limits_{k = 1}^N {x_k^3} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)B + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)C = \sum\limits_{k = 1}^N {{y_k}x_k^{}} \\ \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)B + NC = \sum\limits_{k = 1}^N {{y_k}}
(k=1∑Nxk4)A+(k=1∑Nxk3)B+(k=1∑Nxk2)C=k=1∑Nykxk2(k=1∑Nxk3)A+(k=1∑Nxk2)B+(k=1∑Nxk)C=k=1∑Nykxk(k=1∑Nxk2)A+(k=1∑Nxk)B+NC=k=1∑Nyk 以上系数表达式即为进行二阶多项式拟合的各阶项系数表达式。
其中，性质Ⅰ描述了由分段三次多项式构成的
S
(
x
)
S(x)
S(x)。性质
Ⅱ
Ⅱ
Ⅱ描述了对给定数据点集的分段三次插值。性质
Ⅲ
Ⅲ
Ⅲ和性质
Ⅳ
Ⅳ
Ⅳ保证了分段三次多项式函数是一个光滑连续函数。性质
V
V
V保证了函数的二阶导数也是连续的。 是否可能构造一个三次样条满足性质
Ⅰ
Ⅰ
Ⅰ到性质
V
V
V?每个三次多项式
S
(
x
)
S(x)
S(x)有4个未知常数，因此需要求解
4
N
4N
4N个系数。宽松的讲,要确定
4
N
4N
4N个自由度或条件。数据点提供了
N
+
1
N+1
N+1个条件,性质
Ⅲ
Ⅲ
Ⅲ、性质
Ⅳ
Ⅳ
Ⅳ和性质
V
V
V都提供了
N
−
1
N-1
N−1个条件,总共确定了
N
+
1
+
3
(
N
−
1
)
=
4
N
−
2
N+1+3(N - 1)=4N-2
N+1+3(N−1)=4N−2个条件。这样剩下了两个自由度。.可称之为端点约束( end-point constraints)涉及在点
x
0
x0
x0和
x
N
xN
xN处的导数
S
′
(
x
)
S'(x)
S′(x)或
S
"
(
x
)
S"(x)
S"(x)。通过推导，由于二阶导数满足拉格朗日多项式，因此令其满足性质
V
V
V，再通过积分两次，可以求解出
S
(
x
)
S(x)
S(x)的表达式，再利用性质
I
I
II
II求出积分后不定系数的表达式，因此其表达式为
S
k
(
x
)
=
−
m
k
6
h
k
(
x
k
+
1
−
x
)
3
+
m
k
+
1
6
h
k
(
x
−
x
k
)
3
+
(
y
k
h
k
−
m
k
h
k
6
)
(
x
k
+
1
−
x
)
+
(
y
k
+
1
h
k
−
m
k
+
1
h
k
6
)
(
x
−
x
k
)
{S_k}(x) = - \frac{{{m_k}}}{{6{h_k}}}{\left( {{x_{k + 1}} - x} \right)^3} + \frac{{{m_{k + 1}}}}{{6{h_k}}}{\left( {x - {x_k}} \right)^3} + \left( {\frac{{{y_k}}}{{{h_k}}} - \frac{{{m_k}{h_k}}}{6}} \right)\left( {{x_{k + 1}} - x} \right) + \left( {\frac{{{y_{k + 1}}}}{{{h_k}}} - \frac{{{m_{k + 1}}{h_k}}}{6}} \right)\left( {x - {x_k}} \right)
Sk(x)=−6hkmk(xk+1−x)3+6hkmk+1(x−xk)3+(hkyk−6mkhk)(xk+1−x)+(hkyk+1−6mk+1hk)(x−xk) 其中系数
h
k
h_k
hk为第
k
k
k个和第
k
+
1
k+1
k+1个结点的差分值，
m
k
=
S
′
′
(
x
k
)
{m_k} = S''({x_k})
mk=S′′(xk)。那么只剩下一阶导数的性质没有用上。通过推到，利用一阶导数的性质，可以确定
m
k
m_k
mk的递推关系式，即
h
k
−
1
m
k
−
1
+
2
(
h
k
−
1
+
h
k
)
m
k
+
h
k
m
k
+
1
=
u
k
,
u
k
=
6
(
d
k
−
d
k
−
1
)
,
k
=
1
,
2
,
.
.
.
,
N
−
1
{h_{k - 1}}{m_{k - 1}} + 2\left( {{h_{k - 1}} + {h_k}} \right){m_k} + {h_k}{m_{k + 1}} = {u_k},{u_k} = 6({d_k} - {d_{k - 1}}),k = 1,2,...,N - 1
hk−1mk−1+2(hk−1+hk)mk+hkmk+1=uk,uk=6(dk−dk−1),k=1,2,...,N−1 可以发现，可以构造一个有包含
N
+
1
N+1
N+1个未知数具有
N
−
1
N-1
N−1个方程的线性方程组，即端点处的二阶导数值，因此有以下五种策略求解线性方程组。 
因此采用策略可以定出
m
0
m_0
m0的值，则第一个方程可以化简为
2
(
h
0
+
h
1
)
m
1
+
h
1
m
2
=
u
1
−
h
0
m
0
,
w
h
e
n
k
=
1
2\left( {{h_0} + {h_1}} \right){m_1} + {h_1}{m_2} = {u_1} - {h_0}{m_0},{\rm{ }}when{\rm{ }}k = 1
2(h0+h1)m1+h1m2=u1−h0m0,whenk=1 同理，已知
m
N
m_N
mN的值也可以知道最后一个方程
h
N
−
2
m
N
−
2
+
2
(
h
N
−
2
+
h
N
−
1
)
m
N
−
1
=
u
N
−
1
−
h
N
−
1
m
N
,
w
h
e
n
k
=
N
−
1
{h_{N - 2}}{m_{N - 2}} + 2\left( {{h_{N - 2}} + {h_{N - 1}}} \right){m_{N - 1}} = {u_{N - 1}} - {h_{N - 1}}{m_N},{\rm{ }}when{\rm{ }}k = N - 1
hN−2mN−2+2(hN−2+hN−1)mN−1=uN−1−hN−1mN,whenk=N−1 这样，剩下的系数就可以用一个带状的线性方程组来表示，即 
表示为
H
M
=
V
HM=V
HM=V，利用第三章的几种求解线性方程组的方法，即可确定所有的结点二阶导数值，将这些值带回原插值函数，得到最终的插值函数多项式。到此，三次样条函数插值已经全部构造完成。 对于五种不同的策略，其求解的 有所不同。具体为：【本文地址】