9第九章插值与拟合

1、第九章插值与拟合插值：求过已知右限个数据点的近似函数。拟合：已知有限个数据点.求近似函数.不耍求过己知数据点只耍求在某种意义 卞它在这些点上的总偏差最小。插值和拟合都是耍根据一组数据构造一个函数作为近似由近似的耍求不同.二 者的数学方法上是完全不同的。而而对一个实际问题，究竟应该用插值还是拟合，仃时 容易确定.有时则并不明显。§ 1插值方法卜而介绍几种基本的.常用的插值：拉格朗口多项式插值、牛顿插值、分段线性插 值、Hernute插值和三次样条插值。1.1拉格朗口多项式插值1.1.1插值多项式用多项式作为研究插值的工几，称为代数插值。其基本问题是：已知函数/(X)在 区间s，b上&q

2、uot;+ 1个不同点x°z，,x”处的函数frty, =/(x,)(/= 0,1,-,/?),求一个 至多"次多项式久(丫)=為+ %丫+心*(1)使其在给定点处与f(x)同値，即满足插值条件0” (兀)=/(兀)=Z (/ = 0,M)(2)0”(x)称为插值多项式，兀。=0丄,)称为插值节点，简称节点，eb称为插值区间。从儿何上看."次多项式插值就是过刀+ 1个点(”/(兀)(心0丄/),作一条 多项式曲线y = 0”(x)近似曲线y = f(x)。畀次多项式(1)有畀+1个待定系数，由插值条件(2)恰好给出"+ 1个方程a9 + atxQ + a

3、2x + + anx = yQ(3)aQ + alxl +。工+ + a用=帆+勺兀+冬兀+心疋=儿记此方程组的系数矩阵为A，则1det(A) =1是范德蒙(Vandermonde)行列式。当兀刑,心 互不相同时，此行列式値不为零。因 此方程组(3)有唯一解。这表明,只要刃+ 1个节点互不相同,满足插值要求(2)的 插值多项式(1)是唯一的。插值多项式与被插函数Z间的差&(x) = /(x)-久(X)-175-R=x) = f(x)-Ln(x) = 0w(x)，g e W、b)( + l)!其屮 (x) = fl(x-©)。j.Q1.1.2拉格朗口插值多项式实际上比较方便的作

4、法不是解方程(3)求待定系数，而是先构造一组阜换数 (X)=(X_X°)(X_.£t)(X_X田)(百)'(£ 一兀)(兀一x_)(兀一兀+J(兀一入)-177-#- 丿知 J = liMn次多项式，满足. 0令(4)i0i"O上式称为畀次Lagrange插值多项式，由方程(3)解的唯一性，叶+1个节点的"次 Lagrange插值多项式存在唯-。1.1.3 用 Matlab 作 Lagrange 插值Matlab中没仃现成的Lagrange插值函数，必须编打个M文件实现Lagrange插值。 设个节点数据以数组x0,y0输入(注意MiH

