第九章插值与拟合.docx
- 文档编号:14094434
- 上传时间:2023-06-20
- 格式:DOCX
- 页数:30
- 大小:186.91KB
第九章插值与拟合.docx
《第九章插值与拟合.docx》由会员分享,可在线阅读,更多相关《第九章插值与拟合.docx(30页珍藏版)》请在冰点文库上搜索。
第九章插值与拟合
2)恰好给出n-1个方程
(3)
第九章插值与拟合
插值:
求过已知有限个数据点的近似函数。
拟合:
已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义
下它在这些点上的总偏差最小。
插值和拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二
者的数学方法上是完全不同的。
而面对一个实际问题,究竟应该用插值还是拟合,有时容易确定,有时则并不明显。
§1插值方法
下面介绍几种基本的、常用的插值:
拉格朗日多项式插值、牛顿插值、分段线性插
值、Hermite插值和三次样条插值。
1.1拉格朗日多项式插值
1.1.1插值多项式
用多项式作为研究插值的工具,称为代数插值。
其基本问题是:
已知函数f(x)在
区间[a,b]上n1个不同点Xo,X1,…,Xn处的函数值y^f(x)(^0,1/,n),求一个至多n次多项式
:
:
n(x)•a/anXn
(1)
使其在给定点处与f(x)同值,即满足插值条件
叫(G=f(Xi)=yi(i=0,1,…,n)
(2)
n(X)称为插值多项式,Xi(i=0,1/,n)称为插值节点,简称节点,[a,b]称为插值区间。
从几何上看,n次多项式插值就是过n1个点(片,f(xj)(i=0,1,…,n),作一条多项式曲线y=n(x)近似曲线y=f(x)。
n次多项式
(1)有n•1个待定系数,由插值条件(
玄+a1X°+a2X0+…+anX0=y
』a()+aM1+a2X:
+…+anX;=y1
是范德蒙特(Vandermonde)行列式。
当x0,x1/,xn互不相同时,此行列式值不为零。
因此方程组(3)有唯一解。
这表明,只要n1个节点互不相同,满足插值要求
(2)的插值多项式
(1)是唯一的。
插值多项式与被插函数之间的差
Rn(X)=f(X)-n(X)
称为截断误差,又称为插值余项。
当f(x)充分光滑时,
f(n书化)
Rn(X)=f(X)-Ln(X)nl(X),(a,b)
(n+1)!
n
其中r1(X):
|丨(X—Xj)。
y
1.1.2拉格朗日插值多项式
实际上比较方便的作法不是解方程(3)求待定系数,而是先构造一组基函数
|.(x)二(X—Xo)(X—XiJ(X—XiJ(X—Xn)
1(Xi—Xo)…(Xi—XiJ(Xi—XiG…(X—Xn)
h(x)是n次多项式,满足
©li(Xj)二I
nn
Ln(X)='yili(X)='yi
i=0i=0
上式称为n次Lagrange插值多项式,由方程(3)解的唯一性,n1个节点的n次Lagrange插值多项式存在唯一。
1.1.3用Matlab作Lagrange插值
Matlab中没有现成的Lagrange插值函数,必须编写一个M文件实现Lagrange插值。
设n个节点数据以数组x0,y0输入(注意Matlat的数组下标从1开始),m个插值
点以数组x输入,输出数组y为m个插值。
编写一个名为lagrange.m的M文件:
functiony=lagrange(x0,y0,x);
n=length(x0);m=length(x);fori=1:
m
z=x(i);
s=0.0;
fork=1:
np=1.0;
forj=1:
n
ifj〜=kp=p*(z-x0(j))/(x0(k)-x0(j));end
ends=p*y0(k)+s;endy(i)=s;
end
1.2牛顿(Newton)插值
在导出Newton公式前,先介绍公式表示中所需要用到的差商、差分的概念及性质。
1.2.1差商
定义设有函数f(X),X0,X,,X2,…为一系列互不相等的点,称
f(Xi)—f(Xj)
--(i=j)为f(X)关于点Xi,Xj一阶差商(也称均差)记为f[Xi,Xj],即
Xi_Xj
f[Xi,Xj]-f[Xj,Xk]
Xi-Xk
为f(X)关于点Xi,Xj,Xk的二阶差商,记为f[Xi,Xj,Xk]。
一般地,称f[Xo,Xi,,Xk」]-f[X!
X2,,Xk]
Xo-Xk
为f(X)关于点XoX,,Xk的k阶差商,记为
flXo’X!
,Xk」]-f[X1)X2/,Xk]f[Xo,X!
,Xk]:
Xo—Xk
容易证明,差商具有下述性质:
f[Xi,Xj]二f[Xj,Xi]
f[Xi,Xj,Xk]=f[Xi,Xk,Xj]二f[Xj,Xi,Xk]
1.2.2Newton插值公式
线性插值公式可表成
®1(X)=f(Xo)+(X-Xo)f[Xo,X1]
称为一次Newton插值多项式。
一般地,由各阶差商的定义,依次可得
f(x)二f(Xo)(x-Xo)f[X,Xo]
f[X,Xo]=f[Xo*](X-Xj口纳“]
f[X,Xo,X1]二
f[Xo,X1,X2](X-X2)f[X,Xo,X1,X2〕
f[X,Xo,,Xn4Hf[Xo,X1,,Xn](X-Xn)f[X,X°,,Xn]
将以上各式分别乘以1,(x—Xo),(x-Xo)(x-xj,,(X-Xo)(X-Xj(X-Xn",然后相加并消去两边相等的部分,即得
f(x)=f(Xo)+(X—Xo)f[Xo,X1]+…
(x-Xo)(X-Xj(X-XnGf[Xo,X1,,Xn]
(X-X°)(X-Xj(X-Xn)f[X,X°,X1,X]
记
Nn(X)二f(X。
)(X-Xo)f[Xo,Xi]
Xn]
(x-Xo)(X-Xi)(X-Xn」)f[X°,Xi,
Rn(X)=(X—Xo)(X—Xi)(X—Xn)f[X,X°,Xi,,Xn]
二ni(X)f[X,X°,Xi;,Xn]
多项式。
这种形式的插值多项式称为Newton插值多项式。
Rn(x)称为Newton插值余
项。
Newton插值的优点是:
每增加一个节点,插值多项式只增加一项,即
Nni(x)=Nn(X)(X-Xo)(X-Xn)f[Xo,Xi,,X1]
因而便于递推运算。
而且Newton插值的计算量小于Lagrange插值。
由插值多项式的唯一性可知,Newton插值余项与Lagrange余项也是相等的,即
Rn(X)=外1(X)f[X,X°,X1,,Xn]
f(n1)()■
-:
n1(x)(a,b)
(n1)!
由此可得差商与导数的关系
f(n)(®
f[X°,Xi,,Xn]
n!
其中(:
,J,:
二min{xi}「二max{Xj}。
0<<0<<
1.2.3差分
当节点等距时,即相邻两个节点之差(称为步长)为常数,Newton插值公式的形
式会更简单。
此时关于节点间函数的平均变化率(差商)可用函数值之差(差分)来表示。
定义设有等距节点xk=x0+kh(k=0,1,…,n),步长h为常数,fk=f(xk)。
称相邻两个节点xk,xk1处的函数值的增量fk彳-fk(k=0,1,…,n-1)为函数f(x)在点xk处以h为步长的一阶差分,记为■■:
fk,即
%7(k=0,1厂,n)
类似地,定义差分的差分为高阶差分。
如二阶差分为
A2f^Afk+-Afk(k=0,1,…,n-2)
-般地,m阶差分为
左壮十一A""(k=2,3,…),
上面定义的各阶差分又称为向前差分。
常用的差分还有两种:
▽fk=fk-fk4
上式称为Newton向前插值公式。
1.3分段线性插值
1.3.1插值多项式的振荡
用Lagrange插值多项式Ln(x)近似f(x)(a乞x岂b),虽然随着节点个数的增加,Ln(x)的次数n变大,多数情况下误差|尺&)|会变小。
但是n增大时,Ln(x)的光滑性变坏,有时会出现很大的振荡。
理论上,当n—・,在[a,b]内并不能保证Ln(x)处处收敛于f(x)。
Runge给出了一个有名的例子:
1
f(x)2,x[-5,5]
1+x
对于较大的IXI,随着n的增大,J(x)振荡越来越大,事实上可以证明,仅当|x匸3.63
时,才有limLn(x)二f(x),而在此区间外,Ln(x)是发散的。
n^Cnn
高次插值多项式的这些缺陷,促使人们转而寻求简单的低次多项式插值。
1.3.2分段线性插值
简单地说,将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性
函数(i=0,1,…,n)。
In(X)可以表示为
n
ln(X)=Th(X)
izz0
In(x)有良好的收敛性,即对于X[a,b]有,
limIn(x)=f(x)。
n_・
用ln(x)计算x点的插值时,只用到x左右的两个节点,计算量与节点个数n无关。
但n越大,分段越多,插值误差越小。
实际上用函数表作插值计算时,分段线性插值就足够了,如数学、物理中用的特殊函数表,数理统计中用的概率分布表等。
1.3.3用Matlab实现分段线性插值
用Matlab实现分段线性插值不需要编制函数程序,Matlab中有现成的一维插值函
数interpi。
y=interp1(x0,y0,x,'method')
method指定插值的方法,默认为线性插值。
其值可为:
'nearest'最近项插值
'linear'线性插值
'spline'立方样条插值
'cubic'立方插值。
所有的插值方法要求x0是单调的。
当x0为等距时可以用快速插值法,使用快速插值法的格式为'*nearest'、’*linear'、
'*spline'、'*cubic'。
1.4埃尔米特(Hermite)插值
1.4.1Hermite插值多项式
如果对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有相同的一
阶、二阶甚至更高阶的导数值,这就是Hermite插值问题。
本节主要讨论在节点处插值
函数与函数的值及一阶导数值均相等的Hermite插值。
设已知函数y=f(x)在n1个互异节点x0,x1,…,xn上的函数值y^f(xi)(i=0,1/,n)和导数值y:
=f'(xj(i=0,1,…,n),要求一个至多2n1次的多项式H(x),使得
H(x)=yiH'(x)=y'i(i=0,1,…,n)
满足上述条件的多项式H(X)称为Hermite插值多项式。
Hermite插值多项式为
Matlab中没有现成的Hermite插值函数,必须编写一个M文件实现插值。
设n个节点的数据以数组x0(已知点的横坐标),y0(函数值),y1(导数值)输入(注意Matlat的数组下标从1开始),m个插值点以数组x输入,输出数组y为m个插值。
编写一个名为hermite.m的M文件:
functiony=hermite(x0,y0,y1,x);
n=length(x0);m=length(x);fork=1:
m
yy=0.0;
fori=1:
n
h=1.0;
a=0.0;
forj=1:
n
ifj~=i
h=h*((x(k)-x0(j))/(x0(i)-x0(j)))A2;
a=1/(x0(i)-x0(j))+a;
endend
yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));
end
y(k)=yy;
end
1.5样条插值
许多工程技术中提出的计算问题对插值函数的光滑性有较高要求,如飞机的机翼外
形,内燃机的进、排气门的凸轮曲线,都要求曲线具有较高的光滑程度,不仅要连续,而且要有连续的曲率,这就导致了样条插值的产生。
1.5.1样条函数的概念
所谓样条(Spline)本来是工程设计中使用的一种绘图工具,它是富有弹性的细木条或细金属条。
绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线),并
使连接点处有连续的曲率。
三次样条插值就是由此抽象出来的。
数学上将具有一定光滑性的分段多项式称为样条函数。
具体地说,给定区间[a,b]
的一个分划
a=x0:
:
:
为:
:
:
…:
:
:
xn4:
:
:
xn=b
如果函数S(x)满足:
(i)在每个小区间[xi,xi』(i=0,1,,n-1)上S(x)是m次多项式;
(ii)S(x)在[a,b]上具有m_1阶连续导数。
则称S(x)为关于分划的m次样条函数,其图形为m次样条曲线。
显然,折线是一次样条曲线。
1.5.2三次样条插值
禾U用样条函数进行插值,即取插值函数为样条函数,称为样条插值。
例如分段线性插值是一次样条插值。
我们只介绍三次样条插值,即已知函数y二f(x)在区间[a,b]
上的n1个节点
a=Xo:
:
:
Xi:
:
:
…:
:
:
Xn」:
:
:
Xn二b
上的值y=f(xj(i=0,1,,n),求插值函数S(x),使得
(i)S(xJ二y(i=0,1,,n);(5)
(ii)在每个小区间[Xj,Xji](j=0,1,…,n-1)上S(x)是三次多项式,记为Sj(x);
(iii)S(x)在[a,b]上二阶连续可微。
函数S(x)称为f(x)的三次样条插值函数。
由条件(ii),不妨将S(x)记为
S(x)二{Sj(x),x[Xj,Xj1],j=0,1,n_1
Sj(x)二ajx3bjx2cjxdj
其中aj,bj,Cj,dj为待定系数,共4n个。
由条件(iii)
Sj(xj十)=Sj十(xj十)
S''j(XjHt)=S'j十(Xj41)
容易看出,(5)、(6)式共含有4n-2个方程,为确定S(x)的4n个待定参数,尚需再给出2个条件。
常用的三次样条函数的边界条件有3种类型:
(i)S(a)=y'°,S(b)二y'n。
由这种边界条件建立的样条插值函数称为f(x)的完备三次样条插值函数。
特别地,y;=y'n二。
时,样条曲线在端点处呈水平状态。
如果f'(X)不知道,我们可以要求S'(X)与f'(X)在端点处近似相等。
这时以
x0,x1,x2,x3为节点作一个三次Newton插值多项式Na(x),以xn,xn^,xn/,焉」作一个三次Newton插值多项式Nb(x),要求
S(a)=N;(a),S(b)=N'b(b).
由这种边界条件建立的三次样条称为f(x)的Lagrange三次样条插值函数。
(ii)S'(a)二y''°,S'(b)=y''n。
特别地y''^y''^0时,称为自然边界条件。
(iii)S'(aP)=S'(b-0),S''(a•0)=S''(b-0),此条件称为周期条件。
1.5.3三次样条插值在Matlab中的实现
在Matlab中数据点称之为断点。
如果三次样条插值没有边界条件,最常用的方法,
就是采用非扭结(not-a-knot)条件。
这个条件强迫第1个和第2个三次多项式的三阶导数相等。
对最后一个和倒数第2个三次多项式也做同样地处理。
Matlab中三次样条插值也有现成的函数:
y=interp1(x0,y0,x,'spline');
y=spline(x0,y0,x);
pp=csape(xO,yO,conds),
pp=csape(xO,yO,conds,valconds),y=ppval(pp,x)。
其中x0,y0是已知数据点,x是插值点,y是插值点的函数值。
对于三次样条插值,我们提倡使用函数csape,csape的返回值是pp形式,要求插
值点的函数值,必须调用函数ppvalo
pp=csape(x0,y0):
使用默认的边界条件,即Lagrange边界条件。
pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件,其值可为:
'complete'边界为一阶导数,一阶导数的值在valconds参数中给出,若忽略
valconds参数,则按缺省情况处理。
'not-a-knot'非扭结条件
'periodic'周期条件
'second'边界为二阶导数,二阶导数的值在valconds参数中给出,若忽略
valconds参数,二阶导数的缺省值为[0,0]。
'variational'设置边界的二阶导数值为[0,0]。
对于一些特殊的边界条件,可以通过conds的一个12矩阵来表示,conds元素的
取值为0,1,2o
conds(i)=j的含义是给定端点i的j阶导数,即conds的第一个元素表示左边界的条件,第二个元素表示右边界的条件,conds=[2,1]表示左边界是二阶导数,右边界是一阶
导数,对应的值由valconds给出。
详细情况请使用帮助helpcsapeo
例1机床加工
待加工零件的外形根据工艺要求由一组数据(x,y)给出(在平面情况下),用程控
铳床加工时每一刀只能沿x方向和y方向走非常小的一步,这就需要从已知数据得到加
工所要求的步长很小的(x,y)坐标。
表中给出的x,y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变0.1
时的y坐标。
试完成加工所需数据,画出曲线,并求出x=0处的曲线斜率和13乞x空15范围内y的最小值。
""x03__57911121314_15
y01.21.72.02.12.01.81.21.01.6
要求用Lagrange、分段线性和三次样条三种插值方法计算。
解编写以下程序:
x0=[035791112131415];
y0=[01.21.72.02.12.01.81.21.01.6];
x=0:
0.1:
15;
y1=lagrange(x0,y0,x);%前面编写的Lagrange插值函数
y2=interp1(x0,y0,x);
y3=interp1(x0,y0,x,'spline');
pp1=csape(x0,y0);
y4=ppval(pp1,x);
pp2=csape(x0,y0,'second');
y5=ppval(pp2,x);
[X',y1',y2',y3',y4',y5']
subplot(2,2,1)
plot(x0,y0,'+',x,y1)
title('Lagrange')
subplot(2,2,2)plot(xO,yO,'+',x,y2)
title('Piecewiselinear')
subplot(2,2,3)
plot(x0,y0,'+',x,y3)title('Spline1')
subplot(2,2,4)
plot(xO,yO,'+',x,y4)
title('Spline2')
dx=diff(x);
dy=diff(y3);dy_dx=dy./dx;
dy_dxO=dy_dx⑴ytemp=y3(131:
151);ymin=min(ytemp);
index=find(y3==ymin);xmin=x(index);
[xmin,ymin]
计算结果略。
可以看出,拉格朗日插值的结果根本不能应用,分段线性插值的光滑性较差(特别
是在x=14附近弯曲处),建议选用三次样条插值的结果。
1.6二维插值
前面讲述的都是一维插值,即节点为一维变量,插值函数是一元函数(曲线)。
若
节点是二维的,插值函数就是二元函数,即曲面。
如在某区域测量了若干点(节点)的高程(节点值),为了画出较精确的等高线图,就要先插入更多的点(插值点),计算这
些点的高程(插值)。
1.6.1插值节点为网格节点
已知mxn个节点:
仅,yj,Zj)(i=12…,m;j=1,2,…,n),且x^L : : yn。 求点(x,y)处的插值z。 Matlab中有一些计算二维插值的程序。 如 z=interp2(x0,y0,z0,x,y,'method') 其中xO,yO分别为m维和n维向量,表示节点,z0为nm维矩阵,表示节点值,x,y为一维数组,表示插值点,x与y应是方向不同的向量,即一个是行向量,另一个是列 向量,z为矩阵,它的行数为y的维数,列数为x的维数,表示得到的插值,'method' 的用法同上面的一维插值。 如果是三次样条插值,可以使用命令 pp=csape({xO,yO},zO,conds,valconds),z=fnval(pp,{x,y}) 其中xo,yo分别为m维和n维向量,z0为mn维矩阵,z为矩阵,它的行数为x的维数,列数为y的维数,表示得到的插值,具体使用方法同一维插值。 例2在一丘陵地带测量高程,x和y方向每隔100米测一个点,得高程如下表,试 插值一曲面,确定合适的模型,并由此找出最高点和该点的高程。 100200300400500 100 )36 697 624 478 450 200 *98 712 630 478 420 300 680 674 598 412 400 400662626552334310 解编写程序如下: clear,clc x=100: 100: 500; y=100: 100: 400; z=[636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310]; pp=csape({x,y},z') xi=100: 10: 500;yi=100: 10: 400 cz1=fnval(pp,{xi,yi}) cz2=interp2(x,y,z,xi,yi','spline') [i,j]=find(cz1==max(max(cz1))) x=xi(i),y=yi(j),zmax=cz1(i,j) 1.6.2插值节点为散乱节点 已知n个节点: (洛,%,乙)01,2,,n),求点(x,y)处的插值z。 对上述问题,Matlab中提供了插值函数griddata,其格式为: ZI=GRIDDATA(X,Y,Z,XI,YI) 其中X、Y、Z均为n维向量,指明所给数据点的横坐标、纵坐标和竖坐标。 向量XI、 YI是给定的网格点的横坐标和纵坐标,返回值ZI为网格(XI,YI)处的函数值。 XI 与YI应是方向不同的向量,即一个是行向量,另一个是列向量。 例3在某海域测得一些点(x
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第九 章插值 拟合