曲线拟合的线性最小二乘法及其MATLAB程序

3.1 曲线拟合的线性最小二乘法及其MATLAB程序

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

表3-1 例3.1.1的一组数据(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）在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）编写下列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, 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*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 105．

得到关于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/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+(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/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

解线性方程组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('例3.1.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.

Ew = E1 = E2 =

3.105 4 0.903 4 1.240 9

106．

3.2 函数rk(x)的选取及其MATLAB程序

例3.2.1 给出一组实验数据点(xi,yi)的横坐标向量为

x=(-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6)，纵横坐标向量为y=(459.26,52.81,198.27,165.60,59.17,41.66,25.92, 22.37,13.47, 12.87, 11.87,6.69,14.87,24.22)，试用线性最小二乘法求拟合曲线，并估计其误差，作出拟合曲线.

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

>>x=[-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6];

y=[459.26,52.81,198.27,165.60,59.17,41.66,25.92,

22.37,13.47, 12.87, 11.87,6.69,14.87,24.22];

plot(x,y,'r*'),legend('实验数据(xi,yi)') xlabel('x'), ylabel('y'),

title('例3.2.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（3）编写下列MATLAB程序计算f(x)在(xi,yi)处的函数值，即输入程序

>> syms a b

x=[-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5,-2.7,-3.6]; fi=a.*exp(-b.*x)

运行后屏幕显示关于a和b的线性方程组

fi =

[ a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b),

a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(18/5*b), a*exp(17/5*b), a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)]

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

>>y=[459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47,12.87, 11.87, 6.69,14.87,24.22];

fi =[ a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b), a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(18/5*b), a*exp(17/5*b), a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)]; fy=fi-y; fy2=fy.^2; J=sum(fy.^2)

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

J =

(a*exp(17/2*b)-22963/50)^2+(a*exp(87/10*b)-5281/100)^2+(

a*exp(71/10*b)-19827/100)^2+(a*exp(34/5*b)-828/5)^2+(a*exp(51/10*b)-5917/100)^2+(a*exp(9/2*b)-2083/50)^2+(a*exp(18/5*b)-648/25)^2+(a*exp(17/5*b)-2237/100)^2+(a*exp(13/5*b)-1347/100)^2+(a*exp(5/2*b)-1287/100)^2+(a*exp(21/10*b)-1187/100)^2+(a*exp(3/2*b)-669/100)^2+(a*exp(27/10*b)-1487/100)^2+(a*exp(18/5*b)-1211/50)^2

为求a,b使J达到最小，只需利用极值的必要条件，得到关于a,b的线性方程组，这可以由下面的MATLAB程序完成，即输入程序

>> syms a b

J=(a*exp(17/2*b)-22963/50)^2+(a*exp(87/10*b)-5281/100)^2

+(a*exp(71/10*b)-19827/100)^2+(a*exp(34/5*b)-828/5)^2+(a*exp(51/10*b)-5917/100)^2+(a*exp(9/2*b)-2083/50)^2+(a*exp(18/5*b)-648/25)^2+(a*exp(17/5*b)-2237/100)^2+(a*exp(13/5*b)-1347/100)^2+(a*exp(5/2*b)-1287/100)^2+(a*exp(21/10*b)-1187/100)^2+(a*exp(3/2*b)-669/100)^2+(a*exp(27/10*b)-1487/100)^2+(a*exp(18/5*b)-1211/50

107．

)^2;

Ja=diff(J,a); Jb=diff(J,b); Ja1=simple(Ja), Jb1=simple(Jb),

运行后屏幕显示J分别对a,b的偏导数如下

Ja1 =

2*a*exp(3*b)+2*a*exp(17*b)+2*a*exp(87/5*b)+2*exp(68/5*b)

*a+2*exp(9*b)*a+2*a*exp(34/5*b)-669/50*exp(3/2*b)-1487/50*exp(27/10*b)-2507/25*exp(18/5*b)-22963/25*exp(17/2*b)-5281/50*exp(87/10*b)-19827/50*exp(71/10*b)-2237/50*exp(17/5*b)-1656/5*exp(34/5*b)-1347/50*exp(13/5*b)-5917/50*exp(51/10*b)-1287/50*exp(5/2*b)-2083/25*exp(9/2*b)-1187/50*exp(21/10*b)+4*a*exp(36/5*b)+2*a*exp(26/5*b)+2*a*exp(71/5*b)+2*a*exp(51/5*b)+2*a*exp(5*b)+2*a*exp(21/5*b)+2*a*exp(27/5*b)