5、lai的数组卜标从1开始)，加个插值 点以数组兀输入，输出数组y为7个插值。编写一个名为lagrange.m的M文件： function y=lagrange(x0,y0 r x);n=length(xO);m=length(x);for i=l:mz=x(i)；s=n.o； for k=l:np=l.0;for j=l:nif j=kp-p*(z-xO(j)/(xO(k)-xO(j)；endends=p*y0(k)+s;endy(i)=s；end1.2牛顿（Newton）插值在导IB Newton公式前，先介绍公式农示中所需耍用到的差商、差分的概念及性质。1.2.1差商定义设有函数/（X）,

6、兀州，X，为一系列互不相等的点，称'°一 /（»）（i丰j）为/（X）关J：点兀，心一阶差商（也称均差）记为/兀,X，即 兀一 ®.m=-兀一勺 称一阶差商的差商门兀,X-门®，无兀一无为/（X）关丁点Xi9xj9xk的二阶差商，记为fxj9x-,xko 一般地，称X。一耳为/（x）关于点儿,兀,，兀t的R阶差商，记为fr v i /"丫0，兀，,五/ I/。,兀,，兀 J =X。- H 容易证明，差商具有卜述性质：/比，勺=/【厂,兀fxirxk = fxnxk,xj = fxj9xi,xk1.2.2 Newton插值公式 线性插值公

7、式可表成叭（X） = /U0） + （x-xQ）fxQ, xj称为一次Newton插值筋项式。一般地，由各阶差商的定义，依次叫得/（X）=/（X。）+（X - 兀）/口,兀fg X。 = /X。, X J +（X - X J/ x, XQ,兀fx9 X。，X = / X。, “, X J +（X - £ ）/ x, x0, X, x2/，兀円=/1入内，x” +(x-心)/'x,x°,，兀 将以上各式分别乘以 l(X-Xo),(X-Xo)(X-“),，(x- x0)(x-) (x- xM_1),然后相加并消去两边相等的部分，即得/W = /Uo) +(X - Xo)

8、/x°,xJ + + (兀 _Xo)(X_Xj(X_X”T)/Xo，X,，X+ (X_Xo)(X_Xj(X_X”)/x,Xo, “,£-179-N“ （x） = /（x0） + （x -x0）/x0 £ + + （X Xo）（x 一兀） （X X”.,x” R" （x） = U- x°）（x - 兀）（x - X” ）/x,兀,兀,，兀= ©+】（X）/x，Xo, “，心显然，N“（x）是至多"次的多项式，且满足插值条件，因而它是/（X）的舁次插值 多项式。这种形式的插值多项式称为Newton插值多项式。&©

9、;）称为Newton插值余 项。Newton插值的优点是：每增加一个节点，插值务项式只增加一项，即 g （x） = Nn （x） + （x - X。）（x - x ）f 口。, X,兀+ 因而便递推运算。而且Newton插值的计算昴小J-Lagrange插值。由插值多项式的唯- *性可知，Newton插值余项与Lagrange余项也是相等的，即R”（X）=叫+1（X）/ X, X。, X,X” (" + 1)!©+2)处 S,b)-#-由此可得差商与导数的关系叫)n-#-其屮 § w （a, 0）,a = minx,/7 = maxfx,.。0</<w

10、O£i<n1.2.3差分当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton插值公式的形 己会更简单。此时关J：节点间西数的'卜均变化率（差商）町用函数值Z差（差分）来表 示。定义 设有等距节点勺=x° + M伙=0丄，“），步长/?为常数，齐=/（血）。称相邻两个节点无，如处的函数值的增吊丁叶-£伙=0丄"一 1）为函数/（X）在点丑处以/】为步长的一阶差分，记为皿，即AA = A+i - fk伙=o丄/）类似地，定义差分的差分为高阶差分。如二阶差分为A2/* =乩伙=0丄丿一2） 一般地，加阶基分为yfk=/'fz-

11、y'fk 伙=2,3,),上而定义的齐阶差分又称为向前差分。常用的差分述仃两种:称为/（X）在兀处以力为步长的向后差分：h_2；称为/（切在母处以力为步长的屮心差分。一般地，加阶向后差分与加阶中心差分公 式为w* = vw,-7, - v-7smfk = smlf 丄一 "TA+差分具有以下性质：'(I) 各阶差分均町农成换数値的线性组合，例如+myAV, = z(-Dy . /,戶0J /r =工(-1)丿./,_yj=0J >(II) 各种差分之间可以互化。向后差分与中心差分化成向前差分的公式如下:r=%”，m八k->1.2.4等距节点插值公j如杲插值

12、卩点是等葩的，则插值公式叮用差分衣示。设L1知卩点 xk = x0 + kh (k = 0,1,2,/i),则有N” (x) = /U0) + /x0, x J(x - x°)+ .+/5，",£(尤_兀)(乳_坷)心_£_)=人+牛(兀)+船(兀)(“)心一£丿若令x = x0 + th,则上式乂可变形为N” (兀+巾)=人+。+心-y + 1) ”人 n上式称为Newton向前插值公式。13分段线性插值13.1插值多项式的振荡用Lagrange插值多项式厶(x)近似f(x)(a<x<h).虽然随着节点个数的增加.匚(x)的次数“

13、变人多数情况下误差|恥)|会变小但是"増大时.-(X)的光滑 性变坏，仃时会出现很人的振荡。理论上，当"TS，在“上内并不能保证厶(X)处处收敛P/(x)o Runge给出了一个右名的例子:/« =1l + x2x g -5.5对r较人的| x I,随着的增人，Ln(x)振荡越來越人，事实上可以证明，仅为I x |< 3.63 时，才有lmiLn(x) = /(x),而在此区间外，厶”(x)是发散的。高次插值寒项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。1.3.2分段线性插值简单地说，将每两个相邻的卩点丿U口线连起来，如此形成的条折线就是分段线性

14、插值函数，记作/n(x),它满足仁(兀)=且人(X)在每个小区间兀,上是线性-181-函数(j = o入/) l(x)可以表示为l-OX X-_, x G x、兀(/ = 0时舍去) 兀一 GY Y/,U) = 5 -4-» X e xf.,xi+l (/= “时舍去)兀-兀+10,其它/n(x)有良好的收敛性,即对于xW0,b有,lmi/n(x) = /(x)on->x用/”(切计算X点的插值时，只用到X左右的两个节点，计算最与节点个数无关。 但越人，分段越多，插值谋差越小。实际上用说数农作插值计算时，分段线性插值就 足够了，如数学、物理屮用的特殊苗数表，数理统计中用的概率分

15、布表等。133 用Matlab实现分段线性插值用Matlab实现分段线性插值不需耍编制函数程序，Matlab中冇现成的一维插值西 数 mterpl ov=mterp I(x0.v0,x；method,)method指定插值的方法，默认为线性插值。其值可为：Clearest* 最近项插值"Imeaf线性插值'splmQ 逐段3次样条插值cubic1 保凹凸性3次插值。所有的插值方法要求x0是单调的。当x0为等距时町以用快速插值法，使用快速插值法的格式为neaiesf. -*lmear spline cubic'1.4埃尔米特(Henmte)插ffi1.4.1 Hermi

16、te 插 fff 多项式如果对插值函数不仅耍求它在节点处与函数同值，而比耍求它与函数右相同的一 阶、二阶及至更高阶的导数值，这就是Herinite插值问题。本节主耍讨论在节点处插值 曲数与换数的值及一阶导数值均相等的Hemute插伯。设已知函数y = /(x)在 + 1个互异节点兀,“,X”上的函数值 = /(xj(i = 0,1,)和导数值=(i = 0,l,)，要求一个至多2 + 1次的多项式H(x),使得H(xJ=片 H,(xi)=y,i (/ = 0,1,/?)满足上述条件的多项式H(x)称为Hemike插值多项式。Hennite插值多项式为T80-=工&(兀-x)(2aty.

17、 -yj/«0其中九=11/=oX;-X-181-1.4.2 用 Matlab 实现 Heinute 插值Matlab屮没有现成的Hemute插值肉数，必须编写一个M文件实现插值。设个节点的数据以数31 AO (已知点的横坐标)，yO (函数值)，yl (导数值) 输入(注意Mmlat的数组卜标从1开始)，m个插值点以数组X输入，输出数组y为7 个插值。编写一个名为hemute.m的M文件：function y-hermite(xO z yO x yl,x); n=length(xO);m=length(x);for k=l:myy=O0; for i=l:nh=1.0;a=00;f

18、or j=l:n if j=ih=h*(x(k)-xO(j)/(x0(i)-xO(j)2;a=l/(xO(i)-xO(j)+a;endendyy=yy+h*(xO(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);endy(k)=yy;end1.5样条插值许多工程技术中提出的计算问题对插值函数的光滑性仃较高耍求，如飞机的机翼外 形,内燃机的进、排气门的凸轮曲线,都要求曲线具有较高的光滑程度,不仅要连续, 血IL耍灯连续的曲率，这就导致了样条插值的产生。1.5.1样条函数的概念所谓样条(Spline)本來是工程设计中使用的种绘图工几 它是富仃弹性的细木 条或细金属条。绘图员利用它把一

19、些己知点连接成一条光滑曲线(称为样条曲线)，并 使连接点处何连续的曲率。数学上将貝有一定光滑性的分段多项式称为样条函数。II体地说，给定区的一个分划A: a = xQ <xt < <xn=b如果函数$(x)满足：(i)在每个小区间“,兀(心0丄,打-1)上此丫)是R次多项式；(11)s(x)在。上上具有R-1阶连续导数。則称s(x)为关F分划的R次样条函数，其图形称为R次样条曲线。心兀,，兀称为 样条节点，兀，兀，,兀日称为内节点兀，兀称为边界点，这类样条两数的全体记做S八.£）,称为R次样条隨数空间。显然，折线是一次样条曲线。若$（兀）詁3）,则$（x）是关于分划

20、的R次多项式样条函数。R次多项式样 条函数的-般形式为<0-y-i 代苴中y（i = ol,*）和00 = 1.2,/一 1）均为任意常数，而(x-xy/,x>0, x < Xj>（丿=12,一1）-#-在实际屮址常用的是R = 2和3的情况，即为二次样条函数和三次样条函数。二次样条函数：对Jgb上的分划<<£,=/ 则a"°】q、(5)(6)4(X)= G。+ qx + 寸T + 工寸(X-G SpQ,2),/;=1 /77,(丿=12M-l)。0、 x<x三次样条两数：对J ci.h上的分划:a = Xo <兀&

21、lt; <兀=b,则$«x) = ct0 + qx +守v + 守"+e 5P(A,3)，x |Cv-x.）x>x. 其中（x-xj+=g,（八 1,2,/一1）。P X < x利用样条函数进彳插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值 是一次样条插值。下而我们介绍二次、三次样条插值。1.5.2二次样条函数插值首先，我们注意到鼻wSp（“2）中含彳M + 2个特定常数，故应需耍卄+ 2个插 值条件，因此，二次样条插值问题可分为两类：问题（1）:已知插ffi节点兀和相应的函数值儿Q = o丄）以及端点兀（或兀）处的导数值几（或y人求归（x）

22、 w Sp （A,2）使得严（x,）=y,（i = 0,l,2,/）（7）（几（兀）=儿（或几（£）=儿）问题（2）：已知插值节点兀和相应的导数值y0 = 0.1,2,/）以及端点X。（或x”）处的函 数值y0（或儿），求s2（.v） e SP（A,2）使得（兀）=*（心°2 /）氐兀）=儿（或归（兀）=儿）事实上,可以证明这两类插值问题都是唯一可解的。 对于问题（1）,由条件（7）$2（兀）=% + a$o + 工=儿,、 1t （xj = “ + a 內 +旺=r 1 鬥r归（心）= a° + a$j + -a2x; + 牙工03 一 兀）'=儿（J

23、= 2,3,/）ZL i-iM （x0） = a1 + a2x0 = y,0引入记号X为未知向量，0 =（儿,）1，儿，儿）为己0000*（旺一為）0亍（斗一和亍（斗_兀円）o o是，问题转化为求方程组AX = C的解X = gq，s%eJ的问题,即可得到二次样条两数t（x）的表达式。对于问题（2）的情滉类似。1.53三次样条函数插值由e（x）wS/,Q,3）中含冇 + 3个待定系数，故应需要 + 3个插值条件，己知插值节点兀和相应的两数值/（兀）=;（i =（M,2,/）,这里提供了 +1个条件，还 需要2个边界条件。常用的三次样条函数的边界条件有3种类型：（i）叫（0=儿,叫）=儿由这种边

24、界条件建立的样条插值函数称为/（入-）的 完备三次样条插值函数。特别地，几=冗=0时，样条曲线在端点处呈水平状态。如果/'（X）不知道，我们町以要求几（x）与广（Q在端点处近似相等。这时以 x0,x：,x3为节点作一个三次Newton插值多项式N（l（x）,以暫心十兀宀心乍作一 个三次Newton插值多项式Nh（x）,耍求sa） = N （a）, sb） = N （b）由这种边界条件建立的三次样条称为/（X）的Lagrange三次样条插值函数。（11） s”3（a）=y”o$；（b）=y”,。特别地y” =» = o时，称为口然边界条件。(hi) s(a + 0) = s(b

25、-0),s'(a + 0) = s'(fr-0),(这里要求53(a + 0) = 4少一0)此条件称为周期条件。1.5.4三次样条插值在Matlab屮的实现在Matbb中数据点称Z为断点。如果三次样条插值没仃边界条件，故常用的方法, 就是采用非扭结(notaknot)条件。这个条件强迫第1个和第2个三次名项式的三阶 导数相等。对瑕后一个和倒数第2个三次多项式也做同样地处理。Matlab中三次样条插值也有现成的两数：y=inteip 1 (x0,y0,x/spluie*):y=spline(xO.yO,x): pp=csape(xO,vO,conds), v=ppval(pp,

26、x)。其屮xO.yO是己知数据点，x是插值点，y是插值点的函数值。対三次样条插值，我们提也使用函数csape, csape的返回值是pp形式，耍求插 值点的曲数值，必须调用函数ppvaLpp=csape(xO,yO)：使用默认的边界条件，即Lagrange边界条件。 ppNCsape(xO,vO，conds)屮的conds扌旨定插值的边界条件，苴值可为： complete*边界为一阶导数，即默认的边界条件lot-a-kiiot*非扭结条件periodic*周期条件'second*边界为二阶导数，二阶导数的值0, 0。vanationaV设置边界的二阶导数值为0.0 <>对J

27、 一些特殊的边界条件，可以通过conds的-个1x2矩阵來表示，conds 的 取值为b 2。此时，使用命令pp-csape(x0,v0_extxonds)其屮vO_ext=left. y0, right,衣里left表示左边界的取值,right表示右边界的取值。 conds-j的禽义是给定端点i的丿阶导数，即conds的第一个尤素表示左边界的条 件，第二个尤素表示右边界的条件，conds=2.1表示左边界是二阶导数，右边界是一阶 导数，对应的值由left和light给出。详细情况请使用帮助help csapeo例1机床加工待加丁零件的外形根据匸艺要求由一组数据(兀刃给出(在平而情况卜-),用

28、程控 洗床加工时毎一刀只能沿x方向和y方向走非常小的一步,这就需耍从已知数据得到加 工所要求的步长很小的(x, y)坐标。农1屮给出的数抑;位丁机翼断面的卜轮刖线上，假设盂耍得到x坐标每改变0.】时的y坐标.试完成加丁所需数据，画出曲线，并求出x = 0处的曲线斜率和13 515范用内y的最小值。我1X035791112131415y01.21.72.02.12.01.81.21.01.6耍求用Lagrange.分段线性和三次样条三种插值方法讣算。 解编写以下程序：clC/dearx0=035791112131415;yO=O 1.21.72.02.12.01.81.21.01.6);x=0:

29、0.1:15;yl=lagrange ( xO，yO，x ) ;%调用询面编写的Lagrange插值函数y2=interpl(xO r yO r x); y3=interpl(xO r yO z x,'spline 1); ppl=csape(xO r yO); y4=ppval(ppiz x); pp2=csape(xO,yO # 1 second 1); y5=ppval(pp2 fx); fprintf(»比较一 F不同插值方法和边界条件的结果：n') fprintf(1x yly2y3y4 y5nf)xianshi-tx1l1fy29fy39/y4,ry5;f

30、printf(,xianshi1)subplot(2 92,1), plot(xO z yO r1 + 1z x r yl), title( * Lagrange')subplot(2,2,2), plot(xO r yO rf r x,y2)r title(Piecewise linear 1) subplot(2,2,3), plot(xO/yO, + 1#xzy3)z title( 1Splinel1 )subplot(2 r 2,4)r plot(xO # yO r1 + 1# x r y4)z title( 1Spline21 ) dyxO=ppval ( f nder (p

31、pi) , xO (1) )%求灭0处的导数ytemp=y3(131:151);index=f ind(ytemp=min(ytemp); xymin-x(130+index),ytemp(index)计算结果略。可以看出，拉格朗口插值的结呆根本不能应用，分段线性插值的光滑性较差(特别 是在x=14附近弯曲处)，建议选用三次样条插值的结果。1.6 B样条两数插值方法1.6.1磨光函数实际屮的许筋问题，往往是既耍求近似函数(曲线或曲Iftl) 足够的光滑恂 乂要 求与实际函数有相同的凹凸性.一般插值函数和样条函数都不JV仃这种件质。如果对F 一个特殊皈数进行磨光处理生成磨光函数(多项式).则用辭

32、光皈数构造出样条西数作 为插值函数，既仃足够的光滑性，而且也几仃较好的保凹凸性，因此磨光函数在一维插 值(曲线)和二维插值(曲面)问题中有着广泛的应用由积分理论可知，对于可积函数通过积分会提高函数的光滑度，因此 我们可以利 用积分方法对函数进行磨光处理。定义/W为可积函数,对于力>0,则称积分(0 =制倉山为/(X)的一次磨光函数，力称为磨光宽度。同样的，可以定义/(兀)的R次磨光函数为=不 1(人"df ( k >1 )审实上，憔光函数fk l,(x)比/(X)的光滑程度耍高，乩当磨光宽度力很少时fk h(x) 很接近于/(X)o1.6.2等距B样条曲数对r任意的函数/

33、(Q ,定义It步长为1的中心差分算子方如F：歹3 = /(“+*) -/(x-在此取f(x) = £,则是一个单位方波换数(如图1),记G°(x) = &：。并取力=1,对G°(x)进行一次磨光 得启y+471 ° (1 o_|昭(兀)=£_应0 =L' / +刁_Py 山=一=(X +1)+ - 2兀+ + (x -1)+显然是连续的(如图1)。类似地町得到£次磨光两数为4+1広3=工(一1)，产0实际上，町以证明：Qfc(x)是分段£次多项式，H只旳1阶连续导数，其R阶£ + 导数何R + 2

34、个间断点，记为® = j(丿=0丄2,R + 1)。从而uj知G&(x)是 对应丁分划A:-od<x0<x1 <.<x叶<+s的R次多项式样条函数，称之为基本样条曲数，简称为k次E样条。由J：样条节点为丫产j-尹(；= 0,1,2,-,/:+ 1)是 等距的，故Qa(X)又称为R次等距E样条函数。对J：任意隨数/(X)的R次磨光函数，由归纳法町以得到：局(X)= *匸g#Wx+£)特别地，当f(x) = 1时，有+匚£(宇” =1,从而匚仏住皿=1,且当 k>l时有递推关系1）X+ 2八2R + 1)x +21.6.3

35、一维等距B样条两数插值 等距B样条函数与通常的样条有如卜的关系：h o定理 设有区间“上的均匀分划: x = xQ + jh ( y = 0J,2<-,n ), /?=nJ«n-1则对任意R次样条函数sk(x)eSp(k)都可以衣示为B样条西数族n/ R + 1的线性组合。根据定理，如果已知曲线上一组点(x.,y.),其中 m + jh (/?>0,; = 0J,2,-,n),则町以构造出一条样条麟光曲线(即为B样条西数族的线 性组合)»（x）=jkXX° ;建中C/( j = k,_k + L,n-)为待定常数。用它來逼近曲线，既有较好的精度, 又有

36、良好的保凸性。实际中，垠常用的是R = 3的情况，即一般形式为 喀也(宁7)直中川+ 3个待定系数勺(J = -lA-,n + l)可以由插值条件确定。 对于插值条件k()=yyU=04,2, /») $； (®)=y：C/ = o“)ftrr+1（兀0）= 7工勺0；（-力=儿n /I（9）n+l片（xj= C7Q3（Z-J） = x，f = 0,1,2,/?/-11 zH-1几（Xn） = EC7Q，3（n-=3?，n9注意到Q3(x)的局部非零性及梵换数値：Q3(0) = -, Q3(±l) = -,当Lvl > 2时36Q3(x) = O； II由0

37、；(尤)=2(尤 + *)-2住一*)知，Q(0) = 0, G；(土 1)=干*， 当|x|>2时G；(x) = 0。则(9)式的每一个方程只有三个非零系数，具体为(10)- j + C = 21 叭Ci+4c,+c出=6%,一 J + j = 2hyt由方程组（10）容易求解出-（丿 = 一1,0, + 1）,即可得到三次样条函数片（x） 表达式。1.6.4二维等距B样条函数插值设有空间曲面Z = /（x,y）（未知），如果已知二维等距节点（兀，儿）= （x04-z/j,y0 4- jr） （/,r>0 ）上的值为 ® （, = 0丄,"；J = 0丄,?）

38、，则相应的 B样条磨光曲面的一般形式为刃 T m1$(兀刃=工工5仏<=-1/=-/苴中C“（i = 0丄/;丿=0丄，加）为待定常数，人/可以取不同值，常用的也是k,l = 2和3的情形.这是一种具有良好保凸性的光滑曲面（函数）,在工程设计中是常 用的.但只能使用均匀划分或近似均匀划分的情况。1.7二维插值前面讲述的都是一维插值,即节点为一维变量,插值函数是一元函数（曲线）。若 节点是二维的,插值函数就是二元函数,即曲面如在某区域测屋了若干点（节点）的 高程（节点值人为了画出较粘确的等高线图，就耍先插入更多的点（插值点人计算这 些点的高程（插值）。1.7.1插值节点为网格节点已知7X”

39、 个节点：（兀，儿，令）（i =j =）,且< <xnt:y <<y o求点（x,刃处的插值z。Matlab中有一些计算二维插值的程序。如Z=mteip2（x0.v0.z0,x,y/uietliod,）其屮xOyO分别为加维和n维向最,表示节点,z0为nxm维矩阵,表示节点值,x,y 为一维数组,表示插值点,x与y应是方向不同的向量:,即一个是行向量，另一个是列 向杲，z为矩阵，它的行数为y的维数，列数为X的维数，表示得到的插值，method1 的用法同上面的一维插值。如果是三次样条插值，可以使用命令pp=csape（x0,v0,zO,conds,valconds）&#

40、187; z=fhval（pp.x,v）其中xO,yO分别为加维和"维向量,z0为mxn维矩阵,z为矩阵,它的行数为x的维 数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。例2在一丘陵地带测量高程，才和y方向每隔100米测一个点,得高程如2表,试插 值一曲面，确定合适的模型，并由此找出最高点和该点的高程。解编写程序如卜：clear zclcx=100:100:500;y=100:100:400;z=636697624478450698712630478420680674598412 400662626552334 310;pp=csape(xzy)fz9)xi=100:1

41、0:500;yi=100:10:400czl = fnval(pp,xi yi) cz2=interp2(x/y/z/xi/yi,1 spline 1) i z j=find(czl=max(max(czl) x=xi(i),y=yi(j),zmax=czl(i f j)表201002003001005001006366976244784502006987126304784203006806745984121001006626265523343101.7.2插值节点为散乱节点已知个节点：（x,ynzf）（/ = 1,2,，力）求点（X,刃处的插值z。对上述问题，Matlab中提供了插值换数gn

42、ddata, K格式为：ZI = GRIDDATA（X 丫 乙 XI,YI）其屮X. Y、Z均为n维向彊，指明所给数据点的横坐标、纵坐标和竖坐标。向最XI、 YI是给定的网格点的横坐标和纵坐标.返回値ZI为网格（XL YI）处的两数值。XI 与YI应是方向不同的向最，即一个是行向最.另一个是列向最。例3在某海域测得一些点（x,y）处的水深z由卜表给出，在矩形区域（75,200） x（-5O,15O）内画出海底曲面的图形。X129140103.5SS185.519510515.5107.57781162162117.5y7.5141.52314722.5137.585.5-6 5-81356.5

43、-66 581-33.5Z4S6S6889988949衣3解编写程序如卜：x=129 110 103.5 88 185.5 y=7.5 141.5 23147 22.5zT4 8 68 6195 105 157.5 107.5 77137.5 85.5 -6.5 -8138899881 162 162 117.5;56. 5 -66. 5 84 -33. 5;8949;xi=75:1:200;yi二-50:1:150; zi-griddata(x,y, z,xi5yi*cubic,) subplot( 1,2,1), plot(x.y；*t) subplot( 1,2,2), mesh(xi,

44、yi,zi)§2曲线拟合的线性最小二乘法2.1线性最小二乘法曲线拟合问题的提法是，已知一组（二维）数据，即'F面上的打个点（兀,yj，r = 1,2,/?, x,互不相同，寻求一个两数（曲线）y = f（x）,使/（兀）在某种准则卜' 与所有数据点最为接近，即曲线拟合得最好。-191-线性最小二乘法是解决曲线拟介最常用的方法，某本思路是，令/（X）=勺斤（X） +（x） + + amrm（X）（11）其屮7；（X）是事先选定的一组线性无关的函数，cik是待定系数伙=1,2,加"V“）。 拟介准则是使,心129/（兀）的距离Q的平方和最小,称为最小二乘准则。2

45、.1.1系数冬的确定记丿 a，"”）=!> = £/（£）-），丁<12）t«lf=l为求©，心使丿达到垠小，只需利用极值的必耍条件 = 0（ = 1,/?）,得到关于,4,”的线性方程组nm工 ri (兀)工绞以(兀)一 X = 0, (/ =】, m) /»1A-1(13)mnn工【工匚（兀）（兀）=工匚（兀）兀，（丿=1,"） *1/!/!人（忑）几（忑）nxiwa=u，,订，丫=（儿，,儿广方程组（13）可表为(14)RrRA = RTY当01（x），J；”（x）线性无关时，人列满秩，疋尺可逆，丁是方程组

46、（14）冇唯 一解A = （RtR）'RtY2.1.2函数4 （x）的选取面对一组数据=,/?,用线性址小二乘法作曲线拟合时，首耍的、也是关键的一步是恰当地选取7；（x）,匚（切。如果通过机理分析，能够知道y与xZ 间应该有什么样的函数关系，则“（兀）,匚（兀）容易确定。若无法知道y与兀之间的費 系，通常川以将数据（兀）丿=12/作图，苴观地判断应该用什么样的曲线去作 拟合。人们常用的曲线有（i）直线 二勺兀+（11）多项式 y = aYxm + + amx + am+l （般加= 2,3,不宜A高）（in）双曲线（一支），=鱼+心 叫90(iv)指数曲线 y = ateazX对丁指数

47、曲线，拟介前需作变帚代换，化为对©，的线性两数.已知一组数据，用什么样的曲线拟合最好，町以在占观判断的基础匕选儿种曲线 分别拟合，然后比较，看哪条曲线的最小二乘扌呂标丿最小。2.2最小二乘法的Matlab实现2.2.1解方程组方法在上面的记号下，丿，吗Matlab中的线性最小二乘的标准型为 呼阿一帖 命令为A = RY o便它与表4所示的数据例：用最小三乘法求二个够 如 y=a0+a1 sin(x)+a2cos(x) 的经验公式，x=1:10,y=rand (10”求 ax=1：10f; y=rand(10,1);matr=on es(10,1),si n(x),cos (X)；cs

48、1=matry例4用垠小二乘法求一个形如y = a + 的经验公式， 拟合。4X 1925313844y19032349073397?8解编写程序如卞x=1925313844 1 ;y=19.032.349.073.397.8 1 ;r«ones(5 #1)z x.A2 ;ab=ryx0=19:0.1:44;y0=ab(1)+ab(2)*x0八2;plot(xryz o* zx0,y0z fr 1 )2.2.2多项式拟介方法如果取9(x),，G+(x) = 1,x,x”',即用加次多项式拟介给定数据,Matlab 中有现成的函数a=polvfit(xO,vOan)比中输入参数

49、x0,v0为耍拟合的数据，m为拟合多项式的次数，输出参数a为拟合多项 式 V=amXm+* +aix+ac 系数 a= am, » a：t ao。多项式在x处的值y可用卜面的函数计算y=polyval(a, x)o例5某乡镇企业1990-1996年的生产利润如表5c衣5年份1990199119921993199419951996利润(万元)70122A MM144152174196202试预测1997年和1998年的利润。解作已知数据的的散点图，x0=1990199119921993199419951996;y0=70122144152174196202;plot(xOyO/* *

50、1 )-191-发现该乡値企业的年生产利润儿乎直线上升。因此，我们可以用v = arv + a0作为 拟合两数来预测该乡镇企业未來的年利润。编写程序如卜： x0-1990199119921993199419951996;y0=70122144152174196202;a=polyf it(xO r yO z1) y97=polyval(a z1997 ) y98=polyval(a z1998 )求得 6=20, a0 =-4.0705 xlO4 , 1997 年的生产利润 y97=233. 4286, 1998 年的 生产利润 y98=253. 9286。§3最小二乘优化在无约束最

51、优化问JB中,有些重要的特殊情形,比如目标函数由若干个函数的平 方和构成。这类函数-般可以写成：mi-1其屮x = (x1,.一般假设m > n o我们把极小化这类丙数的问题：mmm F(x)=X/W1-1称为最小二乘优化问题。最小二乘优化是一类比较特殊的优化问题.在处理这类问题时，Matlab也提供了 一些强人的函数。在Matlab优化匸八箱屮，用J：求解赧小二乘优化问題的函数伺：lsqlin、 lsqcuivefit lsqnonlHK lsqnonneg» 用法介绍如卜。3.1 lsqlm 函数求解 mmCx-dA*x <bs.t. Aeq *x = beq lb

52、<x< ub»其中C, A, Aeq为矩阵,d,b,beq,lb,ub,x为向量.Mat lab中的函数为孑x=lsqlin (C, d, A, b, Aeq, beq, lb, ub, xO)例6用lsqlin命令求解例4.解编写程序如2x二1925313844'y=19.032.349.073.397.8'r二ones (5, 1), x. *2;ab=lsqlin(r, y)xO二19:0. 1:44;yO=ab (1) +ab (2) *x0. 2;plot (x, y,' o', xO, y0,' r )F(x, .dat

53、a) - ydat= £ 工(F(x,xdaa) -平呦)3. 2 lsqcurvefit 丿两数给泄输入输出数列idata,ydata,求参;乳使得 1minx 2Mat lab中的隨数；X二LSQCURVEFIT (FUN, XO, XDATA, YDATAjXUBfOPTIONSk H；中FUN是定义函数F(x,xdma)的M文石U 碩 蔓另例7用卜面表6屮的数据拟介丙数W) =辟中的参数a、b、k。农6tj1002003004005006007008009001000Cj4. 544. 99 5. 355. 655. 906.106.266. 396. 50659解该问题即解

54、绘优化问题：mm F(ci,b,k)=工(a + be°Q2kti - cy):1=1(1) 编写M文件funl. m定义西数F(xytdata): function f=funl(xjdata);fx(l)+x(2)*exp(-0.02*x(3)*tdata); %其中 x(l)-a> x(2)-b, x(3)-k(2) 调用函数lsqcurvefit.编写程序如卜：td=100:100:1000;cd二4. 544.99 5. 35 5. 655. 90 6. 10 6. 26 6. 39 6. 50 6. 59;x0= 0. 2 0. 05 0. 05;x=lsqcurvefit(funlt xO, td, cd) 3. 3 lsqnonlin 函数已知负数向S：F(x) = /(x),，