08第8章--插值与拟合.docx

第8章 插值与拟合在工程和科学实验中，变量之间往往存在着固有的函数关系，但这种关系经常很难有明显的解析表达式，通常只能由观察与测试得到一些离散数值。即使有时能给出解析表达式，却因结构过于复杂，不仅不便于使用而且不易于进行计算与理论分析。解决这类问题的方法通常是寻求固有函数的近似逼近，而近似函数的产生办法则因观测数据与背景要求的不同而不同，最常用的两种方法是数据插值与数据拟合。8.1 插值方法8.1.1 一般多项式插值1.多项式插值的一般提法对于未知函数，已知它在区间上个观测点，要求一个至多次多项式，使其在给定点处与有相同的值，即满足条件：，(8.1)称为插值多项式，称为插值节点，简称节点，称为插值区间，式称为插值条件。 从几何上看，次多项式插值就是过个点，作一条多项式曲线近似代替未知函数曲线。2. Lagrange插值公式（1）一阶Lagrange插值设已知及，为不超过一次多项式且满足插值条件，。在几何上为过点和的直线，从而得到.(8.2)为了便于推广到高阶情形，将式变形为对称形式：.（2）二阶Lagrange插值设已知及，为不超过二次的多项式，且满足，和。经计算得到其二阶Lagrange插值多项式，又称为抛物线插值多项式：.（3）阶Lagrange插值按和的求解方法，进行推广可以得到阶Lagrange插值的公式：.3. Newton插值在导出Newton插值公式前，先介绍公式表示中所需要用到的差商概念。（1）函数的差商设有函数及一系列相异的节点，则称为函数关于节点的一阶差商，记为，即.称一阶差商的差商为关于点的二阶差商，记为。一般地，称为关于点的阶差商，记为.（2）Newton插值公式由于关于两节点的线性插值多项式为，可将其表示成，称为一次Newton插值多项式。一般地，由各阶差商的定义，依次可得， ，将以上各式分别乘以，然后相加并消去两边相等的部分，即得 记，显然，是至多次的多项式，且满足插值条件，因而它是的次插值多项式。这种形式的插值多项式称为Newton插值多项式。称为Newton插值余项。Newton插值的优点是：每增加一个节点，插值多项式只增加一项，即，因而便于递推运算。而且Newton插值的计算量小于Lagrange插值。8.1.2 分段线性插值1.插值多项式的振荡用Lagrange插值多项式近似，虽然随着节点个数的增加，的次数变大，多数情况下误差会变小。但是增大时，的光滑性变坏，有时会出现很大的振荡。理论上，当，在内并不能保证处处收敛于。Runge给出了一个有名的例子：，对于较大的，随着的增大，振荡越来越大，事实上可以证明，仅当时，才有，而在此区间外，是发散的。由于高次插值多项式的这些缺陷，也说明并不是插值次数越高效果越好，这就促使人们转而寻求简单的低次多项式插值方法，分段线性插值方法就是一种有效的方法。2.分段线性插值对于未知函数，已知它在区间上个观测点，这里，要寻求一个近似函数。将上的每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作。它满足插值条件，且在每个小区间上是线性函数。事实上，分段线性插值函数可以表示为，其中， 有良好的收敛性，即对于，有。用计算点的插值时，只用到左右的两个节点，计算量与节点个数无关。但越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等都是用分段线性插值计算出来的。8.1.3 样条插值在工程技术实践中，所提出的很多问题都是对插值函数的光滑性有较高要求。如飞机的机翼外形，内燃机的进、排气门的凸轮曲线等，都要求曲线（曲面）具有较高的光滑程度。即不仅要曲线（曲面）连续，而且要求有连续的曲率，即要有一定的光滑性，这就是样条插值要解决的问题。1.样条函数的概念所谓样条原本是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线，并使连接点处有连续的曲率，即要有充分的光滑性。在数学上，将具有一定光滑性的分段多项式称为样条函数。具体地说，给定区间上的一个分划.如果函数满足：（1）在每个小区间上是次多项式；（2）在上具有阶连续导数。则称为关于分划的次样条函数，其图形为次样条曲线。显然，折线是一次样条曲线。2.一维数据的三次样条插值对于区间上的（未知）函数，已知它在个节点上的函数值为，求插值函数使得（1）；（2）在每个小区间上是三次多项式，记为；（3）在上二阶连续可微。则函数称为的三次样条插值函数。由条件（2），不妨将记为，其中，其中为待定系数，共个。由条件（3）.(8.3)容易看出，为确定的个待定参数，尚需再给出2个边界条件。于是有三种常用的三次样条函数的边界条件：（1），。由这种边界条件建立的样条插值函数称为的完备三次样条插值函数。特别地，时，样条曲线在端点处呈水平状态。如果不知道，我们可以要求与在端点处近似相等。这时以为节点作一个三次Newton插值多项式，以作一个三次Newton插值多项式，要求，.由这种边界条件建立的三次样条称为的Lagrange三次样条插值函数。（2），。特别地时，称为自然边界条件。（3），此条件称为周期条件。3.二维数据的双三次样条插值对于二维数据的插值，首先要考虑两个问题：一是二维区域是任意区域还是规则区域；二是给定的数据是有规律分布的还是散乱、随机分布的。第一个问题比较容易处理，只需将不规则区域划分为规则区域或扩充为规则区域来讨论即可。对于第二个问题，当给定的数据是有规律分布时，方法较多也较成熟；而给定的数据是散乱、随机分布时，没有固定的方法，但一般的处理思想是：从给定的数据出发，依据一定的规律恢复出规则分布点上的数据，转化为数据分布有规律的情形来处理。当给定的二维数据在规则区域上有规律分布时，一种常用的插值方法是所谓的双三次样条插值。其基本思想是，对于给定函数的二维数据如表8.1所示。它们对轴和轴的分割：可以导出平面上矩形区域的一个矩形网格分割，如图8.1所示。表8.1 二维规则数据图8.1 矩形网格分割如果令，则所谓双三次样条插值，就是求一个关于和都是三次的多项式，使其满足（1）插值条件：；（2）在整个上，函数的偏导数都是连续的；并称多项式为双三次样条插值函数，其插值公式为.(8.4)实际上，双三次样条插值函数是由两个一维三次样条插值函数作直积产生的。对任意固定的，是关于的三次样条函数；同理，对任意固定的，是关于的三次样条函数。8.1.4 用MATLAB求解插值问题1.一维插值（1）interp1函数MATLAB中一维插值函数interp1的调用格式为：vq=interp1(x0,y0,xq,method,extrapolation)其中x0为已知的插值节点，y0是对应于x0的函数值，xq是欲求函数值的节点坐标，返回值vq是求得的节点xq处的函数值，method指定插值的方法，默认为线性插值，常用的取值有：nearest 最近邻插值；linear 线性插值；spline 三次样条插值，函数是二次光滑的；cubic 立方插值，函数是一次光滑的。extrapolation是外插策略。（2）三次样条插值函数csape三次样条插值还可以函数csape，csape的返回值是pp形式。其调用格式如下：pp=csape(x0,y0)使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件，其值可为complete 边界为一阶导数，一阶导数的值在valconds参数中给出。not-a-knot 非扭结条件。 periodic 周期条件。second 边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为0, 0。variational 设置边界的二阶导数值为0,0。对于一些特殊的边界条件，可以通过conds的一个矩阵来表示，conds元素的取值为0，1，2。conds(i)=j的含义是给定端点的阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。利用pp结构的返回值，要求插值点的函数值，需要调用函数fnval。还可以计算返回值函数的导数和积分，命令分别为fnder，fnint，这两个函数的返回值也是pp结构。pp结构相关函数的功能总结于表8.2。表8.2 一维插值的相关函数调用格式函数功能pp1=csape(x0,y0)计算插值函数pp2=fnder(pp1)计算pp1对应函数的导数，返回值pp2也是pp结构pp3=fnint(pp1)计算pp1对应函数的积分，返回值pp3也是pp结构y=fnval(pp1,x)计算pp1对应的函数在x点的取值例8.1 已知欧洲一个国家的地图，为了算出它的国土面积和边界长度，首先对地图作如下测量：以由西向东方向为轴，由南向北方向为轴，选择方便的原点，并将从最西边界点到最东边界点在轴上的区间适当地分为若干段，在每个分点的方向测出南边界点和北边界点的坐标和，这样就得到了表8.3的数据（单位：mm）。表8.3 某国国土地图边界测量值（单位：mm）7.010.513.017.534.040.544.548.056.0444547505038303034445970729310011011011061.068.576.580.591.096.0101.0104.0106.5363441454643373328117118116118118121124121121111.5118.0123.5136.5142.0146.0150.0157.0158.0326555545250666668121122116838182868568根据地图的比例我们知道18mm相当于40km，试由测量数据计算该国国土边界的近似长度和近似面积，并与国土面积的精确值41288km2比较。解 该地区的示意图见图2。图2 区域边界示意图若区域的下边界和上边界曲线的方程分别为，则该地区的边界线长为，计算时用数值积分即可。计算该区域的面积，可以把该区域看成是上、下两个边界为曲边的曲边四边形，则区域的面积.计算相应的数值积分就可求出面积。为了提高计算精度，分别对上、下边界用分段线性插值和三次样条插值计算。利用分段线性插值计算时，得到边界长度的近似值为1030.3km；区域面积的近似值为42407 km2，与其准确值41288km2的相对误差为2.71%。利用三次样条插值计算时，得到边界长度的近似值为1161.0km，区域面积的近似值为42476km2，与其准确值41288km2的相对相差为2.88%。计算分段线性插值的MATLAB程序如下：clc, clear, close all, S0=41288;a=load(gdata8_1.txt);x0=a(1:3:end,:); x0=x0; x0=x0(:); %提取点的横坐标y10=a(2:3:end,:); y10=y10; y10=y10(:); %提出下边界的纵坐标y20=a(3:3:end,:); y20=y20; y20=y20(:);plot(x0,y10,*-) %画下边界曲线hold on, plot(x0,y20,.-) %画上边界曲线x1=min(x0):max(x0); %插值节点y1=interp1(x0,y10,x1); y2=interp1(x0,y20,x1); %分段线性插值L1=trapz(x1,sqrt(1+gradient(y1,x1).2) %计算下边界的长度L2=trapz(x1,sqrt(1+gradient(y2,x1).2) %计算上边界的长度L=L1+L2; %计算地图上边界的长度LL1=L/18*40 %计算实际的边界长度S1=trapz(x1,y2-y1); %计算地图上的近似面积SS1=S1/182*1600 %计算实际面积的近似值delta1=abs(SS1-S0)/S0 %计算面积的相对误差计算三次样条插值的MATLAB程序如下：clc, clear, S0=41288;a=load(gdata8_1.txt);x0=a(1:3:end,:); x0=x0; x0=x0(:); %提取点的横坐标y10=a(2:3:end,:); y10=y10; y10=y10(:); %提出下边界的纵坐标y20=a(3:3:end,:); y20=y20; y20=y20(:);pp1=csape(x0,y10); pp2=csape(x0,y20); %计算三次样条插值函数dp1=fnder(pp1); dp2=fnder(pp2); %求三次样条插值函数的导数L2=quad(x)sqrt(1+fnval(dp1,x).2)+sqrt(1+fnval(dp2,x).2),x0(1),x0(end)L2=L2/18*40 %换算成边界的实际长度S2=quad(x)fnval(pp2,x)-fnval(pp1,x),x0(1),x0(end) %计算地图上面积SS2=S2/182*1600 %计算实际面积的近似值delta2=abs(SS2-S0)/S0 %计算面积的相对误差2.二维插值（1）插值节点为散乱节点已知个节点，求点处的插值。对上述问题，MATLAB中提供了插值函数scatteredInterpolant，其调用格式为：Fz=scatteredInterpolant(x0,y0,z0,Method,ExtrapolationMethod);其中返回值Fz是结构数组，相当于给出了插值函数的表达式；x0,y0,z0分别为已知个点的坐标；Method是插值方法；ExtrapolationMethod是区域外部节点的外插方法，ExtrapolationMethod的值可以取为linear和nearest，默认值为linear。要计算插值点（x,y）处的值，调用函数Fz即可，z=Fz(x,y)就求出所要求点的插值。（2）插值节点为网格节点已知个节点，且；，要求在点处的插值。MATLAB提供了一些计算二维插值的命令。如z=interp2(x0,y0,z0,x,y,method)其中x0，y0分别为维和维向量，表示节点，z0为维矩阵，表示节点值，x，y为一维数组，表示插值点，x与y应是方向不同的向量，即一个是行向量，另一个是列向量，z为矩阵，它的行数为y的维数，列数为x的维数，表示得到的插值，method的用法同上面的一维插值。如果是三次样条插值，可以使用命令pp=csape(x0,y0,z0,conds,valconds)，z=fnval(pp,x,y)其中x0，y0分别为维和维向量，z0为维矩阵，为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。例8.2 已知平面区域，的高程数据见表8.4（单位：m）。求该区域地表面积的近似值，并画出该区域的三维表面图。表8.4 高程数据表12001100100090080070060050040030020010001350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 1501370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 2101380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 3501420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 5001430 1450 1460 1500 1550 1600 1550 1600 1600 1600 1550 1500 1500 1550 1550950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 870 900 936 950830 980 1180 1320 1450 1420 400 1300 700 900 850 810 380 780 750740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550650 760 880 970 1020 1050 1020 830 800 700 300 500 550 480 350510 620 730 800 850 870 850 780 720 650 500 200 300 350 320370 470 550 600 670 690 670 620 580 450 400 300 100 150 2500 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400解 原始数据给出的网格节点上的高程数据，为了提高计算精度，我们利用双三次样条插值，得到给定区域上网格节点上的高程。利用分点把剖分成140个小区间，利用分点把剖分成120个小区间，把平面区域，剖分成个小矩形，对应地把所计算的三维曲面剖分成个小曲面进行计算，每个小曲面的面积用对应的三维空间中4个点所构成的两个小三角形面积的和作为近似值。计算三角形面积时，使用海伦公式，即设的边长分别为，则的面积。利用MATLAB求得的地表面积的近似值为m2，所画的三维表面图如图8.2所示。图8.2 三维表面图计算和画图的MATLAB程序如下：clc, clear, close allz0=load(gdata8_2.txt);x0=0:100:1400; y0=1200:-100:0;x=0:10:1400; y=1200:-10:0; pp=csape(x0,y0,z0); %双三次样条插值,注意这里z0要转置z=fnval(pp,x,y); z=z; surf(x,y,z) %画三维表面图m=length(x); n=length(y); s=0; for i=1:m-1 for j=1:n-1 p1=x(i),y(j),z(j,i); p2=x(i+1),y(j),z(j,i+1); p3=x(i+1),y(j+1),z(j+1,i+1); p4=x(i),y(j+1),z(j+1,i); p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3); L1=(p12+p23+p13)/2;s1=sqrt(L1*(L1-p12)*(L1-p23)*(L1-p13); L2=(p13+p14+p34)/2; s2=sqrt(L2*(L2-p13)*(L2-p14)*(L2-p34); s=s+s1+s2; endends %显示面积的计算值8.2 数据的拟合方法在实际中，如果已知函数的有限个数据点，并且实际需要求一个近似函数，但不要求过所有的数据点，只要求在某种意义下，近似函数与这些数据点最接近，即在这些数据点上的总偏差最小，则此类问题就是数据拟合问题。插值和拟合都是要根据一组数据构造一个近似函数，对一个实际的问题，究竟应该用插值方法，还是用拟合方法，一般需要依据实际的要求来确定。数据拟合问题的形式多种多样，解决的方法也有许多。在此，只简单介绍基于最小二乘原则的一维数据拟合方法的数学原理，及数据拟合MATLAB工具箱的使用。8.2.1 最小二乘拟合已知一组二维数据，即平面上的个点，互不相同，要寻求一个函数（曲线），使在某种准则下与所有数据点最为接近，即曲线拟合得最好。记，则称为拟合函数在点处的偏差（或残差）。为使在整体上尽可能与给定数据最为接近，可以采用“偏差的平方和最小”作为判定准则，即通过使(8.5)达到最小值。这一原则称为最小二乘原则，根据最小二乘原则确定拟合函数的方法称为最小二乘法。一般来讲，拟合函数应是自变量和待定参数的函数，即.(8.6)因此，按照关于参数的线性性，最小二乘法也分为线性最小二乘法和非线性最小二乘法两类。1.线性最小二乘法给定一个线性无关的函数系，如果拟合函数以其线性组合的形式(8.7)出现，例如，或者.则就是关于参数的线性函数。将式代入式，则目标函数是关于参数的多元函数。由，亦即.可得.(8.8)于是式形成了一个关于的线性方程组，称为正规方程组。记，则正规方程组可表示为.(8.9)由代数知识可知，当矩阵是列满秩时，是可逆的。于是正规方程组有唯一解，即(8.10)为所求的拟合函数的系数，就可得到最小二乘拟合函数。2.非线性最小二乘拟合对于给定的线性无关函数系，如果拟合函数不能以其线性组合的形式出现，例如或者，则就是关于参数的非线性函数。将代入式中，则形成一个非线性函数的极小化问题。为得到最小二乘拟合函数的具体表达式，可通过适当的变换转化为线性最小二乘问题，或者用非线性优化方法求解出参数。3.拟合函数的选择数据拟合时，首要也是最关键的一步就是选取恰当的拟合函数。如果能够根据问题的背景通过机理分析得到变量之间的函数关系，那么只需估计相应的参数即可。但很多情况下，问题的机理并不清楚。此时，一个较为自然的方法是先做出数据的散点图，从直观上判断应选用什么样的拟合函数。一般来讲，如果数据分布接近于直线，则宜选用线性函数拟合；如果数据分布接近于抛物线，则宜选用二次多项式拟合；如果数据分布特点是开始上升较快随后逐渐变缓，则宜选用双曲线型函数或指数型函数，即用或拟合。如果数据分布特点是开始下降较快随后逐渐变缓，则宜选用，或等函数拟合。常被选用的非线性拟合函数有对数函数，S形曲线函数等。8.2.2 数据拟合的MATLAB实现MATLAB提供了两种方法进行数据拟合。一种是以命令行形式进行函数拟合；另一种是以图形界面的形式拟合一元或二元函数。MATLAB中拟合命令很多，我们这里只介绍fit，lsqcurvefit；另外介绍使用用户图形界面的拟合命令cftool。1.fittype和fit函数函数fit需要和函数fittype配合使用，fittype用于定义拟合的函数类，fit进行函数拟合。fit既可以拟合一元或二元线性函数，也可以拟合非线性一元或二元函数，我们这里首先介绍一下这两个函数的调用格式，然后举几个拟合的例子. fittype的调用格式为：aFittype=fittype(libraryModeName) %利用库模型定义函数类aFittype=fittype(expression,Name,Value) %利用字符串定义函数类aFittype=fittype(linearModeTerm,Name,Value) %利用基函数的线性组合定义函数类aFittype=fittype(anonymousFunction,Name,Value) %利用匿名函数定义函数类函数fit的调用格式为：fitobject=fit(x,y, aFittype) %x和y分别为自变量和因变量的观测值，必须为列向量；返回值fitobject为拟合函数的信息。fitobject=fit(x,y,z,aFittype) %x,y为自变量的观测值的两列矩阵，z为因变量的观测值列向量，这里是拟合二元函数。fitobject=fit(x,y, aFittype,fitOptions) % fitOptions指明拟合的算法。fitobject, gof=fit(x,y, aFittype,Name,Value) %返回值gof为结构数组，给出了模型的一些检验统计量；Name是设置某个属性，Value为属性的对应取值。（1）使用库模型拟合函数三次多项式函数：f=fittype(poly3)f = Linear model Poly3: f(p1,p2,p3,p4,x) = p1*x3 + p2*x2 + p3*x + p4一阶Fourier级数：f=fittype(fourier1)f = General model Fourier1:f(a0,a1,b1,w,x) =a0 + a1*cos(x*w) + b1*sin(x*w)例8.3 利用函数，产生一组数据，再在上添加服从标准正态分布的随机数。利用产生的数据拟合一阶Fourier级数和二阶Fourier级数。 clc, clearx0=linspace(1,20,50); x0=x0; y0=1+2*cos(3*x0)+3*sin(3*x0);ya=y0+normrnd(0,1,50,1); %添加随机扰动%利用库模型实际上不需要使用fittypef1,g1=fit(x0,ya,fourier1) %拟合一阶Fourier级数f2,g2=fit(x0,ya,fourier2) %拟合二阶Fourier级数MATLAB工具箱库模型中的函数类是很丰富的，就不一一举例了。可以在命令窗口中输入doc fittype打开帮助页面，再点击“Model Names and Equations”链接，就可以看到库函数中的所有函数类。（2）利用字符串定义函数类例8.4 用模拟数据拟合函数. clc, clearx0=linspace(1,10,20); y0=linspace(3,10,20);z0=x0.2+2*x0.*y0+3*y0.2; %产生z的无扰动数据za=z0+normrnd(0,1,1,20); %产生扰动数据g=fittype(a*x2+b*x*y+c*y2,dependent,z,independent,x, y)f1,g1=fit(x0,y0,z0,g,Start,rand(1,3) %利用无扰动数据拟合f2,g2=fit(x0,y0,za,g,Start,rand(1,3) %利用扰动数据拟合注8.1 定义拟合函数类型时，只需说明自变量属性independent和因变量属性dependent。自变量属性independent的默认值为x，可以不说明；因变量属性dependent的默认值y也可以不说明。例8.5 利用模拟数据拟合函数。我们要拟合一个如下的4阶傅里叶级数，其中。拟合时通过参数下界属性“Lower”和上界属性“Upper”的限制，使得等于常数1。clc, clearx0=1:50; y0=2*cos(2*x0)+3*sin(4*x0);lb=-inf*ones(1,9),1; ub=inf*ones(1,9),1; %w参数在最后一个位置，约束w=1f,g=fit(x0,y0,fourier4,Lower,lb,Upper,ub) （3）利用匿名函数定义函数类例8.6（MATLAB工具箱的帮助例程）利用给定数据拟合分段线性函数clc, clearg=(a,b,c,d,k,x)(a+b*x).*(x=k); %定义匿名函数g=fittype(g) %生成fittype类型的函数类, ,independent的默认属性值x可不说明x=0.81;0.91;0.13;0.91;0.63;0.098;0.28;0.55;0.96;0.96;0.16;0.97;0.96;y=0.17;0.12;0.16;0.0035;0.37;0.082;0.34;0.56;0.15;-0.046;0.17;-0.091;-0.071;f=fit(x,y,g,Start,1, 0, 1, 0, 0.5) %函数拟合plot(f,x,y) %画图注8.2 匿名函数的定义中，必须把自变量放在匿名函数输入参数的最后一个。2.lsqcurvefit函数非线性拟合的主要准则也是最小二乘准则，我们以一元函数为例说明。对于未知函数，其中为未知的参数向量，函数关于是非线性的。给定的一些观测值，拟合参数的最小二乘准则，就是所确定参数的值要使在观测点上误差平方和(8.11)最小，这里。 求式的极小值问题有很多算法，我们这里就不介绍了，例如MATLAB工具箱fit函数使用的默认算法是Levenberg-Marquardt算法，还有Trust-Region算法。MATLAB非线性拟合的主要命令有fit（要用fittype定义函数类），lsqcurvefit，nlinfit等命令. fit函数使用很方便，但只能拟合一元和二元函数，lsqcurvefit可以拟合任意多个自变量的函数，并且可以约束未知参数的下界和上界；nlinfit函数无法约束参数的界限，我们这里就不介绍了。下面首先给出lsqcurvefit的用法说明，再举几个用lsqcurvefit拟合函数的例子。要拟合函数，给定的观测值，的观测值，求参数,使得误差平方和最小，即.MATLAB 中的函数为theta=lsqcurvefit(fun,theta0,xdata,ydata,lb,ub,options)其中fun是定义函数的M文件或者是定义的匿名函数返回值，theta0是的初始值，xdata是自变量的观测值，ydata是函数的观测值，lb是参数的下界，ub是参数的上界，options参数可以对计算过程的一些算法等属性进行设置，返回值theta是所求参数的值。例8.7 用模拟数据拟合一元函数。clc, clearx0=linspace(0,10,20); y0=4*exp(-0.03*x0).*(2*cos(x0)+1); %取a=4,b=0.03,c=2,d=1y=(c,x)c(1)*exp(-c(2)*x).*(c(3)*cos(x)+c(4);o=optimset(MaxFunEvals,4000,MaxIter,4000);c0=lsqcurvefit(y,rand(4,1),x0,y0,o)注8. 3 （1）程序中定义匿名函数时，必须使用“.*”，否则会出错。（2）对于一些非线性拟合，有时感觉像在凑数，MATLAB每次运行的答案是不一样的。例8.8 用模拟数据拟合三元函数。clc, clearx=linspace(0,10,20); y=linspace(2,10,20); z=linspace(0,5,20);w0=4*exp(-0.02*x).*(3*sin(y)+6*cos(z); %取a=4,b=0.02,c=3,d=6w=(c,t)c(1)*exp(-c(2)*t(:,1).*(c(3)*sin(t(:,2)+c(4)*cos(t(:,3);o=optimset(MaxFunEvals,40000,MaxIter,4000);c0=lsqcurvefit(w,rand(1,4),x,y,z,w0,o)3.曲线和曲面拟合的用户图形界面解法可以使用cftool进行曲线和曲面拟合。具体执行步骤如下：（1）把数据导入到工作空间；（2）运行cftool，打开用户图形界面窗口；（3）选择适当的模型进行拟合；（4）把计算结果输出到MATLAB工作空间。可以通过帮助（运行doc cftool）熟悉该命令的使用细节。例8.9 利用数据x=1:50; y=2*cos(2*x)+3*sin(4*x); 拟合4阶Fourier级数，(8.12)即待定参数。使用用户图形界面进行拟合。（1）第一步在命令窗口中运行x=1:50; y=2*cos(2*x)+3*sin(4*x);即把数据加载到工作空间。（2）运行cftool，打开用户图形界面。进行如图8.3所示的操作，首先在左上方选择有关的数据，然后选择拟合的函数为4阶Fo