Jb1 =

1/500*a*(2100*a*exp(21/10*b)^2+8500*a*exp(17/2*b)^2+6800

*a*exp(34/5*b)^2-10035*exp(3/2*b)-40149*exp(27/10*b)-180504*exp(18/5*b)-3903710*exp(17/2*b)-459447*exp(87/10*b)-1407717*exp(71/10*b)-76058*exp(17/5*b)-1126080*exp(34/5*b)-35022*exp(13/5*b)-301767*exp(51/10*b)-32175*exp(5/2*b)-187470*exp(9/2*b)-24927*exp(21/10*b)+7100*a*exp(71/10*b)^2+5100*a*exp(51/10*b)^2+4500*a*exp(9/2*b)^2+7200*a*exp(18/5*b)^2+3400*a*exp(17/5*b)^2+2600*a*exp(13/5*b)^2+2500*a*exp(5/2*b)^2+1500*a*exp(3/2*b)^2+2700*a*exp(27/10*b)^2+8700*a*exp(87/10*b)^2)

用解二元非线性方程组的牛顿法的MATLAB程序求解线性方程组Ja1 =0，Jb1 =0，得

a = b=

2.811 0 0.581 6

故所求的拟合曲线（7.13）为

f(x)?2.8110e?0.5816x.

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

>> xi=[-8.5 -8.7 -7.1 -6.8 -5.10 -4.5 -3.6 -3.4 -2.6 -2.5

-2.1 -1.5 -2.7 -3.6];

y=[459.26 52.81 198.27 165.60 59.17 41.66 25.92 22.37

13.47 12.87 11.87 6.69 14.87 24.22];

n=length(xi); f=2.8110.*exp(-0.5816.*xi); x=-9:0.01: -1; F=2.8110.*exp(-0.5816.*x); 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.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew = 390.141 5，平均误差E1=36.942 2和均方根误差E2=106.031 7及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.

3.3 多项式拟合及其MATLAB程序

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

表3–3 例3.3.1的一组数据(xi,yi) xi yi -2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6 53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88 解 （1）首先根据表3–3给出的数据点(xi,yi)，用下列MATLAB程序画出散点图.

108．

在MATLAB工作窗口输入程序

>> x=[-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; y=[53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12

19.88];

plot(x,y,'r*'), legend('数据点(xi,yi)') xlabel('x'), ylabel('y'),

title('例3.3.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（3）用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数ak

(k?1,2,3).输入程序

>> a=polyfit(x,y,2)

运行后输出（7.16）式的系数

a =

2.8302 -7.3721 9.1382

故拟合多项式为

f(x)?2.8302x2?7.3721x?9.1382.

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

>> xi=[-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6];

y=[53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88]; n=length(xi); f=2.8302.*xi.^2-7.3721.*xi+9.1382 x=-2.9:0.001:3.6;F=2.8302.*x.^2-7.3721.*x+8.79; fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n), plot(xi,y,'r*', x,F,'b-'), legend('数据点(xi,yi)','拟合曲线y=f(x)') xlabel('x'), ylabel('y'),

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

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.

Ew = E1 = E2 =

0.745 7, 0.389 2, 0.436 3

3.4 拟合曲线的线性变换及其MATLAB程序

例3．4.1 给出一组实验数据点(xi,yi)的横坐标向量为x=（7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6），纵横坐标向量为y=（359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22），试用线性变换和线性最小二乘法求拟合曲线，并估计其误差，作出拟合曲线.

解 （1）首先根据给出的数据点(xi,yi)，用下列MATLAB程序画出散点图. 在MATLAB工作窗口输入程序

>> x=[7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7

3.6];

y=[359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87

11.87 6.69 14.87 24.22];

plot(x,y,'r*'), legend('数据点(xi,yi)') xlabel('x'), ylabel('y'),

title('例3.4.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（2）根据数据散点图，取拟合曲线为

bx y?ae (a?0,b?0),

其中a,b是待定系数.令Y?lny,A?lna,B?b，则（7.19）化为Y?A?Bx.在MATLAB工作窗口输入程序

>> x=[7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7

3.6];

109．

y=[359.26 165.60 59.17 41.66 25.92 22.37 13.47

12.87 11.87 6.69 14.87 24.22];

Y=log(y); a=polyfit(x,Y,1); B=a(1);A=a(2); b=B,a=exp(A) n=length(x); X=8:-0.01:1; Y=a*exp(b.*X); f=a*exp(b.*x); plot(x,y,'r*',X,Y,'b-'), xlabel('x'),ylabel('y') legend('数据点(xi,yi)','拟合曲线y=f(x)')

title('例3.4.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形') fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)

bx运行后屏幕显示y?ae的系数b =0.624 1，a =2.703 9,数据(xi,yi)与拟合函数f的最大误差Ew =67.641 9，平均误差E1=8.677 6和均方根误差E2=20.711 3及其数据点(xi,yi)和拟合曲线f(x)?2.7039e

0.6241x的图形（略）.

3.5 函数逼近及其MATLAB程序

最佳均方逼近的MATLAB主程序

function [yy1,a,WE]=zjjfbj(f,X,Y,xx)

m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1); if n~=length(Y)

