第三章 数据拟合
? 用插值的方法对一函数进行近似,要求所得到的
插值多项式经过已知插值节点 ;在 n比较大的情
况下,插值多项式往往是高次多项式,这也就容
易出现振荡现象（龙格现象），即虽然在插值
节点上没有误差,但在插值节点之外插值误差变
得很大,从“整体”上看,插值逼近效果将变得
“很差”。
? 所谓数据拟合是求一个简单的函数,例如是一个
低次多项式,不要求通过已知的这些点,而是要
求在整体上, 尽量好, 的逼近原函数 。 这时,在
每个已知点上就会有误差,数据拟合就是从整体
上使误差,尽量的小一些 。
3.1 多项式拟合
? n次多项式,
? 曲线与数据点的残差为,
? 残差的平方和为,
? 为使其最小化，可令 R关于 的偏导数为零，
即,
1
1 2 1()
nn
ng x c x c x c
?
?? ? ? ?L
( ),1,2,,i i ir y g x i L? ? ? L
2
1
L
i
i
Rr
?
? ?
jc
0,1,2,,1
j
R jn
c
? ? ? ?
?
L
? 或
? 或矩阵形式,
1
2 2 1
1 1 1
( ),1,2,,1
n L L
n j k n k
i j i i
j i i
x c x y k n
?
? ? ? ? ?
? ? ?
? ? ?? ? ? L
2 2 1
1 1 1 1
1
2 1 1 1
2
1 1 1
1
0
1 1 1
.
..
.
.,,,,
..
L L L L
n n n n
i i i i i
i i i i
L L L
n n n
i i i i
i i i
n
L L L
n
i i i
i i i
x x x x y
c
x x x yc
c
x x y
?
? ? ? ?
? ? ?
? ? ?
?
? ? ?
? ? ? ?
? ? ? ?
? ? ? ???
? ? ? ???
? ? ? ???
?
? ? ? ???
? ? ? ???
??? ? ? ???
? ? ? ?
? ? ? ?
? ? ? ?
? ? ? ?
? ? ?
? ? ?
多项式拟合 MATLAB命令,polyfit
格式,p=polyfit(x,y,n)
>> x0=0:.1:1; y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
>> p3=polyfit(x0,y0,3); vpa(poly2sym(p3),10)
% 可以如下显示多项式
ans =
2.839962923*x^3-
4.789842696*x^2+1.943211631*x+.5975248921e-1
例
? 绘制拟合曲线,
>> x=0:.01:1; ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> y1=polyval(p3,x); plot(x,y1,x,ya,x0,y0,'o')
? 就不同的次数进行拟合,
>> p4=polyfit(x0,y0,4); y2=polyval(p4,x);
>> p5=polyfit(x0,y0,5); y3=polyval(p5,x);
>> p8=polyfit(x0,y0,8); y4=polyval(p8,x);
>> plot(x,ya,x0,y0,'o',x,y2,x,y3,x,y4)
? 拟合最高次数为 8的多项式,
>> vpa(poly2sym(p8),5)
ans =
-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-
125.29*x^4+74.450*x^3-
27.672*x^2+4.9869*x+.42037e-6
? Taylor幂级数展开,
>> syms x; y=(x^2-3*x+5)*exp(-5*x)*sin(x);
>> vpa(taylor(y,9),5)
ans =
5.*x-28.*x^2+77.667*x^3-142.*x^4+192.17*x^5-
204.96*x^6+179.13*x^7-131.67*x^8
? 多项式表示数据模型是不唯一的，即是两个多项式
函数完全不同。在某一区域内其曲线将特别近似。
多项式拟合的效果并不一定总是很精确的。
>> x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2);
>> x=-1:.01:1; ya=1./(1+25*x.^2);
>> p3=polyfit(x0,y0,3);
>> y1=polyval(p3,x);
>> p5=polyfit(x0,y0,5);
>> y2=polyval(p5,x);
>> p8=polyfit(x0,y0,8);
>> y3=polyval(p8,x);
>> p10=polyfit(x0,y0,10);
>> y4=polyval(p10,x);
>> plot(x,ya,x,y1,x,y2,'-.',x,y3,'--',x,y4,':')
例
? 用 Taylor幂级数展开效果将更差。
>> syms x; y=1/(1+25*x^2); p=taylor(y,x,10)
p =
1-25*x^2+625*x^4-15625*x^6+390625*x^8
? 多项式拟合效果
>> x1=-1:0.01:1;
>> ya=1./(1+25*x1.^2);
>> y1=subs(p,x,x1);
>> plot(x1,ya,'--‘,x1,y1)
3.2 函数线性组合的曲线拟合方法
该方程的最小二乘解为,
其中
例
>> x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
>> y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;
2.1979;2.5409;2.9627;3.155;3.2052];
>> A=[ones(size(x)),exp(-3*x),cos(-2*x).*exp(-
4*x), x.^2];
>> c=A\y; c1=c'
c1 =
1.2200 2.3397 -0.6797 0.8700
? 图形显示
>> x0=[0:0.01:1.5]';
>> A1=[ones(size(x0)) exp(-3*x0),cos(-
2*x0).*exp(-4*x0) x0.^2];
>> y1=A1*c;
>> plot(x0,y1,x,y,'x')
? 数据分析
>> x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,..,
2.2255,2.4596,2.7183,3.6693];
>> y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,..,
0.2864,0.2532,0.2238,0.1546];
>> plot(x,y,x,y,'*')
例
? 分别对 x,y进行对数变换,
>> x1=log(x); y1=log(y); plot(x1,y1)
>> A=[x1',ones(size(x1'))]; c=[A\y1']‘
c =
-1.2339 -0.2630
>> exp(c(2))
ans =
0.7687
>> x=[0:0.1:1]'; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); n=8;
A=[];
>> for i=1:n+1,A(:,i)=x.^(n+1-i); end
>> c=A\y; vpa(poly2sym(c),5)
ans =
-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-
125.29*x^4+74.450*x^3-
27.672*x^2+4.9869*x+.42037e-6
例
3.3 最小二乘曲线拟合
? 格式,
[a,jm]=lsqcurvefit(Fun,a0,x,y)
例
>> x=0:.1:10;
>> y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
>> f=inline('a(1)*exp(-a(2)*x)+a(3)*…
exp(-a(4)*x).*sin(a(5)*x)','a','x');
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y); xx',res
Optimization terminated successfully,
Relative function value changing by less than
OPTIONS.TolFun
ans =
0.1197
0.2125
0.5404
0.1702
1.2300
res =
7.1637e-007
修改最优化选项,
>> ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff); xx',res
Optimization terminated successfully,
Relative function value changing by less than OPTIONS.TolFun
ans =
0.1200
0.2130
0.5400
0.1700
1.2300
res =
9.5035e-021
? 绘制曲线,
>> x1=0:0.01:10; y1=f(xx,x1); plot(x1,y1,x,y,'o')
例
>> x=0.1:0.1:1;
>> y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,
4.8232,5.1275];
function y=c8f3(a,x)
y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);
>> a=lsqcurvefit('c8f3',[1;2;2;3],x,y); a'
Maximum number of function evaluations exceeded;
increase options.MaxFunEvals
ans =
2.4575 2.4557 1.4437 2.0720
? 绘制曲线,
>> y1=c8f3(a,x); plot(x,y,x,y1)
8.4 B样条函数及其 MATLAB表示
? 格式
S=spapi(k,x,y)
例
>> x0=[0,0.4,1,2,pi];
>> y0=sin(x0);
>> ezplot('sin(t)',[0,pi]);
>> hold on
>> sp1=csapi(x0,y0);
>> fnplt(sp1,'--');
% 三次分段多项式样条插值
>> sp2=spapi(5,x0,y0); fnplt(sp2,':') % 5 次 B 样条插值
y=sin(t)
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> ezplot('(x^2-3*x+5)*exp(-5*x)*sin(x)',[0,1]),hold on
>> sp1=csapi(x,y); fnplt(sp1,'--'); sp2=spapi(5,x,y);
fnplt(sp2,':')
? 用插值的方法对一函数进行近似,要求所得到的
插值多项式经过已知插值节点 ;在 n比较大的情
况下,插值多项式往往是高次多项式,这也就容
易出现振荡现象（龙格现象），即虽然在插值
节点上没有误差,但在插值节点之外插值误差变
得很大,从“整体”上看,插值逼近效果将变得
“很差”。
? 所谓数据拟合是求一个简单的函数,例如是一个
低次多项式,不要求通过已知的这些点,而是要
求在整体上, 尽量好, 的逼近原函数 。 这时,在
每个已知点上就会有误差,数据拟合就是从整体
上使误差,尽量的小一些 。
3.1 多项式拟合
? n次多项式,
? 曲线与数据点的残差为,
? 残差的平方和为,
? 为使其最小化，可令 R关于 的偏导数为零，
即,
1
1 2 1()
nn
ng x c x c x c
?
?? ? ? ?L
( ),1,2,,i i ir y g x i L? ? ? L
2
1
L
i
i
Rr
?
? ?
jc
0,1,2,,1
j
R jn
c
? ? ? ?
?
L
? 或
? 或矩阵形式,
1
2 2 1
1 1 1
( ),1,2,,1
n L L
n j k n k
i j i i
j i i
x c x y k n
?
? ? ? ? ?
? ? ?
? ? ?? ? ? L
2 2 1
1 1 1 1
1
2 1 1 1
2
1 1 1
1
0
1 1 1
.
..
.
.,,,,
..
L L L L
n n n n
i i i i i
i i i i
L L L
n n n
i i i i
i i i
n
L L L
n
i i i
i i i
x x x x y
c
x x x yc
c
x x y
?
? ? ? ?
? ? ?
? ? ?
?
? ? ?
? ? ? ?
? ? ? ?
? ? ? ???
? ? ? ???
? ? ? ???
?
? ? ? ???
? ? ? ???
??? ? ? ???
? ? ? ?
? ? ? ?
? ? ? ?
? ? ? ?
? ? ?
? ? ?
多项式拟合 MATLAB命令,polyfit
格式,p=polyfit(x,y,n)
>> x0=0:.1:1; y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
>> p3=polyfit(x0,y0,3); vpa(poly2sym(p3),10)
% 可以如下显示多项式
ans =
2.839962923*x^3-
4.789842696*x^2+1.943211631*x+.5975248921e-1
例
? 绘制拟合曲线,
>> x=0:.01:1; ya=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> y1=polyval(p3,x); plot(x,y1,x,ya,x0,y0,'o')
? 就不同的次数进行拟合,
>> p4=polyfit(x0,y0,4); y2=polyval(p4,x);
>> p5=polyfit(x0,y0,5); y3=polyval(p5,x);
>> p8=polyfit(x0,y0,8); y4=polyval(p8,x);
>> plot(x,ya,x0,y0,'o',x,y2,x,y3,x,y4)
? 拟合最高次数为 8的多项式,
>> vpa(poly2sym(p8),5)
ans =
-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-
125.29*x^4+74.450*x^3-
27.672*x^2+4.9869*x+.42037e-6
? Taylor幂级数展开,
>> syms x; y=(x^2-3*x+5)*exp(-5*x)*sin(x);
>> vpa(taylor(y,9),5)
ans =
5.*x-28.*x^2+77.667*x^3-142.*x^4+192.17*x^5-
204.96*x^6+179.13*x^7-131.67*x^8
? 多项式表示数据模型是不唯一的，即是两个多项式
函数完全不同。在某一区域内其曲线将特别近似。
多项式拟合的效果并不一定总是很精确的。
>> x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2);
>> x=-1:.01:1; ya=1./(1+25*x.^2);
>> p3=polyfit(x0,y0,3);
>> y1=polyval(p3,x);
>> p5=polyfit(x0,y0,5);
>> y2=polyval(p5,x);
>> p8=polyfit(x0,y0,8);
>> y3=polyval(p8,x);
>> p10=polyfit(x0,y0,10);
>> y4=polyval(p10,x);
>> plot(x,ya,x,y1,x,y2,'-.',x,y3,'--',x,y4,':')
例
? 用 Taylor幂级数展开效果将更差。
>> syms x; y=1/(1+25*x^2); p=taylor(y,x,10)
p =
1-25*x^2+625*x^4-15625*x^6+390625*x^8
? 多项式拟合效果
>> x1=-1:0.01:1;
>> ya=1./(1+25*x1.^2);
>> y1=subs(p,x,x1);
>> plot(x1,ya,'--‘,x1,y1)
3.2 函数线性组合的曲线拟合方法
该方程的最小二乘解为,
其中
例
>> x=[0,0.2,0.4,0.7,0.9,0.92,0.99,1.2,1.4,1.48,1.5]';
>> y=[2.88;2.2576;1.9683;1.9258;2.0862;2.109;
2.1979;2.5409;2.9627;3.155;3.2052];
>> A=[ones(size(x)),exp(-3*x),cos(-2*x).*exp(-
4*x), x.^2];
>> c=A\y; c1=c'
c1 =
1.2200 2.3397 -0.6797 0.8700
? 图形显示
>> x0=[0:0.01:1.5]';
>> A1=[ones(size(x0)) exp(-3*x0),cos(-
2*x0).*exp(-4*x0) x0.^2];
>> y1=A1*c;
>> plot(x0,y1,x,y,'x')
? 数据分析
>> x=[1.1052,1.2214,1.3499,1.4918,1.6487,1.8221,2.0138,..,
2.2255,2.4596,2.7183,3.6693];
>> y=[0.6795,0.6006,0.5309,0.4693,0.4148,0.3666,0.3241,..,
0.2864,0.2532,0.2238,0.1546];
>> plot(x,y,x,y,'*')
例
? 分别对 x,y进行对数变换,
>> x1=log(x); y1=log(y); plot(x1,y1)
>> A=[x1',ones(size(x1'))]; c=[A\y1']‘
c =
-1.2339 -0.2630
>> exp(c(2))
ans =
0.7687
>> x=[0:0.1:1]'; y=(x.^2-3*x+5).*exp(-5*x).*sin(x); n=8;
A=[];
>> for i=1:n+1,A(:,i)=x.^(n+1-i); end
>> c=A\y; vpa(poly2sym(c),5)
ans =
-8.2586*x^8+43.566*x^7-101.98*x^6+140.22*x^5-
125.29*x^4+74.450*x^3-
27.672*x^2+4.9869*x+.42037e-6
例
3.3 最小二乘曲线拟合
? 格式,
[a,jm]=lsqcurvefit(Fun,a0,x,y)
例
>> x=0:.1:10;
>> y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);
>> f=inline('a(1)*exp(-a(2)*x)+a(3)*…
exp(-a(4)*x).*sin(a(5)*x)','a','x');
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y); xx',res
Optimization terminated successfully,
Relative function value changing by less than
OPTIONS.TolFun
ans =
0.1197
0.2125
0.5404
0.1702
1.2300
res =
7.1637e-007
修改最优化选项,
>> ff=optimset; ff.TolFun=1e-20; ff.TolX=1e-15; % 修改精度限制
>> [xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff); xx',res
Optimization terminated successfully,
Relative function value changing by less than OPTIONS.TolFun
ans =
0.1200
0.2130
0.5400
0.1700
1.2300
res =
9.5035e-021
? 绘制曲线,
>> x1=0:0.01:10; y1=f(xx,x1); plot(x1,y1,x,y,'o')
例
>> x=0.1:0.1:1;
>> y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,
4.8232,5.1275];
function y=c8f3(a,x)
y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);
>> a=lsqcurvefit('c8f3',[1;2;2;3],x,y); a'
Maximum number of function evaluations exceeded;
increase options.MaxFunEvals
ans =
2.4575 2.4557 1.4437 2.0720
? 绘制曲线,
>> y1=c8f3(a,x); plot(x,y,x,y1)
8.4 B样条函数及其 MATLAB表示
? 格式
S=spapi(k,x,y)
例
>> x0=[0,0.4,1,2,pi];
>> y0=sin(x0);
>> ezplot('sin(t)',[0,pi]);
>> hold on
>> sp1=csapi(x0,y0);
>> fnplt(sp1,'--');
% 三次分段多项式样条插值
>> sp2=spapi(5,x0,y0); fnplt(sp2,':') % 5 次 B 样条插值
y=sin(t)
>> x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
>> ezplot('(x^2-3*x+5)*exp(-5*x)*sin(x)',[0,1]),hold on
>> sp1=csapi(x,y); fnplt(sp1,'--'); sp2=spapi(5,x,y);
fnplt(sp2,':')
