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

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

一．曲线拟合 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');

可以绘制出如下图形：

例3.从三次多项式()32

32345P x x x x =-+-上找出21个点(),,1

,2,,21i i x y i = ，然后对这21个点进行“差错处理”，得到新的21个点，根据新的21个点拟合一个新的3次多项式函数，然后和原函数进行比较。

解：在MATLAB 命令窗口输入： p3=inline('2.*x.^3-3.*x.^2+4.*x-5'); x=-10:10; y=p3(x);

e=randn(1,length(x))*80; y=y+e; p=polyfit(x,y,3); xi=-10:0.01:10; yi=polyval(p,xi); plot(xi,yi,x,y, 'o');

hold on

fplot(p3,[-10,10],'r');

5.利用MATLAB的多项式拟合命令polyfit来实现多项式的插值例1.过随机6个数据点，构造5次多项式函数。

解：在MATLAB命令窗口输入：

x=1:6;

y=round(10*randn(1,6));

p=polyfit(x,y,length(x)-1);

xi=1:0.01:6;

yi=polyval(p,xi);

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

可以得到以下图形：

6.利用最小二乘法解超定方程组

例1．解下列超定方程组

24113532627

x y x y x y x y +=??-=?

?

+=??+=? 解：设超定方程的解为00,x x y y ==。

方法一：点()00,x y 到4条直线的距离平方分别为：

()

2

00212

2

241124

x y d +-=

+，()

()

2

00222

2

35335x y d --=

+-，()

2

00232

2

2612

x y d +-=

+，()

2

002422

2721x y d +-=

+

设()4

2001

,i i Q x y d ==∑，根据

00

0,0Q Q

x y ??==??，有： ()()()()000000002223222

24113532627203455x y x y x y x y ???+-+--++-++-=0 ()()()()000000002425222

24113532627203455

x y x y x y x y ???+----++-++-=0 化简有：

00002491299120

1294319450x y x y +-=??

+-=?

解得002.99, 1.30x y == 方法二：最小二乘法：

点()00,x y 关于4条直线的残差平方和为：

()2

002411Q x y =+-()2

00353x y +--()2

0026x y ++-+()2

0027x y +-

根据

00

0,0Q Q

x y ??==??，有： ()()()()00000000424116353226427x y x y x y x y +-+--++-++-=0 ()()()()000000008241110353426227x y x y x y x y +----++-++-=0 化简有：

0000183510

346480x y x y --=??

-+-=?

解得003.04, 1.24x y == 用MATLAB 命令有： syms x0 y0

f1=4*(2*x0+4*y0-11)+2*3*(3*x0-5*y0-3)+2*(x0+2*y0-6)+2*2*(2*x0+y0-7) f2=8*(2*x0+4*y0-11)-2*5*(3*x0-5*y0-3)+4*(x0+2*y0-6)+2*(2*x0+y0-7) 解得： f1 =

36*x0-6*y0-102 f2 f2 =

-6*x0+92*y0-96

继续在MATLAB 命令窗口输入： A=[36,-6,102;-6,92,96]; B=rref(A) x0=B(1,3) y0=B(2,3) 解得： x0 =

3.04029304029304 y0 =

1.24175824175824

方法三：最小二乘法（矩阵运算）

针对方程组Ax b =的最小二乘近似解即为方程组T T A Ax A b =的解 于是，在MATLAB 命令窗口输入： A=[2,4;3,-5;1,2;2,1]; b=[11;3;6;7]; x=inv(A'*A)*A'*b 计算结果为： x =

3.04029304029304

1.24175824175824

方法四：用MATLAB左除命令“\”

在MATLAB命令窗口输入：

A=[2,4;3,-5;1,2;2,1];

b=[11;3;6;7];

x=A\b

即可以得到答案

x =

3.04029304029304

1.24175824175824

可以看出用MATLAB的左除“\”命令计算得到的答案与最小二乘法得到的答案是一致的。其实，MATLAB的左除“\”命令就是按照最小二乘法的原来来编写的。

另外，可以用MATLAB的ezplot命令绘制四条直线的图形

ezplot('2*x+4*y=11');

hold on

ezplot('3*x-5*y=3');

ezplot('x+2*y=6');

ezplot('2*x+y=7');

plot(2.99,1.30,'o');

A=[2,4;3,-5;1,2;2,1];

b=[11;3;6;7];

x=A\b

plot(x(1),x(2),'*');

绘制图形如下：

二．函数逼近

问题，已知函数f(x)，求一个多项式函数()n P x 在区间[a,b]上逼近f(x)。 解决方法：函数的最佳平方逼近。令

()()2

00,,n b

j n j a j Q a a a x f x dx =??=- ???

∑? ，使Q 最小，则有0,0,1,,k Q k n a ?==?

例1.求一次多项式()1P x 在0,2π??

????

上逼近函数()()sin f x x =。 解：构造直线为：y ax b =+，()()()2

20

,sin Q a b ax b x dx π

=+-?，

0,0.Q Q

a b

??==??则有 ()()()()20202sin 0

2sin 0

ax b x xdx ax b x dx π

π

?+-=????+-=??

??，

22220002220

00sin()sin()a x dx b xdx x x dx

a xdx

b dx x dx

πππ

πππ

?+=????+=???????? 32

2

1248

182a b a b ππππ?+=????+=?? 解得：

a=0.6644389，b=0.1147707 在MATLAB 命令窗口输入： xi=0:0.01:pi/2; yi=sin(xi); p=polyfit(xi,yi,1); pi=polyval(p,xi); plot(xi,yi,xi,pi); 可以绘制以下图形：

作业：

（1）用最小二乘法求一个形如2y a bx =+的经验公式，使它与下列数据表拟合。

解：方法一，最小二乘法；方法二，用解超定方程组的思路来解题。

()5

2

21

i

i i Q a bx y ==+-∑，根据

0,0Q Q

a b

??==??，有： ()()5

2

1

5

221

2020i i i i i i i a bx y a bx y x ==?+-=????+-=??∑∑ 在MATLAB 命令窗口输入： x=[19,25,31,38,44]; y=[19,32.3,49,73.3,97.8];

A=[5,sum(x.^2),sum(y);sum(x.^2),sum(x.^4),x.^2*y'];

B=rref(A);

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

xi=min(x):0.01:max(x);

yi=polyval(p,xi);

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

绘制图形如下：

（2）已知数据表如下，试用二次多项式拟合。

（3）求一个形如bx

y ae

（a，b为常数，a>0）的经验公式，使它能和下表数据拟合。

解：公式bx

=可以变为：lny=lna+bx，进一步可以写为Y=A+bx。其中Y=lny，A=lna，

y ae

对应表格为：

（4）求函数y=在区间[1/4,1]上的最小一次式。