数值计算方法教案曲线拟合与函数逼近.docx
《数值计算方法教案曲线拟合与函数逼近.docx》由会员分享,可在线阅读,更多相关《数值计算方法教案曲线拟合与函数逼近.docx(15页珍藏版)》请在冰豆网上搜索。
数值计算方法教案曲线拟合与函数逼近
第三章曲线拟合与函数逼近
一.曲线拟合
1.问题提出:
已知多组数据,由此预测函数的表达式。
数据特点:
(1)点数较多。
(2)所给数据存在误差。
解决方法:
构造一条曲线反映所给数据点的变化总趋势,即所谓的“曲线拟合”。
2.直线拟合的概念
设直线方程为y=a+bx。
则残差为:
,其中。
残差是衡量拟合好坏的重要标志。
可以用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);
holdon
fori=1:
6
plot([i,i],[y(i),y1(i)],'r');
end
可以绘制出如下图形:
三个准则:
(1)最小
(2)最小
(3)最小
3.最小二乘法的直线拟合
问题:
对于给定的数据点,求一次多项式y=a+bx,使得总误差Q最小。
其中。
根据
故有以下方程组(正则方程):
例1.给定数据表,求最小二乘拟合一次多项式
xi
165
123
150
123
141
yi
187
126
172
125
148
解:
N=5,=702,=758,=99864,=108396。
则有方程组
解得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.最小二乘法的多项式拟合
问题:
对于给定的数据点,求m次多项式(m<其中。
根据
故有正则方程:
当m=2时,有
例2.求数据表的最小二乘法拟合的二次多项式函数
xi
-1
-0.75
-0.5
-0.25
0
0.25
0.5
0.75
1
yi
50
40
25
20
18
21
35
56
66
在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.从三次多项式上找出21个点,然后对这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');
holdon
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.解下列超定方程组
解:
设超定方程的解为。
方法一:
点到4条直线的距离平方分别为:
,,,
设,根据,有:
=0
=0
化简有:
解得
方法二:
最小二乘法:
点关于4条直线的残差平方和为:
+
根据,有:
=0
=0
化简有:
解得
用MATLAB命令有:
symsx0y0
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
方法三:
最小二乘法(矩阵运算)
针对方程组的最小二乘近似解即为方程组的解
于是,在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');
holdon
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),求一个多项式函数在区间[a,b]上逼近f(x)。
解决方法:
函数的最佳平方逼近。
令
,使Q最小,则有
例1.求一次多项式在上逼近函数。
解:
构造直线为:
,,则有
,
解得:
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)用最小二乘法求一个形如的经验公式,使它与下列数据表拟合。
xi
19
25
31
38
41
yi
19.0
32.3
49.0
73.3
97.8
解:
方法一,最小二乘法;方法二,用解超定方程组的思路来解题。
,根据,有:
在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)已知数据表如下,试用二次多项式拟合。
xi
0
1
2
3
4
5
6
yi
15
14
14
14
14
15
16
(3)求一个形如(a,b为常数,a>0)的经验公式,使它能和下表数据拟合。
xi
1
1.25
1.5
1.75
2
yi
5.1
5.79
6.53
7.45
8.46
解:
公式可以变为:
lny=lna+bx,进一步可以写为Y=A+bx。
其中Y=lny,A=lna,对应表格为:
xi
1
1.25
1.5
1.75
2
yi
5.1
5.79
6.53
7.45
8.46
Yi
1.629
1.756
1.867
2.008
2.135
(4)求函数在区间[1/4,1]上的最小一次式。