MATLAB实现非线性曲线拟合最小二乘法

非线性曲线拟合最小二乘法

一、问题提出

设数据（i i y x ,）,(i=0,1,2,3,4).由表3-1给出，表中第四行为i i y y =ln ,可以看出数学模型为bx ae y =,用最小二乘法确定a 及b 。 i

0 1 2 3 4 i x 1.00 1.25 1.50 1.75 2.00 i y

5.10 5.79

6.53

7.45

8.46 i y

1.629

1.756

1.876

2.008

2.135

二、理论基础

根据最小二乘拟合的定义:在函数的最佳平方逼近中],[)(b a C x f ∈，如果f(x)只在一组离散点集{i x ,i=0,1,…,m}，上给定，这就是科学实验中经常见到的实验数据{（i i y x ,）, i=0,1,…,m}的曲线拟合，这里)(i i x f y =，i=0,1,…,m,要求一个函数)(*x S y =与所给数据{（i i y x ,）,i=0,1,…,m}拟合，若记误差

i i i y x S -=)(*δ,i=0,1,…,m,T m ),,(10δδδδ， =,设)(,),(),(10x x x n ??? 是]

,[b a C 上线性无关函数族，在)}(,),(),({10x x x span n ???? =中找一函数)(*x S ，使误差平方和

∑∑∑===∈

-=-==m i m

i m

i i

i

x S i i i

y x S y x S 0

2

)(2

*

2

22

])([])([min

?

δδ

,

这里

)()()()(1100x a x a x a x S n n ???+++= (n

这就是一般的最小二乘逼近，用几何语言来说，就称为曲线拟合的最小二乘

法。

在建模的过程中应用到了求和命令（sum ）、求偏导命令(diff)、化简函数命令（simple ）、用迭代方法解二元非线性方程组的命令（fsolve ），画图命令（plot ）等。

三、实验内容

用最小二乘法求拟合曲线时，首先要确定S(x)的形式。这不单纯是数学问题，

还与所研究问题的运动规律及所得观测数据（i i y x ,）有关；通常要从问题的运动规律及给定数据描图，确定S(x)的形式，并通过实际计算选出较好的结果。S(x)的一般表达式为线性形式，若)(x k ?是k 次多项式，S(x)就是n 次多项式，为了使问题的提法更有一般性，通常在最小二乘法中22

δ

都考虑为加权平方和

∑=-=m

i i i i x f x S x 0

222)]()()[(ωδ

.

这里0)(≥x ω是[a,b]上的权函数，它表示不同点))(,(i i x f x 处的数据比重不同。

用最小二乘法求拟合曲线的问题，就是求形如S(x)的一个函数)(*x S y =，使∑=-=m

i i i i x f x S x 0

222

)]()()[(ωδ

取得最小。它转化为求多元函数

20

10)]()([)(),,,(i n

j i j j m i i n x f x a x a a a I -=∑∑==?ω

的极小点),,,(**1*0n a a a 问题。再由求多元函数极值的必要条件，有

∑∑===-=??m i i k n

j i i j j i k x x f x a x a I

00

0)(])()()[(2??ω (k=0,1,…,n) 此题中假设1)(=x ω，由已知所给数据点（i i y x ,）画出图形，根据离散点的位置观察出它们所拟合的曲线图形应类似于指数函数的曲线图形，故设拟合曲线的函数为bx ae y =。

本题编程过程中，令f=y,z1=a,z2=b,

令拟合曲线中对应i x 的函数值与i y 的差的平方和为J,即J=sum(fy.^2);分别求J 关于z1,z2的偏导，简化后并令其分别为0得一关于z1,z2的二元非线性方程组，

最后利用fsolve 命令求得z1,z2的值分别为

z1=3.0751 z2=0.5052

故得到拟合曲线为

x e y 5052.00751.3=

为证明曲线拟合的正确性，我们将离散点（i i y x ,）与所得的拟合曲线

x e y 5052.00751.3=

画于同一图形中，图形如下：

1 1.1 1.

2 1.

3 1.4

1.5 1.6 1.7 1.8 1.92

5

5.5

6

6.5

77.5

8

8.5

x 轴

y 轴

例8的数据点(x(i),y(i))和拟合曲线y=f(x)的图形

数据点(x(i),y(i))拟合曲线y=f(x)

四、结果分析

根据实验内容求得拟合曲线bx ae y =中未知数a,b 分别为

a=3.0751 b=0.5052

即拟合曲线为

x e y 5052.00751.3=。

由图形知拟合成功！

参考文献

1·《数值分析》，李庆扬，王能超，易大义，2001，清华大学出版社（第

四版）。 2·《数值方法》，关治，陆金甫，2006，清华大学出版社。 3·《数值分析与实验学习指导》，蔡大用，2001，清华大学出版社。 4·《数值分析与实验》，薛毅，2005，北京工业大学出版社.

附录

程序1：syms z1 z2

x=1.00:0.25:2.00;

y=[5.10,5.79,6.53,7.45,8.46];

f=z1*exp(z2.*x)

fy=f-y;

J=sum(fy.^2);

Ja=diff(J,z1);

Jb=diff(J,z2);

Ja1=simple(Ja),

Jb1=simple(Jb),

程序2：function [y1,y2]=fun(z)

y1=2*z(1)*exp(2*z(2))-51/5*exp(z(2))+2*z(1)*exp(5/2*z(2))-579/50*exp(

5/4*z(2))+2*z(1)*exp(3*z(2))-653/50*exp(3/2*z(2))+2*z(1)*exp(7/2*z(2))

-149/10*exp(7/4*z(2))+2*exp(4*z(2))*z(1)-423/25*exp(2*z(2));

y2=-1/200*z(1)*(-400*z(1)*exp(2*z(2))+2040*exp(z(2))-500*z(1)*exp(5/2

*z(2))+2895*exp(5/4*z(2))-600*z(1)*exp(3*z(2))+3918*exp(3/2*z(2))-700

*z(1)*exp(7/2*z(2))+5215*exp(7/4*z(2))-800*z(1)*exp(4*z(2))+6768*exp

(2*z(2)));

作图程序：x=1.00:0.25:2.00;

y=[5.10,5.79,6.53,7.45,8.46];

f=3.0751*exp(0.5052*x);

plot(x,y,'r*',x,f,'b-');

xlabel('x轴'),ylabel('y轴'),

legend('数据点(x(i),y(i))','拟合曲线y=f(x)'),

title('例8的数据点(x(i),y(i))和拟合曲线y=f(x)的图形')