曲线拟合的应用

曲线拟合的应用

摘要：在实际问题中，常常会从一组数据中筛选出对自己有用的部分，这样的问题可转化为寻找一种函数曲线去拟合这些数据，在解决这类问题的数据处理和误差分析中应用最广泛的是曲线拟合。它不但可以提高数据处理效率，而且还能保证相当的精确度。 关键词：曲线拟合，最小二乘法，应用

1.直线拟合

直线拟合数据点(,)(1,2,

)i i x y i n =的最小二乘法，即找一个一次函数y Ax B =+，使二元函数

21(,)()n

i i i E A B Ax B y ==+-∑

达到最小。

由多元函数取得极值的必要条件知，由方程组：

1

1(,)2()0(,)2()10

n

i i i i n

i i

i E A B x y x A

E A B x y B

==??=+-?=??????=+-?=???∑∑ 化简可得正规方程组：

2

111

1

1()()()n n n i i i i i i i n n

i i i i A x B x x y A x nB y

=====?+=????+=??∑∑∑∑∑ （1-1）

由方程组（1-1）解出,A B ，即得一次函数y Ax B =+为所求的拟合直线.

2.幂函数拟合

在某些情况下的拟合函数M

y Ax =，其中M 是一个已知常数 设}

{

1

(,)n

i i i x y =有n 个点，最小二乘幂函数拟合曲线M

y Ax =，求函数()E A 的最小值？

21

()()n

M

i

i i E A Ax

y ==

-∑

对上式求关于A 的导数： 1

()2()()n

M M i

i i i E A Ax

y x ='=-?∑

令导数等于0，化简得： 21

1

(

)()0n

n M

M i

i i i i A x

x y ==-=∑∑

121

()

n

M

i

i i n M

i

i x

y A x

===

∑∑

即：M

y Ax =为所求的拟合曲线。

3.指数拟合

3.1 求解Ax y Ce =的非线性最小二乘法

设给定一组点集(,)(1,2,

)i i x y i n =，需要拟合指数曲线

采用非线性最小二乘法求下式的最小值： 21

(,)()i

n

Ax i i E A C Ce

y ==

-∑

（3.1-1） 对上式分别求关于的偏导数，并令导数等于0

1

1(,)2()()0(,)2()()0

i i

i i n Ax Ax i i i n

Ax Ax i

i E A C Ce y Ce x A

E A C Ce y e C ==??=-?=??????=-?=???∑∑ （3.1-2） 化简可得正规方程组：

211

211()()0()()0i i

i i n n Ax Ax i i i i i n n

Ax Ax i i i C x e x y e C e y e ====?-=????-=??∑∑∑∑ （3.1-3） 方程（3.1-3）对于未知数A 和C 是线性的，可用牛顿法求解。但这是一个耗时的计算，而且迭代

需要好的A 和C 的初始值。

3.2 求解Ax y Ce =线性化方法

设给定一组点集(,)(1,2,

)i i x y i n =，求指数函数Ax y Ce =的拟合曲线.

对上式两边同取对数得： ln ln y C Ax =+

令： ln Y y =，ln B C =，X x = 得： Y AX B =+

xy 平面上的初始点集(,)i i x y 变成平面上XY 的点集(,)(,ln )i i i i X Y x y =，这个过程称为数据线性化。这样可用最小二乘拟合曲线拟合点集}

