用MATLAB进行数据拟合

1、 拟拟 合合 对于情况较复杂的实际问题（因素不易化简，作用对于情况较复杂的实际问题（因素不易化简，作用机理不详）可直接使用数据组建模，寻找简单的因果变机理不详）可直接使用数据组建模，寻找简单的因果变量之间的数量关系，量之间的数量关系， 从而对未知的情形作预报。这样组从而对未知的情形作预报。这样组建的模型为拟合模型。建的模型为拟合模型。 拟合模型的组建主要是处理好观拟合模型的组建主要是处理好观测数据的误差，使用数学表达式从数量上近似因果变量测数据的误差，使用数学表达式从数量上近似因果变量之间的关系。拟合模型的组建是通过对有关变量的观测之间的关系。拟合模型的组建是通过对有关变量的观测 数据的观察、

2、分析和选择恰当的数学数据的观察、分析和选择恰当的数学 表达方式得到的。表达方式得到的。拟拟 合合2.2.拟合的基本原理拟合的基本原理1. 拟合问题引例拟合问题引例拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C) 20.5 32.7 51.0 73.0 95.7电阻电阻R( ) 765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。2040608010070080090010001100 设设 R=at+ba,b为待定系数为待定系数拟拟 合合 问问 题题 引引 例例 2 2 t (h) 0.25 0.5 1 1.5 2 3 4

3、 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形为待定系数kcectckt,)(002468100101102MATLAB(aa1)曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻

4、求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则下与所在某种准则下与所有数据点最为接近，即曲线拟合得最好。有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离拟合与插值的关系拟合与插值的关系 函数插值与曲线拟合都是要根据一组数据构造一个函数作函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同为近似，由于近似的要求不同，二者的数学方法上是完全不同的。的。 实例：实例：下面数据是某次实验所得，希望得到X和 f之间的关系？x1

5、24791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416180510152025已知数据点spline三次多项式插值024681012

6、1416180510152025已知数据点linest三次多项式插值0246810121416180510152025已知数据点nearest三次多项式插值曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第一步: :先选定一组函数先选定一组函数 r1(x), r2(x), rm(x), mn, 令令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1）其中其中 a1,a2, am 为待定系数。为待定系数。 第二步: 确定确定a1,a2, am 的准则（最小二乘准则）：的准则（最小二乘准则）：使使n个点个点（xi,yi) 与与曲线曲

7、线 y=f(x) 的距离的距离 i 的平方和最小的平方和最小 。记记 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 问题归结为，求问题归结为，求 a1,a2, am 使使 J(a1,a2, am) 最小。最小。线性最小二乘法的求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组)( 221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaarrrrrrR112111211,其中其中超定方程一般是不存在解的矛盾方程组。超

8、定方程一般是不存在解的矛盾方程组。 如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。 212211)(imniimiiyararar线性最小二乘法的求解线性最小二乘法的求解 定理：定理：当当R RT TR R可逆时，超定方程组（可逆时，超定方程组（3 3）存在最小二乘解，）存在最小二乘解，且即为方程组且即为方程组 R RT TRa=RRa=RT Ty y的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y 所以，曲线拟合的最小二乘法要解决的问题，实际上就是所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下

9、超定方程组的最小二乘解的问题。求以下超定方程组的最小二乘解的问题。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）线性最小二乘拟合线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中中函数函数rr1 1(x), (x), r rm m(x)(x)的选取的选取 1. 1. 通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2. 2. 将数据将数据 (xi,yi) i=1, n 作图，通过直观判

10、断确定作图，通过直观判断确定 f(x)：用用MATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1. 1. 作多项式作多项式f(x)=a1xm+ +amx+am+1拟合拟合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2. 2. 对超定方程组对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用yRa3.3.多项式在多项式在x x处的值处的值y y可用以下命令计算：可用以下命令计算： y=polyvaly=poly

11、val（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1, am , am+1 (数组数组) )）输入同长度输入同长度的数组的数组X，Y拟合多项拟合多项式次数式次数即要求即要求 出二次多项式出二次多项式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例1 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合xi0.10.20.40.50.60.70.80.91yi1.9783.286.167.347.669.589.489.3011.21）输入以下命令）输入以下命令：x=0:0.1:1; y=-0.447 1.978 3

