数值曲线拟合最小二乘法课程设计

（一）、曲线拟合的最小二乘法定义：对给定数据 (xi,yi) (i=0,1,…，m)，在取定的函数类中,求p(x),使误差rip(xi)yi(i=0,1,…,m)的平方和最小，即

i0rim2=i02p(x)yiiminm

从几何意义上讲，就是寻求与给定点(xi,yi)(i=0,1,…,m)的距离平方和为最小的曲线yp(x)（图6-1）。函数p(x)称为拟合 函数或最小二乘解，求拟合函数p(x)的方法称为曲线拟合的最小二乘法。

在曲线拟合中，函数类可有不同的选取方法.

（二）、多项式拟合

假设给定数据点(xi,yi)(i=0,1,…,m)，为所有次数不超过n(nm)的多项

pn(x)akxkk0n式构成的函数类，现求一

m,使得

2Ipn(xi)yii02nakxikyimini0k0 (1)

m当拟合函数为多项式时，称为多项式拟合，满足式（1）的pn(x)称为最小二乘拟合多项式。特别地，当n=1时，称为线性拟合或直线拟合。 显然

I(akxikyi)2i0k0mn

为a0,a1,an的多元函数，因此上述问题即为求II(a0,a1,an)的极值 问题。

曲线拟合的最小二乘法 由多元函数求极值的必要条件，得

mnI2(akxikyi)xij0,aji0k0j0,1,,n (2)

即

k0(xi0nmjki)akxijyi,i0mj0,1,,n (3)

（3）是关于a0,a1,an的线性方程组，用矩阵表示为

m1mxii0mxini0xi0mi0mi2xixi0mn1imxyiai0i00mmn1axi1xiyii0i0ammn2nnxiyixii0 (4) i0nim可以证明，方程组（4）的系数矩阵是一个对称正定矩阵，故存在唯一解。从式（4）中解出ak(k=0,1,…，n)，从而可得多项式

pn(x)akxkk0mn (5)

npn(x)为所求的拟合多项式。我们把i0pn(x)的平方误差，记作

r22p(xi)yi2称为最小二乘拟合多项式

pn(xi)yii0m2

由式(2)可得

r22yak(xikyi)2ii0k0i0mnm (6)

（三）、问题的提出

给出一组数据点(xi,yi)列入表7-2中，试用线性最小二乘法求拟合曲线，并估计其误差，作出拟合曲线.

曲线拟合的最小二乘法 表7–2 一组数据(xi,yi)

xi yi -2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6 -192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04 问题的分析

多项式拟合的一般方法可归纳为以下几步：

(1) 由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n；

(2) 列表计算i0xmji(j0,1,,2n)和i0xmjiyi(j0,1,,2n)；

(3) 写出正规方程组，求出a0,a1,an；

pn(x)akxkk0n(4) 写出拟合多项式

。

三、程序及实验结果

（1）在MATLAB工作窗口输入程序：

>> x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];

y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04]; plot(x,y,'r*'),

legend('实验数据(xi,yi)') xlabel('x'), ylabel('y'), title('数据点(xi,yi)的散点图') （2）运行后屏幕显示的散点图如下：

曲线拟合的最小二乘法 （3）编写下列MATLAB程序计算f(x)在(xi,yi)处的函数值，即输入程序

>> syms a1 a2 a3 a4

x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4

运行后屏幕显示关于a1,a2, a3和a4的线性方程组

fi =[ -125/8*a1+25/4*a2-5/2*a3+a4,

-4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4,

a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4,

19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4]

编写构造误差平方和的MATLAB程序：

>> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12

曲线拟合的最小二乘法 6.50 68.04];

fi=[-125/8*a1+25/4*a2-5/2*a3+a4, -4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4,

19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4]; fy=fi-y; fy2=fy.^2; J=sum(fy.^2)

运行后屏幕显示误差平方和如下

J=

(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/100

0*a1+289/100*a2-17/10*a3+a4+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2

为求a1,a2,a3,a4使J达到最小，只需利用极值的必要条件

J0 (k1,2,3,4) ak得到关于a1,a2,a3,a4的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 ：

>> syms a1 a2 a3 a4

J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/10

00*a1+289/100*a2-17/10*a3+a4...+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+

曲线拟合的最小二乘法 a4+843/100)^2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2;

Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3);

Ja4=diff(J,a4);

Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),

Ja41=simple(Ja4),

运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下

Ja11=

56918107/10000*a1+32097579/25000*a2+1377283/25

00*a3+23667/250*a4-8442429/625

Ja21 =

32097579/25000*a1+1377283/2500*a2+23667/250*a3

+67*a4+767319/625

Ja31 =

1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232

638/125

Ja41 =

23667/250*a1+67*a2+18/5*a3+18*a4+14859/25

解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序

>>A=[56918107/10000, 32097579/25000, 1377283/2500,

23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18];

B=[8442429/625, -767319/625, 232638/125, -14859/25]; C=B/A, f=poly2sym(C)

运行后屏幕显示拟合函数f及其系数C如下

C = 5.0911 -14.1905 6.4102 -8.2574 f=716503695845759/140737488355328*x^3 -7988544102557579/562949953421312*x^2 +1804307491277693/281474976710656*x

曲线拟合的最小二乘法 -4648521160813215/562949953421312 故所求的拟合曲线为

f(x)5.0911x314.1905x26.4102x8.2574.

（4）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形.输入程序:

>> xi=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12

6.50 68.04];

n=length(xi);

f=5.0911.*xi.^3-14.1905.*xi.^2+6.4102.*xi -8.2574; x=-2.5:0.01: 3.6;

F=5.0911.*x.^3-14.1905.*x.^2+6.4102.*x -8.2574; fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)

plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off legend('数据点(xi,yi)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),

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

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2

Ew =

3.1054 E1 =

0.9034

曲线拟合的最小二乘法 E2 =

1.2409

以及其数据点(xi,yi)和拟合曲线y=f(x)的图形：

四、结论与总结

通过实验的验证，使我更深刻的理解曲线拟合最小二乘法在数值计算中的重要性，同时，也更深刻的理解了曲线拟合最小二乘法的计算原理。

本次试验我也学会了曲线拟合的方法以及过程，并学会了在计算机上的应用，实现了理论与实践的结合，提高了动手能力。