{

1

(,)n

i i i X Y =，求解A 和B 的正规方程组为：

2

111

1

1()()()n n n

i i i i i i i n n

i i i i A X B X X Y A X nB Y

=====?+=????+=??∑∑∑∑∑ 解出A 和B ，B

C e =.

4.线性最小二乘法的一般形式

一般地，设给定数据组(,)(1,2,)i i x y i n =,01(),(),()m x x x ???为已知的一组[],a b 上线性

无关的函数，选取近似函数为：

0011()()()()m m x a x a x a x ????=++

+

使得： []2

2

110()()n

n

m

i i i k k i i i k y x y a x ??===??

-=-????∑∑∑ （4-1）

[]2

1

min ()n

i i H

i y x ψψ∈==-∑

H 为01(),(),,()m x x x ???的线性组合的全体，这就是线性最小二法的一般形式。

特别地，取()(0,1,,)k k x x k m ?==，就是多项式拟合。

上述问题的正规方程组为:

[]1()()0n

i

i

j

i

i y x x ??=-=∑ (0,1,

,)j m = （4-2）

即：

11

[()()]()m

n

n

k k

i

j

i

i

j

i

k i i a x x y x ???===?=∑∑∑ (0,1,

,)j m = （4-3）

如果记1(,)[()()]n

k j k

i

j

i

k x x ????==

∑，1

(,)()n

j

i

j

i

i y y x ??==∑方程组（4-3）可表示成矩阵形式：

00010001011

11101(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)

(,)(,)m m m m m m m m a y a y a y ???????????????????????????

???????

?????=???????

?????

??????

（4-4） 由()(0,1,

,)k x k m ?=线性无关可导出式（4-4）的系数非奇异，从而保证了方程组的解存在唯一。

5.多项式拟合

对给定的数据组(,)(0,1,

,)i i x y i n =去一个m 次多项式()m n <

2012()m m m x a a x a x a x ?=++++

（5-1） 使得：

2

20121

1

[()](,,,

,)n

n

i

i m i m i i e

y p x F a a a a ===-=∑∑

（5-2） 为最小，即选取参数(0,1,

,)i a i m =，使得

2

201211

(,,,

,)[()]min [()]n

n

m i m i i i H

i i F a a a a y p x y x ψψ∈===-=-∑∑

其中H 为至多m 次多项式集合。这就是数据的m 次多项式拟合，称为这组数据的最小二乘m 次拟合多项式。

由多元函数取得极值的必要条件，得方程组：

10

2()0n m

k j i k i i i k j F

y a x x a ==?=--?=?∑∑ (0,1,

,)i m =

移项得：

1

1

()m

n

n

k j

j k i

i i k i i a x

y x +====∑∑∑

（5-3） 即： 20121111

231

0121

11111220121

1111()()()()()()()()()()()n n n n

m

i i m i i

i i i i n n n n n

m i i i m i i i i i i i i n n n n n m m m m m

i i i m i i i

i i i i i na a x a x a x y a x a x a x a x y x a x a x a x a x y x ====+=====++=====?++++=??

?++++=???++++=??∑∑∑∑∑∑∑∑∑∑∑∑∑∑ （5-4）

这是最小二乘拟合多项式的系数(0,1,

,)k a k m =应满足的方程组，称为正规方程组。由函数组

}{2

1,,,

,m x x x 的线性无关性可证明，方程组（5-4）存在唯一解，且解所对应的多项式（5-1）必定

是已给数据组(,)(1,2,

,)i i x y i n =的最小二乘m 次拟合多项式。

6.矛盾方程组的最小二乘解

设,,m n n m A R x R b R ?∈∈∈方程组为Ax b = （6-1） 其中A 为列满秩矩阵，且方程组（6-1）是矛盾方程组，则

T T Ax b A Ax A b =?=

因为A 列满秩，所以T

A A 正定对称，因而可逆，从而

1()T T x A A A b -=

为矛盾方程组的最小二乘解。

7.求解数据组的最小二乘拟合函数的一般步骤

(1)由给定数据点确定近似函数的表达式，一般可通过描点观察或经验估计得到。

(2)按最小二乘原则确定表达式中的参数，即由偏差平方和最小导出正规方程组，求解参数。 注意：一些简单的非线性最小二乘问题通常需先做变量代换将问题化为线性最小二乘问题再求解。

8.典型应用

（1）已知一组实验数据如下表，求它的拟合曲线？

解： 建立文件w1.m x=[-2,-1,0,1,2]; y=[10,1,0,2,9]; plot(x,y,'o') xlabel('自变量xi') ylabel('函数yi') title('散点图')

画出所给数据(,)i i x y 的散点图

i x -2 -1 0 1 2 i y

10

1

2

9

-2

-1.5

-1

-0.5

00.5

1

1.5

2

012345678910自变量xi

函数y i

散点图

从图可见它像一条抛物线，因而可取抛物线函数2012y a a x a x =++. 将数据带入方程组（5-4）中，得：

0210

251022101

103479

a a a a a +=??

=-??+=? 解得：10.1a =-，2 2.5a =，00.6a =-.

∴拟合曲线为：20.60.1 2.5y x x =--+

（2）设一发射源的发射强度公式形如：00t I I e α-= 现测得与的数据如下表：

i t 0.2 0.3 0.4 0.5 0.6 0.7 0.8 i I

3.16

2.38

1.75

1.34

1.00

0.74

0.56

使用最小二乘法确定0I 与0α.

解：対强度公式两边同取对数得：00ln ln I t I α=-+ 令： ln Y I = ,X t =, 0A α=-, 0ln B I = 得： Y AX B =+

表8-2 数据的代换

i X

0.2 0.3 0.4

0.5 0.6 0.7 0.8

图8-1 数据的散点图

i Y

1.1506 0.8671 0.55962 0.29267 0 -0.30111 -0.57982

将数据带入方程组（1-1）中得：

2.03

3.50.1858

3.57 1.98906

A B A B +=??

+=? 解得： 2.88832A =- , 1.728311B =

∴ 0 2.88832α= , 0exp() 5.631135I B ==

则：强度公式为： 2.8885.631t

I e

-=

(3)某乡镇企业2004-2010年的利润如下表所示，试预测2011和2012年的生产利润？

年份 2004 2005 2006 2007 2008 2009 2010 利润/(万

元)

80

116

144

158

174

196

202

解：由已知数据做一草图，发现该乡镇企业的年生产利润呈直线上升趋势，因此，可用

01()x a a x ?=+作为拟合函数来预测该乡镇企业未来的年生产利润，为简化计算，可把年份记为

2003i i x t =+，相应年份的利润记作，求如表8-3所示数据的线性最小二乘拟合y at b =+.

表8-3 数据的简化

i t 1 2 3 4 5 6 7 i y

80

116

144

158

174

196

202

方法一：

将数据带入方程组（1-1）中得：

140284836

2871070a b a b +=??

+=?

解得： 19.85714a = , 73.42857b =

∴ 拟合曲线为： 19.8571473.42857y t =+

8232.2857y = ， 9252.1428y =

则：2011年的生产利润为232.2857万元，2012年的生产利润为252.1428.

参考文献

[1] John H.Mathews kurtis D.Fink著.数值方法(MATLAB版)(第四版)[M].北京：电子工业出版社，2010.

[2] 杨万利主编.数值分析教程[M].北京：国防工业出版社，2002.

[3] 冯天祥著.数值计算方法理论与实践研究[M].成都：西南交通大学出版社，2005.