12、.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; R=(x.2) x ones(11,1)； A=Ry11 11211121xxxxR此时MATLAB(zxec1)解法解法1 1用解超定方程的方法用解超定方程的方法2）计算结果）计算结果： = -9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf1）输入以下命令：）输入以下命令： x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; A=polyfit(x,y,2) z=

13、polyval(A,x); plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形2）计算结果：）计算结果： = -9.8108 20.1293 -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令MATLAB(zxec2)00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf例例2 已知观测数据点如表所示已知观测数据点如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分别用分别用3次和次和

14、6次多项式曲线拟合这些数据点次多项式曲线拟合这些数据点.x=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3 -2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)编写编写Matlab程序如下程序如下:t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)hold onplot(t,s,r-,linewidth,2)plot(t,s,b-,linewidth,2)gridx=0:0.1:

15、1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3 -2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)例例3 用切削机床进行金属品加工时用切削机床进行金属品加工时, 为了适当地调整为了适当地调整机床机床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量刀在一定的时间测量刀具的厚度具的厚度, 得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627

16、.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm解解: 描出散点图描出散点图, 在命令窗口输入在命令窗口输入:t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)解解: 描出散点图描出散点图, 在命令窗口输入在命令窗口输入:t=0:1:16y=30.0 29.1 28.4 28.1 28.0 2

17、7.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)a = -0.3012 29.3804hold onplot(t, y1), hold offa=polyfit(t,y,1)y1=-0.3012*t+29.3804例例3 用切削机床进行金属品加工时用切削机床进行金属品加工时, 为了适当地调整机为了适当地调整机床床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量刀具在一定的时间测量刀具的厚度的厚度, 得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.432

18、8.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm拟合曲线为拟合曲线为: y=-0.3012t+29.3804例例4 一个一个15.4cm30.48cm的混凝土柱在加压实验中的的混凝土柱在加压实验中的应力应力-应变关系测试点的数据如表所示应变关系测试点的数据如表所示1.55 2/ N/m 2/ / N/m 6500 10 33.103 10 2.4761000 10 32.465 10 2. 9361500 10 31.

19、953 10 3. 0362000 10 31.517 10 已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 2.8962375 10 31.219 10 已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 解解 选取选取指数函数指数函数作拟合时作拟合时, 在拟合前需作变量代换在拟合前需作变量代换, 化为化为 k1, k2

20、 的线性函数的线性函数.于是于是,12lnln kk令令0211ln,ln zakak即即01 zaa在命令窗口输入在命令窗口输入:x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*1.0e-6 2375*1.0e-6y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=1.55 2.47 2.93 3.03 2.89plot(x,w,*)y1=exp(8.3009)*x.*exp( -494.5209*

21、x)plot(x,w,*,x,y1,r-)已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 拟合曲线为拟合曲线为:3 -494.52094.0275 10 e0211ln,ln, zakak01 zaa0211-494.5209,ln8.3009, akak3124.0275 10 ,494.5209kk令令则则求得求得于是于是在实际应用中常见的拟合曲线有在实际应用中常见的拟合曲线有:01ya xa直线直线101 nnnya xa xa多项式多项式一般一般

22、n=2, 3, 不宜过高不宜过高.01ayax双曲线双曲线(一支一支) bxyae指数曲线指数曲线二、人口预测线性模型二、人口预测线性模型对于开始提出的实验问题对于开始提出的实验问题, 代如数据，计算得代如数据，计算得27754,1521 aa从而得到人口数与年份的函数关系为从而得到人口数与年份的函数关系为2775415 xy把把x=1999代如，估算出代如，估算出1999年的人口数为年的人口数为 y=1252.1（百万）（百万）12.52亿亿1999年实际人口数量为年实际人口数量为.亿。亿。线性预测模型线性预测模型 英国人口学家英国人口学家MalthusMalthus根据百余年的人口统计资根

