插值与数据拟合模型
(完整版)Matlab学习系列13.数据插值与拟合

13. 数据插值与拟合实际中,通常需要处理实验或测量得到的离散数据(点)。
插值与拟合方法就是要通过离散数据去确定一个近似函数(曲线或曲面),使其与已知数据有较高的拟合精度。
1.如果要求近似函数经过所已知的所有数据点,此时称为插值问题(不需要函数表达式)。
2.如果不要求近似函数经过所有数据点,而是要求它能较好地反映数据变化规律,称为数据拟合(必须有函数表达式)。
插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数。
区别是:【插值】不一定得到近似函数的表达形式,仅通过插值方法找到未知点对应的值。
【拟合】要求得到一个具体的近似函数的表达式。
因此,当数据量不够,但已知已有数据可信,需要补充数据,此时用【插值】。
当数据基本够用,需要寻找因果变量之间的数量关系(推断出表达式),进而对未知的情形作预测,此时用【拟合】。
一、数据插值根据选用不同类型的插值函数,逼近的效果就不同,一般有:(1)拉格朗日插值(lagrange插值)(2)分段线性插值(3)Hermite(4)三次样条插值Matlab 插值函数实现:(1)interp1( ) 一维插值(2)intep2( ) 二维插值(3)interp3( ) 三维插值(4)intern( ) n维插值1.一维插值(自变量是1维数据)语法:yi = interp1(x0, y0, xi, ‘method’)其中,x0, y0为原离散数据(x0为自变量,y0为因变量);xi为需要插值的节点,method为插值方法。
注:(1)要求x0是单调的,xi不超过x0的范围;(2)插值方法有‘nearest’——最邻近插值;‘linear’——线性插值;‘spline’——三次样条插值;‘cubic’——三次插值;默认为分段线性插值。
例1 从1点12点的11小时内,每隔1小时测量一次温度,测得的温度的数值依次为:5,8,9,15,25,29,31,30,22,25,27,24.试估计每隔1/10小时的温度值。
matlab插值与拟合

matlab插值与拟合
在MATLAB中,插值和拟合都是通过函数来实现的。
插值是通过创建新的数据点来填充在已知数据点之间的空白。
MATLAB提供了几种不同的插值方法,例如分段线性插值、三次样条插值、立方插值等。
具体使用哪种插值方法取决于数据的特性和所需的精度。
插值函数的一般形式是`interp1(x, y, xi, 'method')`,其中`x`和`y`是已知的数据点,`xi`是待插值点的横坐标向量,`method`是插值方法,例如最近邻点插值、线性插值、三次样条插值、立方插值等。
拟合是通过调整一个数学模型来使得该模型尽可能地接近给定的数据点。
在MATLAB中,可以使用`polyfit`函数进行多项式拟合。
该函数的一般形式是`p = polyfit(x, y, n)`,其中`x`和`y`是已知的数据点,`n`是多项式的阶数。
该函数返回一个向量`p`,表示多项式的系数。
可以使用`polyval`函数来评估这个多项式模型在给定数据点上的值。
需要注意的是,插值和拟合都是数学上的近似方法,它们只能尽可能地逼近真实的情况,而不能完全准确地描述数据的变化。
因此,选择合适的插值和拟合方法是非常重要的。
数值计算插值法与拟合实验

dy0=-10.*(1-5.^4)./(1+5.^4).^2;dyn=10.*(1-5.^4)./(1+5.^4).^2;
m=maspline(x1,y1,dy0,dyn,xx);
plot(xx,m,'ok')
2、
程序:
x=[-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5]';
plot(xx,m,'ok')
第二个方程
程序
x=-5:0.2:5;
y=atan(x);
plot(x,y,'r');
hold on
x1=-5:1:5;
y1=atan(x1);
xx=-4.5:0.5:4.5;
yy=malagr(x1,y1,xx);
plot(xx,yy,'+')
dy0=1./(1+25);dyn=1./(1+25);
实验报告三
一、实验目的
通过本实验的学习,各种插值法的效果,如多项式插值法,牛顿插值法,样条插值法,最小二乘法拟合(即拟合插值),了解它们各自的优缺点及插值。
二、实验题目
1、插值效果比较
实验题目:将区间 10等份,对下列函数分别计算插值节点 的值,进行不同类型的插值,作出插值函数的图形并与 的图形进行比较:
y=[-4.45 -0.45 0.55 0.05 -0.44 0.54 4.55]';
plot(x,y,'or');hold on
%三.2:1.5;
y1=p1(1)*x1.^3+p1(2)*x1.^2+p1(3)*x1+p1(4);
数学建模之插值与拟合

matlab中拟合的函数
非线性曲线拟合 Matlab中对于多项式拟合,有现成的函数
c = lsqcurvefit ( ′fun′, x0, xdata, ydata)
matlab中拟合的函数
非线性曲线拟合例题
对下面的x、y进行数据拟合
x=[3.6,7.7,9.3,4.1,8.6,2.8,1.3,7.9,10,5.4]; y=[16.5,150.6,263.1,24.7,208.5,9.9,2.7,163.9,325,54.3];
最小二乘法
线性最小二乘法是解决曲线拟合最常用的方法,基本 思路是,令
f (x) a1r1(x) a2r2 (x) amrm (x) • 其中,rk(x)是事先选定的一组线性无关的函数,ak是待定系
数(k=1,2,...,m,m<n)。拟a合准则是使yi,i=1,2,3...,n,与f (xi )
• 求:利用最小二乘法求得上述拟合函数
求解方法
(1)做散点图,通过散点图判断函数为:y=ax+b
(2)根据最小二乘法原理可知,即使下式中M最小
10
M yi axi b2
i 1
(3)把M看作是自变量为a和b的函数,由多元函数取最值
的条件可知:
M M
a b
a, a,
b b
0 0
M
a
M
b
目录
1
插值法与拟合法
2 matlab中插值的函数
3 matlab中拟合的函数 4 插值与拟合的运用
插值法与拟合法的基本介绍
插值法:求过已知有限个数据点的近似函数。
拟合法:已知有限个数据点,求近似函数,不要求
过已知数据点,只要求在某种意义下它在这些点上 的总偏差最小。
插值与拟合方法

