第九章插值与拟合.

1、2)恰好给出n - 1个方程(3)0 +a1Xn +a2nanXnyn记此方程组的系数矩阵为1det( A)二A，则2X0 X02X1X1nX0nX1d 21 XnXnnXn第九章插值与拟合插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点， 只要求在某种意义下它在这些点上的总偏差最小。插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。 而面对一个实际问题， 究竟应该用插值还是拟合， 有时 容易确定，有时则并不明显。 1插值方法下面介绍几种基本的、常用的插值：拉格朗日多项式插值、牛顿插值、分段线性插值、

2、Hermite插值和三次样条插值。1.1拉格朗日多项式插值1.1.1插值多项式用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数f(x)在区间a,b上n 1个不同点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次多项式插值就是过 n 1个点(片，f(xj) (i =0,1，n)，作一条 多项式曲线y

3、 = n (x)近似曲线y = f (x)。n次多项式(1)有n 1个待定系数，由插值条件(玄 +a1X +a2X0 + +anX0 = ya()+aM1 +a2X： + +anX； =y1是范德蒙特(Vandermonde)行列式。当x0,x1/ ,xn互不相同时，此行列式值不为零。 因 此方程组(3)有唯一解。这表明，只要 n 1个节点互不相同，满足插值要求( 2)的 插值多项式(1)是唯一的。插值多项式与被插函数之间的差Rn(X)= f(X)- n(X)称为截断误差，又称为插值余项。当f(x)充分光滑时，f(n书化)Rn(X)= f(X)- Ln(X)n l(X),(a,b)(n +1)

4、!n其中 r 1(X)： | 丨(X Xj)。y1.1.2 拉格朗日插值多项式实际上比较方便的作法不是解方程(3)求待定系数，而是先构造一组基函数|.(x)二(XXo) (X Xi J(X Xi J (X Xn)1(Xi Xo)(Xi XiJ(Xi XiG(X Xn)n=nj =oj 7XXjXi _Xj(i 71, ,n)-109-h(x)是n次多项式，满足 li(Xj)二 In XXjn Lj=0 X - Xjlj鼻丿nnLn(X) = yili(X)= yii =0i =0上式称为n次Lagrange插值多项式，由方程(3)解的唯一性，n 1个节点的n次 Lagrange插值多项式存在唯

5、一。1.1.3 用 Matlab 作 Lagrange 插值Matlab中没有现成的Lagrange插值函数，必须编写一个M文件实现Lagrange插值。 设n个节点数据以数组 x0, y0输入(注意Matlat的数组下标从1开始)，m个插值点以数组x输入，输出数组 y为m个插值。编写一个名为lagrange.m的M文件:functiony=lagra nge(x0,y0,x);n=le ngth(x0);m=le ngth(x); for i=1:mz=x(i);s=0.0;for k=1: n p=1.0;for j=1: nif j=k p=p*(z-x0(j)/(x0(k)-x0(j)

6、; endend s=p*y0(k)+s; end y(i)=s;end1.2 牛顿(Newton)插值在导出Newton公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质。1.2.1差商定义 设有函数f(X),X0,X,X2,为一系列互不相等的点，称f (Xi) f (Xj)-(i =j)为f (X)关于点Xi,Xj 一阶差商(也称均差)记为fXi,Xj，即Xi _ XjfXi,Xj称一阶差商的差商f (Xi) - f (Xj)Xi _Xjf Xi,Xj - fXj,XkXi -Xk为f (X)关于点Xi,Xj,Xk的二阶差商，记为fXi,Xj,X k。一般地，称 fXo,Xi, ,

7、Xk- f X!,X2, ,XkXo - Xk为f (X)关于点XoX, ,Xk的k阶差商，记为flXoX!, ,Xk- fX1)X2/ ,Xk fXo,X!,Xk：Xo Xk容易证明，差商具有下述性质：fXi,Xj二 f Xj,XifXi ,Xj,Xk =f Xi,Xk,Xj二 fXj,Xi,Xk1.2.2 Newton插值公式线性插值公式可表成1(X) = f(Xo)+(X-Xo)f Xo,X1称为一次Newt on插值多项式。一般地，由各阶差商的定义，依次可得f(x)二 f (Xo) (x - Xo) fX,XofX,Xo =f Xo* (X - Xj 口纳“fX,Xo,X1二fXo,X

8、1,X2 (X -X2)f X,Xo,X1,X2fX,Xo, ,Xn4H fXo,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) fXo,X1 +(x -Xo)(X -Xj (X -XnG fXo,X1, ,Xn(X -X)(X -Xj (X -Xn) fX,X,X1, X记Nn(X)二 f(X。) (X -Xo)fXo,Xi,Xn(x -Xo)(X - Xi) (X-Xn)f X,Xi,Rn(

9、X)=(X Xo)(X Xi)(X Xn)fX,X,Xi, ,Xn二 n i(X)f X,X,Xi； ,Xn显然,Nn(X)是至多n次的多项式，且满足插值条件，因而它是f (x)的n次插值多项式。这种形式的插值多项式称为Newt on插值多项式。Rn(x)称为Newton插值余项。Newton插值的优点是：每增加一个节点，插值多项式只增加一项，即Nn i(x) =Nn(X)(X -Xo) (X-Xn)fXo,Xi, , X 1因而便于递推运算。而且Newton插值的计算量小于 Lagrange插值。由插值多项式的唯一性可知，Newt on插值余项与Lagra nge余项也是相等的，即Rn(X)

10、=外 1(X)f X,X,X1,Xnf(n1)( )-:n 1(x)(a,b)(n 1)!由此可得差商与导数的关系f (n)(fX,Xi, ,Xnn!其中 (：，J,：二 minxi二 maxXj。0 0 1.2.3 差分当节点等距时，即相邻两个节点之差(称为步长)为常数，Newt on插值公式的形式会更简单。此时关于节点间函数的平均变化率(差商)可用函数值之差(差分)来表 示。定义 设有等距节点xk = x0 + kh(k = 0,1,n)，步长h为常数，fk = f (xk)。称相邻两个节点 xk,xk1处的函数值的增量 fk彳- fk(k = 0,1，n-1)为函数f (x)在 点xk处

11、以h为步长的一阶差分，记为 :fk，即%7 (k=0,1厂,n)类似地，定义差分的差分为高阶差分。如二阶差分为A2fAfk+-Afk(k=0,1，n-2)-般地，m阶差分为左壮十一A(k=2,3,)，上面定义的各阶差分又称为向前差分。常用的差分还有两种:fk = fk - fk4(丄h、(h、沐=fXk +一fXk _、2 丿 2丿称为f (x)在兀处以h为步长的中心差分。称为f(x)在Xk处以h为步长的向后差分;式为般地，m阶向后差分与 m阶中心差分公Jfk.mm_jfk Ofk -严1 -；差分具有以下性质：(i)各阶差分均可表成函数值的线性组合，例如m :j丿m、j(ii )各种差分之间

12、可以互化。向后差分与中心差分化成向前差分的公式如下:mfkf mm .21.2.4等距节点插值公式如果插值节点是等距的，则插值公式可用差分表示。设已知节点x - x0 kh (k = 0,1,2, n)，则有Nn(X)工 fgf X,X1(X -X。)+fX,X1,，Xn(XXo)(XXj(X-XQ、+&fo()(、()-(XXo)n(X Xo)(X X1)(XXn)hn! h若令x = xo th，则上式又可变形为Nn(Xo +th) = fo +tM + +t(tfn*1)n!m% 八 Gi)jj出mJfk 八(-1)jj卫fo (x-Xo)k-2k m _j上式称为Newton向前插值公

13、式。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给出了一个有名的例子：1f(x) 2，x -5,51 +x对于较大的I X I，随着n的增大，J (x)振荡越来越大，事实上可以证明，仅当|x匸3.63时，才有lim Ln(x)二f (x)，而在此区间外，Ln(x)是发散的。nC nn高次插值多

14、项式的这些缺陷，促使人们转而寻求简单的低次多项式插值。1.3.2分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作In (X)，它满足In(Xi)二yi，且In(x)在每个小区间Xi,Xj 1上是线性-1O8-n函数(i = 0,1，,n)。I n(X)可以表示为nln(X)= Th (X)i zz0li (x)= “X - xi 二XiXi JX - Xi 1Xi - Xi 10,x XjxJ (i =0时舍去)X Xi,Xi i (i = n时舍去)其它In(x)有良好的收敛性，即对于 X a,b有，lim I n(x) = f (x)。n_

15、用ln(x)计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。 但n越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就 足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。1.3.3 用Matlab实现分段线性插值用Matlab实现分段线性插值不需要编制函数程序，Matlab中有现成的一维插值函数 interpi。y=i nterp1(x0,y0,x,method)method指定插值的方法，默认为线性插值。其值可为：n earest最近项插值li near线性插值spli ne立方样条插值cubic立方插值。所有的插值方法要求 x0是单调的

16、。当x0为等距时可以用快速插值法，使用快速插值法的格式为*nearest、*linear、*spline、*cubic。1.4 埃尔米特(Hermite)插值1.4.1 Hermite插值多项式如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是Hermite插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的Hermite插值。设已知函数y = f (x)在n 1个互异节点x0,x1，xn上的函数值y f(xi) (i =0,1/ , n)和导数值y： =f(xj (i =0,1，n)，要求一个至多2n 1次的多项式 H

17、(x)，使得H(x) = yiH(x)=yi(i=0,1，n)满足上述条件的多项式H (X)称为Hermite插值多项式。Hermite插值多项式为H(x)其中hi = -j =o1.4.2 用 Matlab 实现二為 hi(Xi -x)(2aiyi - yi)yji卫2X XjXi Xj 丿；1aix Xjj-tHermite 插值Matlab中没有现成的Hermite插值函数，必须编写一个M文件实现插值。设n个节点的数据以数组 x0 (已知点的横坐标)，y0 (函数值)，y1 (导数值) 输入(注意Matlat的数组下标从1开始)，m个插值点以数组x输入，输出数组y为m 个插值。编写一个名

18、为 hermite.m的M文件：function y=hermite(x0,y0,y1,x); n=le ngth(x0);m=le ngth(x);for k=1:myy=0.0;for i=1: nh=1.0;a=0.0;for j=1: n if j=ih=h*(x(k)-x0(j)/(x0(i)-x0(j)A2; a=1/(x0(i)-x0(j)+a;endendyy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);end y(k)=yy;end1.5样条插值许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的

19、凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续， 而且要有连续的曲率，这就导致了样条插值的产生。1.5.1样条函数的概念所谓样条(Spline)本来是工程设计中使用的一种绘图工具，它是富有弹性的细木 条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线(称为样条曲线)，并使连接点处有连续的曲率。三次样条插值就是由此抽象出来的。数学上将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间a,b的一个分划a = x0 :为：:：xn 4 : xn =b如果函数S(x)满足：(i) 在每个小区间xi,xi(i =0,1, ,n -1)上S(x)是m次多项式;(ii) S(x)在a,b上

20、具有m_1阶连续导数。则称S(x)为关于分划的m次样条函数，其图形为 m次样条曲线。显然，折线是一次样条曲线。1.5.2三次样条插值-111-禾U用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线 性插值是一次样条插值。我们只介绍三次样条插值，即已知函数y二f (x)在区间a,b上的n 1个节点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)；(ii

21、i) S(x)在a,b上二阶连续可微。 函数S(x)称为f (x)的三次样条插值函数。由条件(ii)，不妨将S(x)记为S(x)二Sj(x),x Xj ,Xj 1, j =0,1, n _1Sj (x)二 ajx3 bjx2 cjx d j其中aj,bj ,Cj,dj为待定系数，共4n个。由条件(iii)Sj (xj 十)=Sj 十(xj 十)Sj (Xj卅)=S j十(Xj十)j=0,1，门-1(6)Sj (XjHt) = Sj 十(Xj41)容易看出，(5)、(6)式共含有4n-2个方程，为确定 S(x)的4n个待定参数，尚需再 给出2个条件。常用的三次样条函数的边界条件有3种类型：(i)

22、 S(a) = y,S(b)二yn。由这种边界条件建立的样条插值函数称为f (x)的 完备三次样条插值函数。特别地，y； =yn二。时，样条曲线在端点处呈水平状态。如果f(X)不知道，我们可以要求S(X)与f (X)在端点处近似相等。这时以x0, x1, x2, x3为节点作一个三次 Newton插值多项式N a (x)，以xn, xn , xn / ,焉作一 个三次Newt on插值多项式Nb(x)，要求S(a)=N； (a),S(b)=Nb(b).由这种边界条件建立的三次样条称为f (x)的Lagrange三次样条插值函数。(ii) S(a)二 y,S(b) =yn。特别地 yy0时，称为

23、自然边界条件。(iii) S(a P) =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,c on ds),pp=csape(xO,yO,c on ds,

24、valc on ds), y=ppval(pp,x)。其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。 对于三次样条插值，我们提倡使用函数csape, csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppvalopp=csape(x0,y0):使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件，其值可为：complete边界为一阶导数，一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。n ot-a-k not非扭结条件periodic周期条件

25、second边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为0, 0。variati on al设置边界的二阶导数值为 0,0。对于一些特殊的边界条件，可以通过conds的一个1 2矩阵来表示，conds元素的取值为0, 1, 2oconds(i)=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条 件，第二个元素表示右边界的条件，con ds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由 valco nds给出。详细情况请使用帮助help csapeo例1机床加工待加工零件的外形根据工艺要求由一组数据(x

26、, y)给出(在平面情况下)，用程控铳床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x, y)坐标。表中给出的x,y数据位于机翼断面的下轮廓线上，假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求出 x = 0处的曲线斜率和 13乞x空15范围内y的最小值。x03_57911121314_15y01.21.72.02.12.01.81.21.01.6要求用Lagrange、分段线性和三次样条三种插值方法计算。解编写以下程序：x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1

27、 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;y1=lagrange(x0,y0,x); %前面编写的 Lagrange 插值函数y2=i nterp1(x0,y0,x);y3=i nterp1(x0,y0,x,spli ne);pp1=csape(x0,y0);y4=ppval(pp1,x);pp2=csape(x0,y0, sec ond );y5=ppval(pp2,x);X,y1,y2,y3,y4,y5subplot(2,2,1)plot(x0,y0, +,x,y1)title(Lagra nge)subplot(2,2,2) plot(xO,yO, +,x,y2)t

28、itle(Piecewise lin ear)subplot(2,2,3)plot(x0,y0, + ,x,y3) title(Spli ne1)subplot(2,2,4)plot(xO,yO, +,x,y4)title(Spli ne2)dx=diff(x);dy=diff(y3); dy_dx=dy./dx;dy_dxO=dy_dx ytemp=y3(131:151); ymin=mi n( ytemp);in dex=fi nd(y3=ymin); xmin=x(i ndex);xmi n,ymi n计算结果略。可以看出，拉格朗日插值的结果根本不能应用，分段线性插值的光滑性较差(特别是

29、在x =14附近弯曲处)，建议选用三次样条插值的结果。1.6二维插值前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数(曲线)。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点(节点)的 高程(节点值)，为了画出较精确的等高线图，就要先插入更多的点(插值点)，计算这些点的高程(插值)。1.6.1插值节点为网格节点已知 mx n个节点：仅,yj,Zj)( i =12,m; j =1,2,n)，且 xL mA = a1,amT，丫 =(y1,yn)T方程组(9)可表为RTRA 二 RTY( 10)当1(X)，,rm(x)线性无关时，R列满秩，RtR可逆，于是方程组(1

30、0)有唯 一解a =(rt r)rty2.1.2 函数rk(x)的选取面对一组数据(人，比)=1,2/ ,n，用线性最小二乘法作曲线拟合时，首要的、也是关键的一步是恰当地选取a(x)，rm(x)。如果通过机理分析，能够知道y与x之间应该有什么样的函数关系，则r1(x)/ , rm(x)容易确定。若无法知道y与x之间的关系，通常可以将数据(Xi,yj,i =1,2/ ,n作图，直观地判断应该用什么样的曲线去作 拟合。人们常用的曲线有(i) 直线a1x a2(ii) 多项式 y =&必亠amx am彳(一般m = 2,3，不宜太高)(iii)双曲线(一支) y =丑 a2x(iv )指数曲线y =

31、 a1ea2x对于指数曲线，拟合前需作变量代换，化为对a1,a2的线性函数。已知一组数据，用什么样的曲线拟合最好，可以在直观判断的基础上，选几种曲线分别拟合，然后比较，看哪条曲线的最小二乘指标J最小。2.2最小二乘法的 Matlab实现2.2.1解方程组方法 在上面的记号下，J(ai, ,am) =|RA-Y|fMatlab中的线性最小二乘的标准型为Min |RAY2 ,命令为 A = R Y o例4用最小二乘法求一个形如y = a bx2的经验公式，使它与下表所示的数据拟合。x1925313844y19.032.349.073.397.8解编写程序如下x=19253138 44;y=19.0

32、 32.3 49.0 73.3 97.8;r=o nes(5,1),x.A2;ab=ryx0=19:0.1:44;yO=ab(1)+ab (2)*xO.A2;plot(x,y, o ,xO,yO, r)2.2.2多项式拟合方法如果取n(x)，,rm 1(x) =1,x,xm,即用m次多项式拟合给定数据，Matlab中有现成的函数a=polyfit(xO,yO,m)其中输入参数x0,y0为要拟合的数据，m为拟合多项式的次数，输出参数a为拟合多项式 y=amxm+ax+a0系数 a= am,a 1, a 0。多项式在x处的值y可用下面的函数计算y=polyval(a,x) 。例5某乡镇企业1990

33、-1996年的生产利润如下表：年份 1990 1991 1992 1993 1994 1995 1996利润(万元)70 122 144 152 174 196 202试预测1997年和1998年的利润。解作已知数据的的散点图，x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;plot(xO,yO, *)发现该乡镇企业的年生产利润几乎直线上升。因此，我们可以用y二a(x a0作为拟合函数来预测该乡镇企业未来的年利润。编写程序如下：x0=1990 1991 1992 1993 1994 1995 1996;y0

34、=70 122 144 152 174 196 202;a=polyfit(x0,y0,1) y97=polyval(a,1997) y98=polyval(a,1998)求得 a1 =20，a0 二-4.0705 104，1997 年的生产利润 y97=233.4286，1998 年的 生产利润 y98=253.9286。3最小二乘优化在无约束最优化问题中，有些重要的特殊情形，比如目标函数由若干个函数的平 方和构成。这类函数一般可以写成：mF(x) = f：(x),x Rnid：其中x =(花,xn)T，一般假设 m _ n。我们把极小化这类函数的问题：mmin F (x)八 fi2 (x)

35、id：称为最小二乘优化问题。最小二乘优化是一类比较特殊的优化问题，在处理这类问题时，Matlab也提供了一些强大的函数。在Matlab优化工具箱中，用于求解最小二乘优化问题的函数有：Isqlin、 Isqcurvefit、Isqnonlin、Isqnonneg,用法介绍如下。3.1 lsqlin 函数: 2求解 min|Cx - d 2x 2A* x _bs.t.丿 Aeq * x = beqJb Ex Eub其中C, A, Aeq为矩阵，d,b, beq,lb,ub, x为向量。Matlab中的函数为：X=LSQLIN(C,d,A,b,Aeq,be q,L B,UB,X0)例6 用lsqli

36、n 命令求解例4。解编写程序如下：x=19253138 44;y=19.0 32.3 49.0 73.3 97.8;r=o nes(5,1),x.A2;ab=lsqli n( r,y)x0=19:0.1:44;yO=ab(1)+ab (2) *x0.A2; plot(x,y,o,xO,yO,r)3.2 Isqcurvefit 函数给定输入输出数列 xdata, ydata ,求参量x,使得1 III2min F (x, xdata) - ydata 2F (x, xdata,) - ydata,x 2 2 ,Matlab中的函数为X=LSQCURVEFIT(FUN,XO,XDATA,YDATA

37、, LB,UB,OPTIONS) 其中FUN是定义函数 F (x, xdata)的M文件。例7用下面一组数据拟合函数c(t)二 a be0.02kt中的参数a, b, k。t j 100 200 300 400 500 600 700 800 900 1000Cj 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59解该问题即解最优化问题：10I ._0.02 kt j2min F(a,b,k)二 (a be-cj)i d：(1)编写M文件fun 1.m定义函数F(x, tdata):fun ctio n f=fun 1(x,tdata);f=x(

38、1)+x(2)*exp(-0.02*x(3)*tdata);%其中 x(1)=a， x(2)=b，x(3)=k(2 )调用函数Isqcurvefit，编写程序如下：td=100:100:1000;cd=4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59;x0=0.2 0.05 0.05;x=lsqcurvefit(fu n1,x0,td,cd)3.3 Isq nonlin函数已知函数向量Fgulfjx)fk (x) T，求x使得1 2叫n o F(x) 2x 2Matlab中的函数为X=LSQNONLIN(FUN,X0,LB,UB,OPTIONS

39、)其中FUN是定义向量函数 F(x)的M文件。例8用lsqnonlin函数求解例7。解这里F(x)二 F(x,t)二 a beQ02kt1 -g a be.02kt10 -c10Tx = a b k 1(1)编写M文件fun2.m如下：fun ctio n f=fun 2(x);td=100:100:1000;cd=4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59; f=x(1)+x(2)*exp(-0.02*x(3)*td)-cd;(2)调用函数Isqnonlin，编写程序如下：x0=0.2 0.05 0.05;x=ls qnon li n(

40、fun 2,x0)minx；|Cxd3.4 lsq nonneg 函数求解非负的x，使得满足Matlab中的函数为X = LSQNONNEG(C,d,X0,OPTIONS)0.03720.68619 已知C0.28690.70710.85870.17810.62330.62450.0747，求x(x _ 0)满足最小。解.0.6344220.6170.0.8405一1min Cx -d x 2编写程序如下:-127-c=0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170; d=0.8587;0.1781;0.0747;0.8405;x

41、=ls qnonn eg(c,d)4曲线拟合与函数逼近前面讲的曲线拟合是已知一组离散数据(Xj, yj,i =1,，n，选择一个较简单的函数f (x)，如多项式，在一定准则如最小二乘准则下，最接近这些数据。如果已知一个较为复杂的连续函数y(x),a,b，要求选择一个较简单的函数f (x)，在一定准则下最接近f (x)，就是所谓函数逼近。与曲线拟合的最小二乘准则相对应，函数逼近常用的一种准则是最小平方逼近，即b 2J = f(x) y(x) dx(11)达到最小。与曲线拟合一样，选一组函数rk(x),k=1,m构造f (x)，即令f (x) pn(x)amrm(x)代入(11)式，求 印,，am

42、使J达到极小。利用极值必要条件可得一小)3(咕，rj(AJm)门1 (y,rj：a(y,rm)_(12)b这里(g,h) g(x)h(x)dx。当方程组(12)的系数矩阵非奇异时，有唯一解。a最简单的当然是用多项式逼近函数，即选r,(x) =1, r2(x) = x, r3(x) = x2,-b并且如果能使 斤(x)(x)dx =0 , (i = j)，方程组(12)的系数矩阵将是对角阵，计 a算大大简化。满足这种性质的多项式称正交多项式。勒让得(Legendre)多项式是在-1,1区间上的正交多项式，它的表达式为Po(x) =1,可以证明R(x)dkkk2 k! dx2k(x -1) ,k=

43、1,2,10，JR(x)Pj(x)dx 才 22i +12k TkPk1(x)=xPk(x) -Px), k =1,2,k +1k +1常用的正交多项式还有第一类切比雪夫(Chebyshev)多项式Tn (x)二cos(narccosx), (x -1,1, n =0,1,2;)和拉盖尔(Laguerre)多项式d nx dn _x24在H =Span1,x ,x 中的最佳平方逼近多Ln(x) =e (x e ), (x 0, ：), n =0,1,2,)。 dx例 10 求 f (x)二 cosx, x项式。解编写程序如下：syms xbase=1,xA2,xA4;y1=base.*base y2=cos(x)*base. r1=i nt(y1,-pi/2,pi/2) r2=i