23、据百余年的人口统计资料，于料，于17981798年提出了著名的人口自然增长的年提出了著名的人口自然增长的指数增指数增长模型长模型。 btaetx 三、人口预测的三、人口预测的Malthus模型模型 0)0(,xxrxdtdxt rextx0)( 基本假设基本假设 : : 人口人口( (相对相对) )增长率增长率 r 是常数是常数设设x(t) 为为t时刻时刻的的人口人口数量数量,t=0时人口数为时人口数为x0 btatx ln指数增长模型指数增长模型实际中，常用实际中，常用1. 1. 由前由前100100年的数据求出美国的人口增长年的数据求出美国的人口增长Malthus模型模型。2. 2. 预测

24、后预测后100100年（每隔年（每隔1010年）的人口状况。年）的人口状况。3. 3. 根据预测的人口状况和实际的人口数量根据预测的人口状况和实际的人口数量, ,讨论人讨论人口模型的改进情况。口模型的改进情况。美国美国17901790年年19801980年每隔年每隔1010年的人口记录年的人口记录226.5204.0179.3150.7131.7123.2106.592.076.062.9人口人口( (百万百万) )1980197019601950194019301920191019001890年份年份50.238.631.423.217.112.99.67.25.33.9人口人口( (百万百

25、万) )1880187018601850184018301820181018001790年份年份例例解：解： btaetx 设设btax ln2101)ln(),(iiixbtabaJ 取得最小值取得最小值. . 其中其中, , iitxx 表示人口数量表示人口数量。it表示年份表示年份, ,解方程组解方程组: 0ln210ln10211011012101101101iiiiiiiiiiitxbtatbJxbtaaJ即得参数即得参数ba,的值的值.ba,使得使得问题转化为求参数问题转化为求参数 prog41.m %prog41.m % % This program is to predict

26、the number of This program is to predict the number of population %population %format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;187t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;0;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;197t2=1890;1900;1910;1920;1930;1940;1950;1960;1970

27、;1980;0;1980;x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;2;x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;3;204.0;226.5;lnx1=log(x1); lnx2=log(x2);lnx1=log(x1); lnx2=log(x2);a12=sum(t1);a11

28、=10;a21=a12;a22=sum(t1.2); d1=sum(lnx1);d2=sum(lnx1.*t1);A=a11,a12;a21,a22;D=d1;d2;ab=inv(A)*D;disp(a=);disp(ab(1);disp(b=);disp(ab(2);forfor i=1:10 i=1:10 xx1(i)=exp(ab(1)+ab(2) xx1(i)=exp(ab(1)+ab(2)* *t1(i);t1(i);endendforfor i=1:10 i=1:10 xx2(i)=exp(ab(1)+ab(2) xx2(i)=exp(ab(1)+ab(2)* *t2(i);t2(

29、i);endendplot(t1,x1,rplot(t1,x1,r* *-,t1,xx1,b+-, -,t1,xx1,b+-, t2,x2,gt2,x2,g* *-,t2,xx2,m+-);-,t2,xx2,m+-);a= -49.79535457790735b=0.02859807120038仿真结果表明：仿真结果表明： 人口增加的指人口增加的指数模型在短期数模型在短期内基本上能比内基本上能比较准确地反映较准确地反映人口自然增长人口自然增长的规律，但长的规律，但长期预测误差很期预测误差很大，需要修正大，需要修正预测模型。预测模型。拟合曲线拟合曲线原始数据曲线原始数据曲线四、人口预测的四、人口

30、预测的Logistic模型模型人口增长到一定数量后，增长率下降的原因：人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用资源、环境等因素对人口增长的阻滞作用且阻滞作用随人口数量增加而变大且阻滞作用随人口数量增加而变大假设假设r固有增长率固有增长率(x很小时很小时)k人口容量（资源、环境能容纳的最大数量）人口容量（资源、环境能容纳的最大数量）)1()(kxrxr r是是x的减函数的减函数krs 0)( kr)1()(mxxrxxxrdtdx 000)()(xexkkxtxtr 例的例的LogisticLogistic模型留给同学们练习模型留给同学们练习用多项式拟合人口模

31、型用多项式拟合人口模型% % This program is to predict the model of population by This program is to predict the model of population by 4-degree polynomial% 4-degree polynomial% % %prog42.m%prog42.m%format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t1=1790;1800;1810;1820;1830;1840;1850;1

32、860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t=t1;t2;t=t1;t2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;2P2=62.9;76.

33、0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;26.5;P=P1;P2;P=P1;P2;n=4; n=4; % The degree of the fitting polynomial% The degree of the fitting polynomial%a,sa,s=polyfit(t1,P1,n);=polyfit(t1,P1,n);y=polyval(a,ty=polyval(a,t););% % a a is the coefficients vector from n-degree to 0- is the coefficie

34、nts vector from n-degree to 0-degree%degree%plot(t,P,rplot(t,P,r* *-,t,y,b-,t,y,b+-);+-);23a =1.0e+006 * -0.00000000000014 0.00000000107892 -0.00000304878595 0.00381927346813 -1.79012132225427仿真结果表明仿真结果表明, , 人口增加的模型用多项式拟合能人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预比较准确地反映人口自然增长的规律，对长期预测具有指导意义。测具有指导意义。例例2: 2

35、: 海底光缆线长度预测模型海底光缆线长度预测模型某一通信公司在一次某一通信公司在一次施工中施工中, ,需要在水面宽需要在水面宽为为20m20m的河沟底沿直线的河沟底沿直线走向铺设一条沟底光走向铺设一条沟底光缆缆. .在铺设光缆之前需在铺设光缆之前需要对沟底的地形做初要对沟底的地形做初Bix2468101214161820986420ADCih0 x20 x探测到一组等分点位置的深度数据如下表所示探测到一组等分点位置的深度数据如下表所示. .25步探测步探测, ,从而估计所需光缆的长度从而估计所需光缆的长度, ,为工程预算为工程预算提供依据提供依据. .基本情况如图所示基本情况如图所示. .10

36、.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度深度(m)9876543210分点分点2121个等分点处的深度个等分点处的深度(1) (1) 预测通过这条河沟所需光缆长度的近似值预测通过这条河沟所需光缆长度的近似值. .(2) (2) 作出铺设沟底光缆的曲线图作出铺设沟底光缆的曲线图. .024681012141618207891011121314解：解： 用用12次多项式函数拟合光缆走势的曲线图如下次多项式函数

37、拟合光缆走势的曲线图如下仿真结果表仿真结果表明明,拟合曲线拟合曲线能较准确地能较准确地反映光缆的反映光缆的走势图走势图.The length of the label is L= 26.3809 (m)假设所铺设的光缆足够柔软假设所铺设的光缆足够柔软,在铺设过程中光缆触地在铺设过程中光缆触地走势光滑走势光滑,紧贴地面紧贴地面,并且忽略水流对光缆的冲击并且忽略水流对光缆的冲击.% % prog45.m This program is to fit the data by polynomial %prog45.m This program is to fit the data by polynom

38、ial %format longformat longt=linspace(0,20,21);t=linspace(0,20,21);x=linspace(0,20,100);x=linspace(0,20,100);P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.28,13.32,12.61,11.29,10.22,9.1

39、5,7.90,7.95,8.86,9.81,10.80,10.93;80,10.93;a,sa,s=polyfit(t,P,12);=polyfit(t,P,12);yy=polyval(a,xyy=polyval(a,x););disp(yy=);disp(yydisp(yy=);disp(yy););plot(x,yy,rplot(x,yy,r* *-,t,P,b-,t,P,b+-);+-);L=0;L=0;forfor i=2:100 i=2:100 L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2); L=L+sqrt(x(i)-x(i-1)2+(yy(i)-

40、yy(i-1)2);endenddisp(The length of the label is L=);disp(Ldisp(The length of the label is L=);disp(L););1. lsqcurvefit1. lsqcurvefit已知数据点数据点： xdata=xdata=（xdata1，xdata2，xdataxdatan n），）， ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n） 用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合 Matlab Matlab的提供了两个求非线性最小二乘拟

41、合的函数：的提供了两个求非线性最小二乘拟合的函数：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。两个命令都要先建立。两个命令都要先建立M-M-文件文件fun.mfun.m，在其中定义函数在其中定义函数f(x)f(x)，但两者定义，但两者定义f(x)f(x)的方式是不同的的方式是不同的, ,可参可参考例题考例题.最小 ),(21niiiydataxdataxF lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1

42、1），），F F（x x，xdataxdatan n）T T中的参变量中的参变量x(x(向量向量),),使得使得 输入格式为输入格式为: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefit (fun,x0,xdata,ydata,options,grad); (4) x, options = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, options,funval = lsqcurv

43、efit (fun,x0,xdata,ydata,); (6) x, options,funval, Jacob = lsqcurvefit (fun,x0,xdata,ydata,);fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata) 的的M-文件文件, 自变量为自变量为x和和xdata说明：x = lsqcurvefit (fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化 lsqnonlin用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 f(x)f(x)=(f=(f1

44、 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的参量中的参量x x，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22221)()()()()(xfxfxfxfxfnT2. lsqnonlin已知数据点：已知数据点： xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）输入格式为：输入格式为： 1） x=lsqnonlin（fun，x0）；

45、2） x= lsqnonlin （fun，x0，options）； 3） x= lsqnonlin （fun，x0，options，grad）； 4） x，options= lsqnonlin （fun，x0，）； 5） x，options，funval= lsqnonlin （fun，x0，）；说明：x= lsqnonlinlsqnonlin （fun，x0，options）；）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化例例6 已知观测数据点如表所示已知观测数据点如表所示xy03.10

46、.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bx2+cx3 与与已知数据点在最小二乘意义上充分接近已知数据点在最小二乘意义上充分接近.首先编写存储拟合函数的函数文件首先编写存储拟合函数的函数文件.function f=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3保存为文件保存为文件 nihehanshu.m例例6 已知观测数据点如表所示已知观测数据

47、点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bx2+cx3 与与已知数据点在最小二乘意义上充分接近已知数据点在最小二乘意义上充分接近.编写下面的程序调用拟合函数编写下面的程序调用拟合函数.xdata=0:0.1:1;ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;x0=0,0,0;x,resnorm=lsqcurvefit(niheha

48、nshu,x0,xdata,ydata)编写下面的程序调用拟合函数编写下面的程序调用拟合函数.xdata=0:0.1:1;ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;x0=0,0,0;x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata)程序运行后显示程序运行后显示x = 3.0022 4.0304 0.9404resnorm = 0.0912例例6 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.05

49、0.78.560.89.690.911.25113.17求三个参数求三个参数 a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bx2+cx3 与与已知数据点在最小二乘意义上充分接近已知数据点在最小二乘意义上充分接近.说明说明: 最小二乘意义上的最佳拟合函数为最小二乘意义上的最佳拟合函数为f(x)= 3ex+ 4.03x2 + 0.94 x3.此时的残差是此时的残差是: 0.0912.f(x)= 3ex+ 4.03x2 + 0.94 x3.拟合函数为拟合函数为:最小二乘拟合：最小二乘拟合：lsqcurvefit()例例7 假设有一组实测数据假设有一组实测数据xi 0.1 0.2 0

50、.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 yi 2.3201 2.647 2.9707 3.2885 3.6008 3.909 4.2147 4.5191 4.8232 5.1275 假设已知该数据可能满足的原型函数为假设已知该数据可能满足的原型函数为debxaxycx2试求出满足下面数据的最小二乘解的试求出满足下面数据的最小二乘解的a,b,c,d的值。的值。 100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc210102. 0),(minjjktc

51、beakbaFj 例例8 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，kktbeatc2 . 0 . 0)(该问题即解最优化问题：该问题即解最优化问题：MATLAB(fzxec1) 1 1）编写）编写M-M-文件文件 curvefun1.mcurvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中其中 x(1)=a; x(2)=b；x(3)=k;2）输入命令）输入命令tdata=100:100:1000tdata=100:100:1000cdata=cdata=1e-03* *

52、4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tktktbeabea),(10102. 002. 0解法解法1 1.

53、 用命令用命令lsqcurvefitlsqcurvefit3 3）运算结果为）运算结果为：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542 x = 0.0063 -0.0034 0.25424）结论）结论：a=0.0063, b=-0.0034, k=0.2542MATLAB(fzxec2)Tktktcb

54、eacbea),(102. 0102. 0101 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）1）编写编写M-M-文件文件 curvefun2.mcurvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*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)*tdata)- cdata2）输入命令）输入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应是已知参数，故应将将cdata tdata的值写在的值写在curvefun2.m中中3 3）运算结果为）运算结果为 f =1.0e-003 f =1.0e-003 * *(0.2322 -0.1243