插值与拟合方法在实际中,常常要处理由实验或测量所得到的一批离散数据.插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻找某个近似函数,使所得到的近似函数与已知数据有较高的拟合精度.插值问题:要求这个近似函数(曲线或曲面)经过所已知的所有数据点.通常插值方法一般用于数据较少的情况.数据拟合:不要求近似函数通过所有数据点,而是要求它能较好地反映数据的整体变化趋势。
共同点:插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法,由于对近似要求的准则不同,因此二者在数学方法上有很大的差异.插值问题的一般提法:已知某函数)(x f y =(未知)的一组观测(或试验)数据),,2,1)(,(n i y x ii⋅⋅⋅=,要寻求一个函数)(x φ,使iiy x =)(φ),,2,1(n i ⋅⋅⋅=,则)()(x f x ≈φ.实际中,常常在不知道函数)(x f y =的具体表达式的情况下,对于i x x =有实验测量值iy y =),,2,1,0(n i ⋅⋅⋅=,寻求另一函数)(x φ使满足:)()(i i i x f y x ==φ),,2,1,0(n i ⋅⋅⋅=称此问题为插值问题,并称函数)(x φ为)(x f 的插值函数,nx x x x ,,,,21⋅⋅⋅称为插值节点,),,2,1,0()(n i y x ii⋅⋅⋅==φ称为插值条件,即)()(iiix f y x ==φ),,2,1,0(n i ⋅⋅⋅=,则)()(x f x ≈φ.(1) 拉格朗日(Lagrange )插值设函数)(x f y =在1+n 个相异点nx x x x ,,,,21⋅⋅⋅上的函数值为ny y y y ,,,,21⋅⋅⋅,要求一个次数不超过n 的代数多项式nnnx a x a x a a x P +⋅⋅⋅+++=221)(使在节点i x 上有),,2,1,0()(n i y x P ii n ⋅⋅⋅==成立,称之为n 次代数插值问题,)(x P n称为插值多项式.可以证明n 次代数插值是唯一的.事实上: 可以得到j n j n i i j in y x x xx x P j i ∑∏==⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫⎝⎛--=≠00)()( 当1=n 时,有二点一次(线性)插值多项式:101001011)(y x x x x y x x x x x P --+--=当n =2时,有三点二次(抛物线)插值多项式:2120210121012002010212))(())(())(())(())(())(()(y x x x x x x x x y x x x x x x x x y x x x x x x x x x P ----+----+----=(2)牛顿(Newton ) 插值牛顿插值的基本思想:由于)(x f y =关于二节点10,x x 的线性插值为)()()()()()()()()(00101000010101x x x x x f x f x p x x x x x f x f x f x p ---+=---+= 假设满足插值条件)2,1,0()()(2===i x p y x f iii的二次插值多项式一般形式为))(()()(1212x x x x c x x c c x p --+-+= 由插值条件可得⎪⎩⎪⎨⎧=--+-+=-+=)())(()()()()(21202202101011000x f x x x x c x x c c x f x x c c x f c 可以解出⎪⎪⎪⎩⎪⎪⎪⎨⎧------=--==020101121220101100)()()()()()(),(x x x x x f x f x x x f x f c x x x f x f c x f c所以))(()())(()()(10211020102x x x x c x p x x x x c x x c c x p --+=--+-+=类似的方法,可以得到三次插值多项式等,按这种思想可以得到一般的牛顿插值公式.函数的差商及其性质对于给定的函数)(x f ,用),,,(10n x x x f ⋅⋅⋅表示关于节点nx x x ,,,1⋅⋅⋅的n 阶差商,则有一阶差商:01011)()(),(x x x f x f x x f --=,121221)()(),(x x x f x f x x f --= 二阶差商:021021210),(),(),,(x x x x f x x f xx x f --=n 阶差商:0110211),,,(),,,(),,,(x x x x x f x x x f x x x f n n n n -⋅⋅⋅-⋅⋅⋅=⋅⋅⋅-差商有下列性质:(1)差商的分加性:∑∏=≠=-=⋅⋅⋅nk nk j j j kk n x xx f xx x f 0)(01)()(),,,(.(2)差商的对称性:在),,,(1nx x x f ⋅⋅⋅中任意调换jix x ,的次序其值不变.牛顿插值公式: 一次插值公式为))(,()()(01001x x x x f x f x p -+=二次插值公式为))()(,,()())()(,,())(,()()(1021011021001002x x x x x x x f x p x x x x x x x f x x x x f x f x p --+=--+-+=于是有一般的牛顿插值公式为)())()(,,,()()())()(,,,())()(,,())(,()()(11010111010102100100----⋅⋅⋅--⋅⋅⋅+=-⋅⋅⋅--⋅⋅⋅+⋅⋅⋅+--+-+=n n n n n n x x x x x x x x x f x p x x x x x x x x x f x x x x x x x f x x x x f x f x p可以证明:其余项为))(())()(,,,,()(11010n n n x x x x x x x x x x x x f x R --⋅⋅⋅--⋅⋅⋅=-实际上,牛顿插值公式是拉格朗日插值公式的一种变形,二者是等价的.另外还有著名的埃尔米特(Hermite )插值等.(3)样条函数插值方法样条,实质上就是由分段多项式光滑连接而成的函数,一般称为多项式样条.由于样条函数的特殊性质,决定了样条函数在实际中有着重要的应用.样条函数的一般概念定义 设给定区间],[b a 的一个分划b x x x a n=<⋅⋅⋅<<=∆1:,如果函数)(x s 满足条件:(1) 在每个子区间),,2,1](,[1n i x x ii ⋅⋅⋅=-上是k 次多项式; (2) )(x s 及直到k -1阶的导数在],[b a 上连续.则称)(x s 是关于分划△的一个k 次多项式样条函数,nx x x ,,,1⋅⋅⋅称为样条节点,121,,,-⋅⋅⋅n x x x 称为内节点,nx x ,0称为边界节点,这类样条函数的全体记作),(k S P∆,称为k 次样条函数空间.若),()(k S x s P∆∈,则)(x s 是关于分划△的k 次多项式样条函数.k 次多项式样条函数的一般形式为∑∑=-=+-+=ki n j k j jii k x x k i x x s 011)(!!)(βα其中),,1,0(k i i=α和)1,,2,1(-=n j jβ均为任意常数,而)1,,2,1(,0,)()(-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x jj kj kj在实际中最常用的是2=k 和3的情况,即为二次样条函数和三次样条函数. 二次样条函数:对于],[b a 上的分划b x x x a n=<⋅⋅⋅<<=∆1:,则)2,()(!2!2)(11222102∆βαααP n j j jS x x x x x s ∈-+++=∑-=+其中)1,2,1(,0,)()(22-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x j j j j . 三次样条函数:对于],[b a 上的分划b x x xa n =<⋅⋅⋅<<=∆10:,则)3,()(!3!3!2)(1133322103∆βααααP n j j jS x x x x x x s ∈-++++=∑-=+其中)1,2,1(,0,)()(33-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x jjj j .1 二次样条函数插值)2,()(2∆∈P S x s 中含有2+n 个待定常数,故应需要2+n 个插值条件,因此,二次样条插值问题可分为两类:问题(1):已知插值节点ix 和相应的函数值),,2,1,0(n i y i⋅⋅⋅=,以及端点0x (或n x )处的导数值0'y (或ny '),求)2,()(2∆∈PS x s 使得⎩⎨⎧'=''='⋅⋅⋅==))(()(),,2,1,0()(20022n n i i y x s y x s n i y x s 或(5.1)问题(2):已知插值节点ix 和相应的导数值),,2,1,0(n i y i⋅⋅⋅=',以及端点0x (或n x )处的函数值0y (或ny ),求)2,()(2∆∈P S x s 使得⎩⎨⎧==⋅⋅⋅='='))(()(),,2,1,0()(20022n n i i y x s y x s n i y x s 或(5.2)事实上,可以证明这两类插值问题都是唯一可解的.对于问题(1),由条件(5.1)⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=+='==-+++==++==++=∑-=00210211222102121211112020201002)(,,3,2,)(2121)(21)(21)(y x x s n j y x x x x x s yx x x s y x x x s j j i i j i jj j ααβααααααααα 引入记号T n ),,,,,(11210-=ββααα X 为未知向量,T nn y y y y ),,,,(10'= C 为已知向量, ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-0010)(21)(21211)(212110211211021212212222211200x x x x x x x x x x x x x x x n n n n n A 于是,问题转化为求方程组C AX =的解Tn ),,,,,(1121-=ββααα X 的问题,即可得到二次样条函数的)(2x s 的表达式.对于问题(2)的情况类似.2.三次样条函数插值由于)3,()(3∆∈P S x s 中含有3+n 个待定系数,故应需要3+n 个插值条件,因此可将三次样条插值问题分为三类: 问题(1):已知插值节点jx 和相应的函数值),,2,1,0(n j y j⋅⋅⋅=,以及两个端点0x ,n x 处的导数值0'y ,ny ',求)3,()(3∆∈PS x s 使满足条件⎪⎩⎪⎨⎧='='⋅⋅⋅==),0()(),,1,0()(33n j y x s n j y x s j j j j(5.3)问题(2):已知插值节点jx 和相应的函数值),,2,1,0(n j y j⋅⋅⋅=,以及两个端点0x ,nx 处的二阶导数值0y '',n y '',求)3,()(3∆∈PS x s 使满足条件⎪⎩⎪⎨⎧=''=''⋅⋅⋅==),0()(),,1,0()(33n j y x s n j y x s j j j j(5.4)问题(3):类似地,求)3,()(3∆∈PSx s 使满足条件⎪⎩⎪⎨⎧=+=-==)2,1,0)(0()0(),,1,0()(0)(3)(33k x s x s n j y x s k n k j j(5.5)这三类插值问题的条件都是3+n 个,可以证明其解都是唯一的〔8〕.一般的求解方法可以仿照二次样条的情况处理方法,在这里给出一种更简单的方法.仅依问题(1)为例,问题(2)和问题(3)的情况类似处理.由于在)3,()(3∆PS x s ∈区间],[b a 上是一个分段光滑,且具有二阶连续导数的三次多项式,则在子区间],[1+j jx x 上)(3x s ''是线性函数,记),,,1,0)((3n j x s d jj =''=为待定常数.由拉格朗日插值公式可得nj x x h h x x d h x x d x s j j j jj j jj j ,,1,0,,)(1113=-=-+-=''+++显然jjj h d dx s -='''+13)(在],[1+j jx x上为常数.于是在],[1+j j x x 上有31233)(6)(2))(()(j jjj j j j j j x x h d d x x d x x x s y x s --+-+-'+=+(5.6)则当1+=j x x 时,由(5.6)式和问题(1)的条件得121231362)()(+++=-++'+=j j jj j j j j j j y h d d h d h x s y x s故可解得)2(6)(113+++--='j j j jjj j d d h h y y x s(5.7)将(5.7)式代入(5.6)式得)1,,1,0](,[,)(6)(2)()2(6)(1312113-=∈--+-+-⎥⎥⎦⎤⎢⎢⎣⎡+--+=++++n j x x x x x h d d x x d x x d d h h y y y x s j j j jj j j jj j j j j j j j(5.8) 在],[1j j x x-上同样的有),,2,1](,[,)(6)(2)()2(6)(131112111111113n j x x x x x h d d x x d x x d d h h y y y x s j j j j j j j j j j j j j j j j =∈--+-+-⎥⎥⎦⎤⎢⎢⎣⎡+--+=------------(5.9) 根据)(3x s的一阶导数连续性,由(5.9)式得)()2(6)0(311113j j j j j j j j x s d d h h y y x s '=++-=-'---- 结合(5.7)式整理得⎪⎪⎭⎫ ⎝⎛---+=++++--+-+----11111111162j j j j j j j j j j j j j j j j j h y y h y y h h d h h h d d h h h 引入记号⎪⎪⎭⎫ ⎝⎛---+=+=--+--111116,j j j j j j j j j j j j j h y y h y y h h c h h h a ,111--+=-j j j j h h h a .则)1,,2,1(,2)1(11-==++-+-n j c d a d d a j j j j j j(5.10)再由边界条件:nny x s y x s '=''=')(,)(33得⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛--'=+⎪⎪⎭⎫ ⎝⎛'--=+----111100010106262n n n n n n n h y y y h d d y h y y h d d(5.11)联立(5.10),(5.11)式得方程组C D A =⋅(5.12)其中⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=----2121212112112200n n n n a a a a a aA ,⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=-n n d d d d 110 D ,⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--'⎪⎪⎭⎫ ⎝⎛'--=----111110001066n n n n n n hy y y h c c y h y y h C 由方程组(6.12)可以唯一解出),,1,0(n j d j=,代入(5.8)式就可以得三次样条函数)(3x s 的表达式.B样条函数插值方法磨光函数实际中的许多问题,往往是既要求近似函数(曲线或曲面)有足够的光滑性,又要求与实际函数有相同的凹凸性,一般插值函数和样条函数都不具有这种性质.如果对于一个特殊函数进行磨光处理生成磨光函数(多项式),则用磨光函数构造出样条函数作为插值函数,既有足够的光滑性,而且也具有较好的保凹凸性,因此磨光函数在一维插值(曲线)和二维插值(曲面)问题中有着广泛的应用.由积分理论可知,对于可积函数通过积分会提高函数的光滑度,因此,我们可以利用积分方法对函数进行磨光处理.定义 若)(x f 为可积函数,对于0>h ,则称积分⎰+-=22,1)(1)(hx h x h dt t f h x f为)(x f 的一次磨光函数,h 称为磨光宽度.同样的,可以定义)(x f 的k 次磨光函数为)1()(1)(22,1,>=⎰+--k dt t f h x f hx h x h k h k事实上,磨光函数)(,x f h k 比)(x f 的光滑程度要高,且当磨光宽度h 很小时)(,x f h k 很接近于)(x f .等距B样条函数对于任意的函数)(x f ,定义其步长为1的中心差分算子δ如下:⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+=2121)(x f x f x f δ在此取0)(+=x x f ,则002121+++⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+=x x x δ是一个单位方波函数(如图5-1),记0)(+=Ωx x δ.并取1=h ,对)(0x Ω进行一次磨光得++++-+++-+++--+-+=-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+==⎰⎰⎰⎰)1(2)1(2121)()(11212100212101x x x dt t dt t dt t t dt t x x xx x x x x x ΩΩ显然)(1x Ω是连续的(如图5-2).)(1x Ωo1-1/2 0 1/2 x -1 0 1 x 图5-1图5-2类似地可得到k 次磨光函数为kk j jk j k j k x k C x ++=+⎪⎭⎫ ⎝⎛-++-=Ω∑21!)1()(11 实际上,可以证明:)(x kΩ是分段k 次多项式,且具有1-k 阶连续导数,其k 阶导数有2+k个间断点,记为)1,,2,1,0(21+⋅⋅⋅=+-=k j k j x j.从而可知)(x kΩ是对应于分划+∞<<⋅⋅⋅<<<-∞∆+110:k x x x 的k 次多项式样条函数,称之为基本样条函数,简称为k 次B样条.由于样条节点为)1,,2,1,0(21+⋅⋅⋅=+-=k j k j xj是等距的,故)(x k Ω又称为k 次等距B样条函数.对于任意函数)(x f 的k 次磨光函数,由归纳法可以得到 [4,8] :⎪⎭⎫⎝⎛+≤≤--Ω=⎰∞+∞--22)()(1)(1,h x t h x dt t f htx h x f k h k 特别地,当1)(=x f 时,有1)(11⎰+∞∞--=-dt htx hk Ω,从而1)(⎰+∞∞-=dx x k Ω,且当k ≥1时有递推关系⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛-Ω⎪⎭⎫ ⎝⎛---⎪⎭⎫ ⎝⎛+Ω⎪⎭⎫ ⎝⎛++=Ω--212121211)(11x x k x k x k x k k k一维等距B样条函数插值等距B样条函数与通常的样条如下的关系: 定理设有区间],[b a 的均匀分划nab h n j jh x x j -=⋅⋅⋅=+=),,,1,0(:0∆,则对任意 k 次样条函数),()(k S x S p k ∆∈都可以表示为B样条函数族1021-=-=⎭⎬⎫⎩⎨⎧⎪⎭⎫⎝⎛+---n j k j k k j h x x Ω的线性组合[14].根据定理 5.1,如果已知曲线上一组点()jjy x ,,其中),,1,0,0(0n j h jh x x j⋅⋅⋅=>+=,则可以构造出一条样条磨光曲线(即为B样条函数族的线性组合)⎪⎭⎫⎝⎛--=∑--=j h x x c x S n kj k j k 01)(Ω 其中)1,,1,(-⋅⋅⋅+--=n k k j c j为待定常数.用它来逼近曲线,既有较好的精度,又有良好的保凸性.实际中,最常用的是3=k 的情况,即一般形式为⎪⎭⎫ ⎝⎛--=∑+-=j h x x c x S n j j 01133)(Ω 其中3+n 个待定系数)1,,0,1(+⋅⋅⋅-=n j c j可以由三类插值条件确定.由插值条件(5.3)得()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=-'='==-='=-'='∑∑∑+-=+-=+-=n n j j n i n j j i n j j y j n c h x S ni y j i c x S y j c h x S 113311330113031)(,,1,0,)(1)(ΩΩΩ(5.13)注意到)(3x Ω的局部非零性及其函数值:61)1(,32)0(33=±=ΩΩ,当2≥x 时0)(3=x Ω;且由)21()21()(223--+='x x x ΩΩΩ知,21)1(,0)0(33=±'='ΩΩ,当2≥x 时0)(3='x Ω.则(5.13)中的每一个方程中只有三个非零系数,具体的为⎪⎩⎪⎨⎧'=+-==++'=+-+-+--n n n i i i i y h c c n i y c c c y h c c 2,,1,0,6421111011(5.14)由方程组(5.14)容易求解出)1,,0,1(+⋅⋅⋅-=n j c j,即可得到三次样条函数)(3x S 表达式.类似地,由插值条件(5.4)得待定系数的)1,,0,1(+⋅⋅⋅-=n j c j所满足的方程组为⎪⎩⎪⎨⎧''=+-==++''=+-+-+--nn n n i i i i y h c c c n i y c c c y h c c c 21111021012,,1,0,642(5.15)由插值条件(5.5)得待定系数的)1,,0,1(+⋅⋅⋅-=n j cj所满足的方程组为⎪⎪⎩⎪⎪⎨⎧==++=-+---=-++-=-+-+-+-+--+--+--ni y c c c c c c c c c c c c c c c c c c c i i i i n n n n n n n n ,,1,0,640)()(2)(0)(0)(0)()(4)(1111011111111011(5.16)方程组(5.15),(5.16)也都是容易求解的.注:上述等距B样条插值公式也适用于近似等距的情形,但在端点0x 和n x 处误差可能较大,实际应用时,为了提高在端点0x 和nx 处的精度,可以适当向左右延拓几个节点.二维等距B样条函数插值设有空间曲面),(y x f z =(未知),如果已知二维等距节点()()τj y ih x y x ji++=0,,)0,(>τh 上的值为),,2,1,0;,,2,1,0(m j n i z ij⋅⋅⋅=⋅⋅⋅=,则相应的B样条磨光曲面的一般形式为⎪⎭⎫ ⎝⎛--⎪⎭⎫⎝⎛--=∑∑--=--=j y y i h x x c y x s l m lj k ij n ki τΩΩ0011),( 其中),,2,1,0;,,2,1,0(m j n i c ij⋅⋅⋅=⋅⋅⋅=为待定常数,l k ,可以取不同值,常用的也是2,=l k 和3的情形.这是一种具有良好保凸性的光滑曲面(函数),在工程设计中是常用的,但只能使用于均匀分划或近似均匀分划的情况.(4) 最小二乘拟合方法最小二乘拟合方法的思想:由于一般插值问题并不总是可解的(即当插值条件多于待定系数的个数时,其问题无解),同时,问题的插值条件本身一般是近似的,为此,只要求在节点上近似地满足插值条件,并使它们的整体误差最小,这就是最小二乘拟合法.最小二乘拟合方法可以分为线性最小二乘拟合方法和非线性最小二乘拟合方法.线性最小二乘拟合方法设{}m k kx 0)(=φ是一个线性无关的函数系,则称线性组合∑==mk k k x a x 0)()(φφ为广义多项式.如三角多项式:∑∑==+=mk k mk kkx b kx ax 0sin cos )(φ.设由给定的一组测量数据),(iiy x 和一组正数),,2,1(n i w i⋅⋅⋅=,求一个广义多项式∑==mk k k x a x 0)()(φφ使得目标函数[]21)(∑=-=ni i i i y x w S φ(5.17)达到最小,则称函数)(x φ为数据),,2,1)(,(n i y x ii⋅⋅⋅=关于权系数),,2,1(n i w i⋅⋅⋅=的最小二乘拟合函数,由于)(x φ关于待定系数ia 是线性的,故此问题又称为线性最小二乘问题. 注意:这里{}m k kx 0)(=φ可根据实际来选择,权系数iw 的选取更是灵活多变的,有时可选取1=i w ,或nw i 1=,对于nw i1=,则相应问题称为均方差的极小化问题.最小二乘拟合函数的求解要使最小二乘问题的目标函数(5.17)达到最小,则由多元函数取得极值的必要条件得),,2,1,0(0m k a Sk==∂∂ 即),,2,1,0(0)()(10m k x y x a w i k ni i m k i k k i ⋅⋅⋅⋅==⎥⎦⎤⎢⎣⎡-∑∑==φφ 亦即),,2,1,0()()()(001m k x y w a x x w n i i k i i j mj n i i k i j i ⋅⋅⋅⋅==⎥⎦⎤⎢⎣⎡∑∑∑===φφφ(5.18)是未知量为ma a a a ,,,,21⋅⋅⋅的线性方程组,称(5.18)式为正规方程组.实际中可适当选择函数系{}m k kx 0)(=φ,由正规方程组解出ma a a a ,,,,210⋅⋅⋅,于是可得最小二乘拟合函数∑==mk kk x a x 0)()(φφ.一般线性最小二乘拟合方法将上面一元函数的最小二乘拟合问题推广到多元函数,即为多维线性最小二乘拟合问题.假设已知多元函数),,,(21nx x x f y ⋅⋅⋅=的一组测量数据);,,,(21iniiiy x x x ⋅⋅⋅),,2,1(m i ⋅⋅⋅=和一组线性无关的函数系{}N k nk x x x 021),,,(=⋅⋅⋅φ,求函数∑=⋅⋅⋅=⋅⋅⋅Nk n k k n x x x a x xx 02121),,,(),,,(φφ对于一组正数mw w w ,,,21⋅⋅⋅,使得目标函数[]2121),,,(∑=⋅⋅⋅-=mi ni i i i i x x x y w S φ达到最小.其中待定系数N a a a a,,,,210⋅⋅⋅由正规方程组),,2,1,0(),(),(0N k y a Nj k j k j⋅⋅⋅==∑=φφφ确定,此处ini i i k mi i k ni i i k mi ni i i j i k j y x x x w y x x x x x x w ),,,(),(),,,(),,,(),(21121121⋅⋅⋅=⋅⋅⋅⋅⋅⋅=∑∑==φφφφφφ注:上面的函数φ关于ia 都是线性的,这就是线性最小二乘拟合问题,对于这类问题的正规组总是容易求解的.如果φ关于ia 是非线性的,则相应的问题称为非线性最小二乘拟合问题.非线性最小二乘拟合方法假设已知多元函数),,,(21nx x x f y ⋅⋅⋅=的一组测量数据);,,,(21iniiiy x x x ⋅⋅⋅),,2,1(m i ⋅⋅⋅=,要求一个关于参数),,2,1,0(N j a j⋅⋅⋅=是非线性的函数),,,;,,,(1021Nn a a a x x x ⋅⋅⋅⋅⋅⋅=φφ对一组正数mw w w ,,,21⋅⋅⋅使得目标函数[]21102110),,,;,,,(),,,(∑=⋅⋅⋅⋅⋅⋅-=⋅⋅⋅mi N ni i i i i N a a a x x x y w a a a S φ达到最小,则称之为非线性最小二乘问题.这类问题属于无约束的最优化问题,一般问题的求解是很复杂的,通常情况下,可以采用共轭梯度法、最速下降法、拟牛顿法和变尺度法等方法求解.实例:黄河小浪底调水调沙问题问题的提出2004年6月至7月黄河进行了第三次调水调沙试验,特别是首次由小浪底、三门峡和万家寨三大水库联合调度,采用接力式防洪预泄放水,形成人造洪峰进行调沙试验获得成功.整个试验期为20多天,小浪底从6月19日开始预泄放水,直到7月13日恢复正常供水结束.小浪底水利工程按设计拦沙量为75.5亿立方米,在这之前,小浪底共积泥沙达14.15亿吨.这次调水调试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪,在小浪底形成人造洪峰,冲刷小浪底库区沉积的泥沙.在小浪底水库开闸泄洪以后,从6月27日开始三门峡水库和万家寨水库陆续开闸放水,人造洪峰于29日先后到达小浪底,7月3日达到最大流量2700立方米/每秒,使小浪底水库的排沙量也不断地增加.下面是由小浪底观测站从6月29日到7月10日检测到的试验数据:表5-1: 试验观测数据单位:水流为立方米/每秒,含沙量为公斤/立方米·84··85·注:以上数据主要是根据媒体公开报道的结果整理而成的,不一定与真实数据完全相符.现在,根据试验数据建立数学模型研究下面的问题:(1) 给出估算任意时刻的排沙量及总排沙量的方法;(2) 确定排沙量与水流量的变化关系.模型的建立与求解对于问题(1),根据所给问题的试验数据,要计算任意时刻的排沙量,就要确定出排沙量随时间变化的规律,可以通过插值来实现.考虑到实际中排沙量应该是随时间连续变化的,为了提高精度,我们采用三次B样条函数进行插值.下面构造三次B样条函数)(x S y =.由试验数据,时间是每天的早8点和晚8点,间隔都是12个小时,共24个点)24,,2,1(⋅⋅⋅=i t i.为了计算方便,令)23,,,1,0(122128⋅⋅⋅=+⎥⎦⎤⎢⎣⎡⋅+-=i i t x i i(5.19)则it 对应于)23,,1,0(1⋅⋅⋅=+=i i x i.于是以)23,,1,0(⋅⋅⋅=i x i为插值节点(等距),步长1=h .其相应的排沙量为)23,,1,0(⋅⋅⋅=i y i 对应关系如下表:·86·表5-2: 插值数据对应关系单位:排沙量为公斤函数)(x S y =所满足的条件为 (1)23,,1,0,)(⋅⋅⋅==i y x S ii;(2) 3500)(,56400)(2223222323231212-=--≈'='=--≈'='x x y y x S y x xy yx S y .取)(x S 的三次B样条函数一般形式为∑-=⎪⎭⎫⎝⎛--=24103)(j j j h x x c x S Ω·87·其中)24,,1,0,1(⋅⋅⋅-=j cj为待定常数,1=h .在这里⎪⎪⎪⎩⎪⎪⎪⎨⎧≥<<+-+-≤+-=Ω2,021,342611,3221)(23233x x x x x x x x x且易知⎪⎪⎪⎩⎪⎪⎪⎨⎧≥±===Ω2,01,610,32)(3x x x x和⎪⎪⎩⎪⎪⎨⎧≥±===Ω'2,01,210,0)(3x x x x 根据B样条函数的性质,)(x S ''在[]23,x x 上连续,则有()∑-=--'='='2413)(j jj xx c x S y Ω由插值条件(1),(2)可得到下列方程组()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=-'=''=-'='⋅⋅⋅==-=∑∑∑-=-=-=23241323024130241323)()(23,,1,0,)(y j c x S y j c x S i y j i c x S j j j j i j j i ΩΩΩ 即⎪⎩⎪⎨⎧'=+-'=+-⋅⋅⋅==++-+-23242311112223,,1,0,64y c c y c c i y c c c i i i i 将232324112,2y c c y c c '+='-=-代入前24个方程中的第一个和最后一个,便可得到方程组F AC =,其中·88·⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⋅⋅⋅⋅⋅⋅=⨯232102424,421410141014124c c c c C A ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡'-'+=3400048000684000458400266626232322100 y y y y y y F显然A 为满秩阵,方程组F AC =一定有解,用消元法求解可得问题的解为56044.39830=c , 4117111.2031=c , 2159510.7882=c , 9189845.6433=c ,1203106.6364=c , 8239727.8115=c ,8249182.1166=c , 1263543.7217=c ,9287842.9988=c , 2302284.2839=c ,4317419.86810=c , 1304836.24311=c ,3307635.15912=c ,6305423.11913=c ,2270672.36214=c ,4240287.43115=c ,0154177.91216=c ,4103000.92017=c ,99818.406218=c , 43725.454719=c ,49279.775020=c ,32155.445221=c , 2098.444222=c ,7450.777923=c ,-450.777924311.2034,2232324011='+=='-=-y c c y c c . 将)24,,1,0,1(⋅⋅⋅-=j c j代入()∑-=--==24131)(j jj x c x S y Ω(5.20)即得排沙量的变化规律.由(5.19)和(5.20)式可得到第i 时间段(12小时为一段)内,任意时刻]12,0[∈t 的排沙量.则总的排沙量为()dt j t c dx x S Y j j⎰∑⎰-=--Ω==284824132411)(经计算可得1110844.1⨯=Y 吨,即从6月29日至7月10日小浪底水库排沙总量大约为1.844亿吨,此与媒体报道的排沙量基本相符.对于问题(2),研究排沙量与水量的关系,从试验数据可以看出,开始排沙量是随着水流量的增加而增长,而后是随着水流量的减少而减少.显然,变化规律并非是线性的关系,为此,我们问题分为两部分,从开始水流量增加到最大值2720立方米/每秒(即增长的过程)为一段,从水流量的最大值到结束为第二段,分别来研究水流量与排沙量的关系.具体数据如表5-3和5-4.表5-3: 第一阶段试验观测数据 单位:水流为立方米/每秒,含沙量为公斤/立方米表5-4: 第二阶段试验观测数据单位:水流为立方米/每秒,含沙量为公斤/立方米对于第一阶段,由表5-3用Matlab作图(如图5-3)可以看出其变化趋势,我们用多项式作最小二乘拟合.·90··91·图5-3设拟合函数为∑==mk kk x a x 1)(φ确定待定常数),,1,0(m k ak=使得211111102])([∑∑∑===⎥⎦⎤⎢⎣⎡-=-=i i i m k k i k i i y x a y x S φ有最小值.于是可以得到正规方程组为m k x y a x mj i k i i j i j k i ,,1,0,0111111 ==⎪⎭⎫⎝⎛∑∑∑===+ 当3=m 时,即取三次多项式拟合,则3,2,1,0,1113111321112111110111==⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛∑∑∑∑∑==+=+=+=k x y a x a x a x a x i k i i i k i i k i i k i i k i求解可得73321108423.1,103172.1,3.1784,-2492.9318--⨯=⨯-===a a a a .于是可得拟合多项式为332213)(x a x a x a a x +++=φ,最小误差为847.72=S ,拟合效果如图所示.·92·图:三次拟合效果,带*号的为拟合曲线.类似地,当4=m 时,即取四次多项式拟合,则正规方程组为4,3,2,1,0111411143111321112111110111==⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛∑∑∑∑∑∑==+=+=+=+=k x y a x a x a x a x a x i ki i i k i i k i i k i i k i i k i求解可得104633210109312.1,1094.1,102626.7,12.0624,-7434.6557---⨯-=⨯=⨯-===a a a a a 于是可得拟合多项式为443322104)(x a x a x a x a a x ++++=φ,最小误差为102.66=S ,拟合效果如图5-5所示.图5-5:四次拟合效果,带*号的为拟合曲线.从上面的三次多项式拟合和四次多项拟合效果来看,差别不大.基本可以看出排沙量与水流量的关系.图5-6:第二段三·93··94· 次多项式拟合效果对于第二阶段,由表5-4可以类似地处理.我们用线性最小二乘法作三次和四多项式拟合.拟合效果如图5-6和5-7所示,最小误差分别为5.459=S 和1.236=S . 从拟合效果来看,显然四次多项式拟合要比三次多项式拟合好的多.图5-7:第二段四次多项式拟合效果。
数学模型第十章插值与拟合方法建模--102数据拟合方法及应用

最小二乘法就是求 c 使得残差平方和
达到最小。
n
Q(c) ( yi f (c, xi )) 2 i 1
2019/9/13
课件
2
1、线性函数
若离散样点 (x1, y1), (x2 , y2 ),,(xn , yn ) 集中在一条直线 的附近,这时可用 y a bx这样的线性函数来拟合。
§2 数据拟合方法及应用
在生产实践和科学研究中,常常有这样的问题:
由实验或测量得到变量间的一批离散样点,要求由
此建立变量之间的函数关系或得到样点之外的数
据。与此有关的另一类问题是数据拟合问题。当原
始数据 (x0 , y0 ), (x1, y1),, (xn , yn ) 有误差时,我们确定的
初等函数 y P(x) 并不要求经过数据点,而是要求在
2019/9/13
课件
13
例 3、某种产品在生产过程中的废品率 y 与它所含
的某种物质量 x 有关,现将试验所得 16 组数据记录
列于下表。
x
34 36 37 38 39 39 39 40
y
1.30 1.00 0.73 0.90 0.81 0.70 0.60 0.50
x
40 41 42 43 43 45 47 48
a
b x
(4) y aebx
2019/9/13
(2)y axb
b
(5) y ae x
课件
(3)y a bln x
(6)
y
a
1 bex
8
b
以
y
ae
x
为例,我们只要作变换
X
数学模型第十章插值与拟合方法建模--101数据插值方法及应用

13.12.2020
课件
6
根据地图的比例,18 mm 相当于 40 km。
根据测量数据,利用 MATLAB 软件对上下边界
进行线性多项式插值,分别求出上边界函数 f2 (x) ,
下边界函数 f1(x) ,利用求平面图形面积的数值积分 方法—将该面积近似分成若干个小长方形,分别求
出这些长方形的面积后相加即为该面积的近似解。
5 188.6
11 191.6
17 458.3
13.12.2020
课件
12
我们将园周展开,借助 MATLAB 软件画出对应的 柱高曲线散点图(下图)。 clear;close; x=linspace(0,2*pi*300,19); y=[502.8 ,525.0,514.3,451.0,326.5,188.6,92.2, 59.6,62.2,102.7,147.1,191.6,236.0,280.5,324.9 ,369.4,413.8,458.3,502.8]; plot(x,y,’o’); axis([0,2000,0,550]);
次样条插值。
13.12.2020
课件
17
例 3、某居民区的自来水是由一个园柱形的水塔提供。 水塔高 12.2 米,直径 17.4 米。水塔由水泵根据塔中水 位高低自动加水,一般每天水泵工作两次。按照设计, 当水塔内的水位降至约 8.2 米时,水泵自动启动加水; 当水位升至约 10.8 米时,水泵停止工作。现在需要了 解该居民区用水规律,这可以通过用水率(单位时间 的用水量)来反映。通过间隔一段时间测量水塔中的 水位来估算用水率。
可以证明当 m n 且 x0 x1 xn 时,这样的多项式 存在且唯一。若要求得到函数表达式,可直接解上 面方程组。
数学建模插值与拟合概要

cz =griddata〔x,y,z,cx,cy,‘method’〕
被插值点 的函数值插值 节点被插值点插值方法
‘nearest’最邻近插值
‘linear’ 双线性插值 ‘cubic’ 双三次插值 'v4'- MATLAB提供的插值方法
缺省时, 双线性插值
要求cx取行向量,cy取为列向量.
▪ %给出〔xi,yj〕点的高程 zij:
▪>>[X,Y]=meshgrid(0:1:20,0:1:20); ▪ % 给出加密的插值坐标网格
第二十五页,共66页。
>>Z=interp2(x,y,z,X,Y,’spline’); %在坐标上进行样条插值
画图: >>clf;%清空图形坐标系中的内容
>>mesh(X,Y,Z) %在网格上画出插值的结果
h=1:0.1:12;
t=interp1(hours,temps,h,'spline'); plot(hours,temps,'+',h,t,hours,temps,'r:')
%作图
xlabel('Hour'),ylabel('Degrees Celsius’)
第十三页,共66页。
第十四页,共66页。
返回
第三十一页,共66页。
%程序一:插值并作海底曲面图
x =[129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5 ];
y =[ 7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5 ];
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
第二讲 插值与数据拟合模型函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二者的数学方法上是完全不同的。
而面对一个实际问题,究竟用插值还是拟合,有时容易确定,有时则并不明显。
在数学建模过程中,常常需要确定一个变量依存于另一个或更多的变量的关系,即函数。
但实际上确定函数的形式(线性形式、乘法形式、幂指形式或其它形式)时往往没有先验的依据。
只能在收集的实际数据的基础上对若干合乎理论的形式进行试验,从中选择一个最能拟合有关数据,即最有可能反映实际问题的函数形式,这就是数据拟合问题。
一、插值方法简介插值问题的提法是,已知1+n 个节点n j y x j j ,,2,1,0),,( =,其中j x 互不相同,不妨设b x x x a n =<<<= 10,求任一插值点)(*j x x ≠处的插值*y 。
),(j j y x 可以看成是由某个函数)(x g y =产生的,g 的解析表达式可能十分复杂,或不存在封闭形式。
也可以未知。
求解的基本思路是,构造一个相对简单的函数)(x f y =,使f 通过全部节点,即),,2,1,0()(n j y x f j j ==,再由)(x f 计算插值,即*)(*x f y =。
1.拉格朗日多项式插值插值多项式从理论和计算的角度看,多项式是最简单的函数,设)(x f 是n 次多项式,记作0111)(a x a x a x a x L n n n n n ++++=-- (1)对于节点),(j j y x 应有n j y x L j j n ,,2,1,0,)( == (2)为了确定插值多项式)(x L n 中的系数011,,,,a a a a n n -,将(1)代入(2),有⎪⎪⎩⎪⎪⎨⎧=++++=++++=++++---n n n n n n nn n n n n n n n n y a x a x a x a y a x a x a x a y a x a x a x a 01110111110001010(3) 记T n T n n n n n n n n n n y y y Y a a a A x x x x x x X ),,,(,),,,(,11110011111100==⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=---- 方程组(3)简写成Y XA = (4) 注意X det 是Vandermonde 行列式,利用行列式性质可得∏≤<≤-=n k j j k x x X 0)(det因j x 互不相同,故0det ≠X ,于是方程(4)中A 有唯一解,即根据1+n 个节点可以确定唯一的n 次插值多项式。
拉格朗日插值多项式实际上比较方便的做法不是解方程(4)求A ,而是先构造一组基函数:n i x x x x x x x x x x x x x x x x x l n i i i i i i n i i i ,,2,1,0,)())(()()())(()()(110110 =--------=+-+- (5) )(x l i 是n 次多项式,满足n j i ji j i x l j i ,,2,1,0,,0,1)( =⎩⎨⎧≠== (6)令∑==n i i i n x l y x L 0)()( (7)显然)(x L n 是满足(2)的n 次多项式,由方程(4)解的唯一性,(7)式表示的)(x L n 的解与(1)式相同。
(5)、(7)称拉格朗日插值多项式,用)(x L n 计算插值称拉格朗日多项式插值。
误差估计插值的误差通过插值多项式)(x L n 与产生节点),(j j y x 的)(x g 之差来估计,记作)(x R n 。
虽然我们可能不知道)(x g 的解析表达式,但不妨设)(x g 充分光滑,具有1+n 阶导数。
利用泰勒展开可以推出,对于任意],[b a x ∈。
),()()!1()()()()(0)1(b a x x n g x L x g x R nj j n n n ∈-+=-=∏=+ξξ (8) 若可以估计1)1(|)(|++≤n n M g ξ (9)则),(,||)!1(|)(|01b a x x n M x R nj j n n ∈-+=∏=+ξ (10) 实际上因为1+n M 常难以确定,所以(10)式并不能给出精确的误差估计。
但是可能看出,n 增加,|)(|x R n 减少;g 越光滑,1+n M 越小,|)(|x R n 越小;x 越接近j x ,|)(|x R n 越小。
例 将区间⎥⎦⎤⎢⎣⎡2,0πn 等分,用x x g y cos )(==产生1+n 个节点,然后作拉格朗日插值多项式。
用)(x L n 计算6cos π(取4位有效数字)。
估计|)(|x R n (取2,1=n )。
解 若1=n ,则)1,0(),(00=y x , ⎪⎭⎫ ⎝⎛=0,2),(11πy x 。
由(5)、(7)式 ππx x x l y l y x L 2102002021)(11001-=--⋅+--⋅=+= 6667.066cos 1=⎪⎭⎫ ⎝⎛=ππL 若2=n ,则)1,0(),(00=y x ,⎪⎭⎫ ⎝⎛=⎪⎭⎫ ⎝⎛=0,2),(,7071.0,4),(2211ππy x y x ,由(5)、(7)式。
⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛---⋅+⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛--⋅+⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⋅=++=4202)4)(0(024042)0(7071.02040241)(2211002ππππππππππππx x x x x x l y l y l y x L ⎪⎭⎫ ⎝⎛-⨯⨯-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-=27071.01624822πππππx x x x8508.066cos 2=⎪⎭⎫ ⎝⎛=ππL 。
估计|)(|x R n :对于x x g cos )(=可设11=+n M ,记节点间隔nh 2π=。
当()1,+∈j j x x x 时 4||||21h x x x x j j <--+ nh h h h x x n j j ⋅⋅⋅⋅<-∏= 324||2于是(10)式给出1112)2)(1(4)1(4324)!1(1|)(|++++=+=⋅⋅⋅⋅+<n n n n n n n h nh h h h n x R π 可以算出6cos π的精确值是0.8660(4位有效数字)⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛6,621ππL L 的误差在|)(|x R n 范围内。
插值多项式的振荡用拉格朗日插值多项式)(x L n 近似))((b x a x g ≤≤虽然随着节点个数的增加,)(x L n 的次数变大,多数情况下误差|)(|x R n 会变小,但n 增加时,)(x L n 的光滑性变坏,有时会出现很大的振荡。
理论上,当∞→n 时,在],[b a 内并不能保证)(x L n 处处收敛于)(x g 。
Runge 给出了一个有名的例子:]5,5[,11)(2-∈+=x x x g 取n j nj x j ,,2,1,0,105 =+-=。
对于 ,6,4,2=n 作)(x L n ,会得到如下图所示的结果。
可以看出,对于较大的||x ,随着n 的增加,)(x L n 的振荡越来越大,事实上可以证明,仅当63.3||≤x 时,才有)()(lim x g x L n n =∞→,而在此区间外,)(x L n 是发散的。
高次插值多项式的这些缺陷,促使人们转而寻求简单的低次数多项式插值。
2. 分段线性插值简单地说,将每两个相邻的节点用直线连起来,如此形成的一条折线就是分段线性插值函数,记作)(x I n ,它满足j j n y x I =)(,且)(x I n 在每个小区间],[1+j j x x 上是线性函数),,1,0(n j =。
)(x I n 可以表示为∑==nj j j n x l y x I 0)()( (12)⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=≤≤--=≤≤--=+++---其它时舍去时舍去,0)(,)0(,)(111111n j x x x x x x x j x x x x x x x x l j j j jj j j j j j j (13))(x I n 有良好的收敛性,即对于],[b a x ∈有,)()(lim x g x I n n =∞→。
用)(x I n 计算x 点的插值时,只用到x 左右的两个节点,计算量与节点个数n 无关。
但n 越大,分段越多,插值误差越小。
实际上用函数表作插值计算时,分段线性插值就足够了,如数学、物理中用的特殊函数表,数理统计中用的概率分布表等。
3. 三次样条插值样条函数的由来分段线性插值虽然简单,n 足够大时精度也相当高。
但是折线在节点处显然不光滑,即)(x I n 在节点处导数不连续。
这影响了它在诸如机械加工等领域(希望插值曲线光滑)中的应用。
所谓样条(Spline),来源于船舶、飞机等设计中描绘光滑外形曲线用的绘图工具。
一根有弹性的细长木条用压铁固定在节点上,其它地方让它自然弯曲,如此画出的曲线称为样条曲线。
因为这种曲线的曲率是处处连续的,所以要求样条函数的二阶导数连续。
人们普遍使用的样条函数是分段三次多项式。
三次样条函数三次样条函数 记作b x a x S ≤≤),(。
要求它满足以下条件:a) 在每个小区间),,1](,[1n i x x i i =-上是3次多项式;b) 在b x a ≤≤上二阶导数连续;c) n i y x S i i ,,1,0,)( ==。
(14) 由条件a ,不妨将)(x S 记为{}n i x x x x S x S i i i ,,1],,[),()(1 =∈=-i i i i i d x c x b x a x S +++=23)( (15)其中i i i i d c b a ,,,为待定系数,共4n 个。
由条件b ,⎪⎩⎪⎨⎧''=''-='='=+++)()(1,,2,1)()()()(111i i i i i i i i i i i i x S x S n i x S x S x S x S (16)容易看出,(14)、(16)式共含有4n -2个方程,为确定)(x S 的4n 个待定参数,尚需再给出2个条件。
最常用的是所谓自然边界条件:0)()(0=''=''n x S x S (17)可以证明,4n 阶线性方程组(14)、(16)、(17)有唯一解,即)(x S 被唯一确定。
但是,这种解法的工作量太大,方程组又常呈病态,所以实际上要设计简便的解法。