3.1 曲线拟合的线性最小二乘法及其MATLAB 程序
例3.1．1 给出一组数据点),(i i y x 列入表3-1中，试用线性最小二乘法求拟合曲线，并估计其误差，作出拟合曲线.
表3-1 例3.1.1的一组数据),(y x
解 （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('例3.1.1的数据点(xi,yi)的散点图')
运行后屏幕显示数据的散点图（略）.
（3）编写下列MA TLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序
>> 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
运行后屏幕显示关于a 1,a 2, a 3和a 4的线性方程组
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, 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/1000*a1+2
89/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*a 2+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
为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0=∂∂k
a J )4,3,2,1(=k ，
得到关于4321,,,a a a a 的线性方程组，这可以由下面的MA TLAB 程序完成，即输入程序
>> syms a1 a2 a3 a4
J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*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+(a 4+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 分别对a 1, a 2 ,a 3 ,a 4的偏导数如下
Ja11=
56918107/10000*a1+32097579/25000*a2+1377283/2500*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-232638/125
Ja41 =
23667/250*a1+67*a2+18/5*a3+18*a4+14859/25
解线性方程组Ja 11 =0，Ja 21 =0，Ja 31 =0，Ja 41 =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
故所求的拟合曲线为
8.25746.410214.19055.0911)(23-+-=x x x x f .
（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('例3.1.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')
运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w ，平均误差E 1和均方根误差E 2及其数据点),(i i y x 和拟合曲线y =f (x )的图形（略）.
Ew = E1 = E2 =
3.105 4 0.903 4 1.240 9