error('X和Y的维数应该相同') end

for j=1:m for k=1:m b(j,k)=0; for i=1:n

b(j,k)=b(j,k)+feval(f(j,:),X(i))*feval(f(k,:),X(i));

end end

c(j)=0; for i=1:n

c(j)=c(j)+feval(f(j,:),X(i))*Y(i); end end a=b\\c; WE=0;

for i=1:n ff=0;

for j=1:m

ff=ff+a(j)*feval(f(j,:),X(i)); end

WE=WE+(Y(i)-ff)*(Y(i)-ff); end

if nargin==3 return; end yy=[]; for i=1:m l=[];

for j=1:length(xx)

l=[l,feval(f(i,:),xx(j))]; end

yy=[yy l']; end

yy=yy*a; yy1=yy'; a=a';WE;

2例3.5.1 对数据X和Y, 用函数y?1,y?x,y?x进行逼近，用所得到的逼近函

110．

数计算在x?6.5处的函数值，并估计误差.其中

X=(1 3 4 5 6 7 8 9); Y=(-11 -13 -11 -7 -1 7 17 29). 解 在MATLAB工作窗口输入程序

>> X=[ 1 3 4 5 6 7 8 9]; Y=[-11 -13 -11 -7 -1 7 17

29];

f=['fun0';'fun1';'fun2']; [yy,a,WE]=zjjfbj(f,X,Y,6.5)

运行后屏幕显示如下

yy =

2.75000000000003 a =

-7.00000000000010 -4.99999999999995 1.00000000000000 WE =

7.172323350269439e-027

y?cosx，y?e，例3.5.2 对数据X和Y，用函数y?1,y?x,y?x2，y?sinx进行逼近，其中X=（0 0.50 1.00 1.50 2.00 2.50 3.00），Y=（0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645）.

解 在MATLAB工作窗口输入程序

>> X=[ 0 0.50 1.00 1.50 2.00 2.50 3.00];

Y=[0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645];

f=['fun0';'fun1';'fun2';'fun3';'fun4';'fun5'];xx=0:0.2:3; [yy,a,WE]=zjjfbj(f,X,Y, xx), plot(X,Y,'ro',xx,yy,'b-')

运行后屏幕显示如下（图略）

yy = Columns 1 through 7

-0.0005 0.2037 0.3939 0.5656 0.7141 0.8348

0.9236

Columns 8 through 14

0.9771 0.9926 0.9691 0.9069 0.8080 0.6766

0.5191

Columns 15 through 16 0.3444 0.1642

a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653 WE = 1.5769e-004

即，最佳逼近函数为

y=0.3828+0.4070*x-0.3901*x^2+0.0765*exp(x) -0.4598*cos(x) +0.5653*sin(x).

x 111．

本文来源：https://www.bwwdw.com/article/q377.html