数值计算方法教案_曲线拟合与函数逼近

第三章 曲线拟合与函数逼近

一．曲线拟合 1.问题提出：

已知多组数据(),,1,2,

,i i x y i N =，由此预测函数()y f x =的表达式。

数据特点：（1）点数较多。（2）所给数据存在误差。

解决方法：构造一条曲线反映所给数据点的变化总趋势，即所谓的“曲线拟合”。 2.直线拟合的概念 设直线方程为y=a+bx 。

则残差为：ˆi i i e y y

=-,1,2,,i N =，其中ˆi i y

a bx =+。 残差i e 是衡量拟合好坏的重要标志。 可以用MATLAB 软件绘制残差的概念。 x=1:6;

y=[3,4.5,8,10,16,20];

p=polyfit(x,y,1); xi=0:0.01:7; yi=polyval(p,xi); plot(xi,yi,x,y, 'o'); y1=polyval(p,x); hold on for i=1:6

plot([i,i],[y(i),y1(i)], 'r');

end

可以绘制出如下图形：

三个准则： （1）max i e 最小 （2）1n

i i e =∑最小

（3）21

N i i e =∑最小

3.最小二乘法的直线拟合

问题：对于给定的数据点(),,1,2,,i i x y i N =，求一次多项式y=a+bx ，使得总误差Q 最

小。其中()2

21

1

N

N

i i i i i Q e y a bx ====-+⎡⎤⎣⎦∑∑。根据

0,0.Q Q

a b

∂∂==∂∂ 2222

1

222N

i i i i i i i Q y a b x y a y x b x ab =⎡⎤=++--+⎣⎦∑

[]()1

2222N

i i i i i Q a y x b Na y b x a =∂=-+=-+∂∑∑∑ ()221

2222N

i i i i i i i i i Q bx y x x a b x x y a x b =∂⎡⎤=-+=-+⎣⎦∂∑∑∑∑ 故有以下方程组（正则方程）：

2

i i

i

i i i aN b x y a x b x x y ⎧+=⎪⎨+=⎪⎩∑∑∑∑∑ 例1.给定数据表，求最小二乘拟合一次多项式

解：N=5，5

1

i i x =∑=702，5

1

i i y =∑=758，5

2

1

i i x =∑=99864，5

1

i i i x y =∑=108396。

则有方程组

570275870299864108396a b b b +=⎧

⎨

+=⎩

解得a=-60.9392，b=1.5138，则一次多项式为y=-60.9392+1.5138b 用MATLAB 计算并画图如下： x=[165,123,150,123,141]; y=[187,126,172,125,148];

A(1,1)=5;A(1,2)=sum(x);A(1,3)=sum(y); A(2,1)=sum(x);A(2,2)=sum(x.^2);A(2,3)=x*y'; B=rref(A); a=B(1,3);b=B(2,3); p=[b,a];

%以上四行，可以用一行命令 p=polyfit(x,y,1); 替代。 xi=min(x)-1:0.01:max(x)+1; yi=polyval(p,xi); plot(xi,yi,x,y, 'o'); 绘制如下图形

4.最小二乘法的多项式拟合

问题：对于给定的数据点(),,1,2,

,i i x y i N =，求m 次多项式0

m

j j j y a x ==∑（m<

，使得总误差Q 最小。其中2

2

110N

N

m j i i j i i i j Q e y a x ===⎛⎫==- ⎪⎝⎭

∑∑∑。根据0,0,1,,k Q k m a ∂==∂

1020N

m j k i j i i i j y a x x ==⎛⎫

--= ⎪⎝⎭∑∑ 故有正则方程：

011m

N N

j k k j i i i j i i a x y x +===⎛⎫= ⎪⎝⎭∑∑∑ 当m=2时，有

201223

0122342012i i i i i i i i i

i i i i a N a x a x y a x a x a x y x a x a x a x y x ⎧++=⎪++=⎨⎪++=⎩∑∑∑∑∑∑∑∑∑∑∑

例2.求数据表的最小二乘法拟合的二次多项式函数

在MATLAB命令窗口输入：

x=-1:0.25:1;

y=[50,40,25,20,18,21,35,56,66];

A(1,1)=length(x);A(1,2)=sum(x);A(1,3)=sum(x.^2);A(1,4)=sum(y);

A(2,1)=sum(x);A(2,2)=sum(x.^2);A(2,3)=sum(x.^3);A(2,4)=y*x';

A(3,1)=sum(x.^2);A(3,2)=sum(x.^3);A(3,3)=sum(x.^4);A(3,4)=y*(x.^2)'; B=rref(A);

p=[B(3,4),B(2,4),B(1,4)];

%以上五行可以用p=polyfit(x,y,2); 替代

xi=min(x)-0.1:0.01:max(x)+0.1;

yi=polyval(p,xi);

plot(xi,yi,x,y, 'o');

可以绘制出如下图形：