【tong】数值建模

1、数 值 建 模数据处理的方法及示例数据处理的方法及示例一、一、插值插值二、二、数据拟合数据拟合三、三、神经网络神经网络方法方法四、四、支持向量机支持向量机讲课大纲插 值造飞机的工厂拟 合姚沁蕾的身高 在工程实践与科学实验中,常常需要从一组试验数据之中找到自变量与因变量之间的关系，一般可用一个近似函数表示。 函数产生的办法因观测数据的要求不同而异，数据插值与拟合是两种常用的方法。Xx0 x1xnY=f(x)y0y1ynxyX0357911 12 13 14 15Y0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.61.901.882.082.26回归分析 1855年 “遗传

2、的身高向平均数方向的回归” 父母越高，子女的身高会比一般儿童高； 矮个子的儿子比其父要高，身材较高的父母所生子女的身高将回降到人的平均身高。高尔顿生物统计学派奠基人施化肥量x15202530354045水稻产量y330345365405445450455施化肥量对水稻产量影响的试验数据：插值与拟合的区别 插值与拟合的插值与拟合的区别：区别： 插值插值试图去通过已知点了解试图去通过已知点了解未知点未知点处的函处的函数值数值； 拟合拟合则在于则在于在整体在整体上用某种已知函数去拟上用某种已知函数去拟合合数据数据点列所在未知函数的性态。点列所在未知函数的性态。 关键关键区别区别在于插值要求在于插值要

3、求必须经过必须经过已知点列，已知点列，拟合只求拟合只求尽量靠近尽量靠近不必经过！不必经过！ 内容：插值是离散函数逼近的重要方法，利用它可通过内容：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处函数在有限个点处的取值状况，估算出函数在其他点处的近似值的近似值 要求：掌握要求：掌握Matlab插值函数，处理插值应用问题插值函数，处理插值应用问题 了解分段线性插值的基本思想了解分段线性插值的基本思想 了解三次样条插值的提法和思路了解三次样条插值的提法和思路 掌握插值函数掌握插值函数 interp interp1 interp2 griddata插 值图像图

4、像插值插值放大放大 数码相机运用插值的方法可以创造出比传感器实际像素数码相机运用插值的方法可以创造出比传感器实际像素更多的图像，这种处理称为更多的图像，这种处理称为“数码变焦数码变焦”。106*40原始图像：原始图像：左边：左边：最近邻插值最近邻插值放大放大450%右边：右边：双三次插值双三次插值放大放大450%插值在图像三维重建中的应用插值 一维插值 二维插值q 最邻近点插值最邻近点插值q 分片线性插值分片线性插值q 分片双线性插值分片双线性插值q 分片双三次样条插值分片双三次样条插值q 最邻近点插值最邻近点插值q 分段分段线性插值线性插值q 三次样条插值三次样条插值q 埃尔米特埃尔米特插值

5、插值用用直线直线(线性线性)连接数据点列上相邻的两点连接数据点列上相邻的两点。比如在比如在两点两点xi-1,xi上线性插值函数上线性插值函数为为11111( ),0,1,2,.,iiiiiiiiiixxxxq xyyxxxxxxxin分段分段线性插值线性插值三次样条(spline)插值 是一个细的、可弯曲的木制或塑料条，在飞机或轮是一个细的、可弯曲的木制或塑料条，在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线船等的设计制造过程中为描绘出光滑的外形曲线( (放样放样) )所用的工具所用的工具 从物理上讲，样条满足插值点的约束，同时使势能从物理上讲，样条满足插值点的约束，同时使势能达到最小达

6、到最小 三次样条本质上是一段一段的三次多项式拼合而成三次样条本质上是一段一段的三次多项式拼合而成的曲线，在拼接处的曲线，在拼接处, ,不仅函数是连续的不仅函数是连续的, ,且一阶和二阶导且一阶和二阶导数也是连续的数也是连续的，19461946年年,Schoenberg,Schoenberg将样条引入数学将样条引入数学, ,即即所谓的所谓的样条函数样条函数q 什么是样条？什么是样条？( )( ) , S xf xa b三三次次样样则则条条称称为为在在上上的的插插值值函函数数q 三次样条插值函数三次样条插值函数(2)( ),( ),( ) , S xS xSxa b 都都在在区区间间上上连连续续1

7、(1)( ),iiS xxx 在在每每个个小小区区间间上上都都是是三三次次多多项项式式( ) , S xa b则则称称为为区区间间上上的的三三次次样样条条函函数数(3)( )S x而而三三次次样样条条函函数数满满足足(),0,1,iif xy in(),0,1,iiS xy in01, , naxxxba b为为区区间间的的一一个个分分割割( ) , :S xa b如如果果函函数数在在区区间间上上满满足足条条件件q 三次样条插值原理三次样条插值原理niyxSii, 1 , 0,)(在在n n个小区间构造个小区间构造S(x)S(x)，共有，共有n n个三次多项式，需确定个三次多项式，需确定4n4

8、n个参数个参数1, 1,)()(limniyxSxSiixxi1, 1),()(limnixSxSixxi1, 1,)()(lim niMxSxSiixxi在所有节点上在所有节点上n+1个方程个方程在除端点外的节点上在除端点外的节点上3(n-1)个方程个方程这样就得到这样就得到4n - 2个方程，为保证待定参数的唯一性，还差两个个方程，为保证待定参数的唯一性，还差两个方程。为此，常用的方法是对边界节点除函数值外附加要求，方程。为此，常用的方法是对边界节点除函数值外附加要求，这就是所谓的边界条件。根据实际问题的不同，三次样条插值这就是所谓的边界条件。根据实际问题的不同，三次样条插值常用到下列常用

9、到下列三类边界条件三类边界条件。 m边界条件：边界条件： 即给定端点处的一阶导数值。即给定端点处的一阶导数值。00() , ()nnS xyS xy M边界条件：边界条件： 即给定端点处的二阶导数即给定端点处的二阶导数值（值（自然边界条件）。00() , ()nnSxySxy 周期边界条件周期边界条件： 当当y = g(x)是以是以 b-a= x0 - xn 为周期的周期函数时，要求为周期的周期函数时，要求S(x)也是周期函数，故端点处要满足也是周期函数，故端点处要满足S(a + 0) = S(b - 0), S(a + 0) = S(b - 0)，此条件称为周期条件此条件称为周期条件。 Ma

10、tlab命令：命令： 用用Matlab实现分段线性插值不需要编制函数程序，实现分段线性插值不需要编制函数程序，Matlab中有中有现成的一维插值函数现成的一维插值函数interp1。y=interp1(x0,y0,x, method)method指定插值的方法，默认为线性插值。其值可为：指定插值的方法，默认为线性插值。其值可为： nearest 最近点等值方式插值最近点等值方式插值 linear 线性插值线性插值 spline 三次样条插值三次样条插值 cubic（新版本改为（新版本改为pchip） 三次三次Hermite插值插值或立方插值。或立方插值。所有的插值方法要求所有的插值方法要求x0

11、是单调的。是单调的。样条插值还可以这样：样条插值还可以这样： y=spline(x0,y0,x)q 分段线性插值与样条插值分段线性插值与样条插值一元插值函数 MATLAB中的一元插值函数为interp1( )， 它的功能是一维数据插值（表格查找）。该命令对数据点之间进行计算内插值，它出一元函数f(x)在中间点的数值，其中函数f(x)由所给数据决定。语法形式说明y=interp1(x,Y,xi)由已知点集(x,Y）插值计算xi上的函数值y=interp1(x,Y,xi)相当于x=1:length(Y)的interp(x,Y,xi)y=interp1(x,Y,xi,method)用指定插值方法计算

12、插值点xi上的函数值y=interp1(x,Y,xi,method,extrap)对xi中超出已知点集的插值点用指定插值方法计算函数值y=interp1(x,Y,xi,method,extrap,extrapval)用指定方法插值xi上的函数值，超出已知点集处函数值取extrapvaly=interp1(x,Y,xi,method,pp)用指定方法插值，但返回结果为分段多项式method方法描述nearest最邻近插值：插值点处函数值取与插值点最邻近的已知点的函数值liner分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，MATLAB中interp1的默认方法spline样条

13、插值：默认为三次样条插值。可用spline函数代替pchip三次Hermite多项式插值。可用pchip函数代替cubic同pchip，三次Hermite多项式插值 （1）nearest方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑； （2）spline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方法小，但当已知数据点不均匀分布时可能出现异常结果。 （3）cubic三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比spline方法略快，但占用内存最

14、多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。例 用interp1对sin函数进行分段线性插值。解：在MATLAB命令窗口中输入以下命令： x=0:2*pi; y=sin(x); xx=0:0.5:2*pi yy=interp1(x,y,xx); plot(x,y,s,xx,yy)例 ：给出数据表格如下：试求三次样条函数满足自然边界条件 解：在matlab工作区输入如下代码：x=28 28.7 29 30;y=4.3 3.9 4.1 3.0;x0=28:0.1:30;yy=spline(x,y,x0);plot(x,y,b*,x0,yy,r);结果如图所示：x2828.72930

15、4.33.94.13.0( )f x(28.7)(30)0SS例：已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。X035791112131415Y01.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6解：在matlab工作区输入如下代码：x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;y=spline(x0,y0,x);plot(x,y,r,x0,y0,b*);title(飞机下轮廓线的三次样条插值);axis e

16、qual;axis(0,15,0,3);例例 假设已知的数据点来自函数。假设已知的数据点来自函数。x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);%原始数据原始数据xdat=0:.02:1; ydat=(xdat.2-3*xdat+5).*exp(-5*xdat).*sin(xdat);%插值点坐标插值点坐标y1=interp1(x,y,xdat); %线性插值结果线性插值结果 y2=interp1(x,y,xdat, pchip); %三次三次Hermite插值结果插值结果y3=interp1(x,y,xdat, spline); %三次样条插值结果三次

17、样条插值结果y4=interp1(x,y,xdat, nearest); %最近点等值方式插值结果最近点等值方式插值结果plot(xdat,y1,y2,y3,y4,:,x,y,o,xdat,ydat)legend(线性插值线性插值, 三次三次Hermite插值插值, 三次样条插值三次样条插值, 最近点等值方式插值最近点等值方式插值, 原始数据原始数据, 插值点插值点)25( )(35)sinxf xxxex产生数据的命令为：产生数据的命令为：x=0:.12:1;y=(x.2-3*x+5).*exp(-5*x).*sin(x);试根据生成的数据进行插值处理，得出较平滑的曲线。试根据生成的数据进行

18、插值处理，得出较平滑的曲线。Matlab程序程序：例例 为了绘制乙醇水二元溶液的气液平衡曲线，查到为了绘制乙醇水二元溶液的气液平衡曲线，查到2424点实验点实验数据如表数据如表4-44-4所示。表中所示。表中x x和和y y分别表示液相和气相中乙醇摩尔分分别表示液相和气相中乙醇摩尔分数。但表中缺少低浓度点的数据，且数据分布不均匀。试用三数。但表中缺少低浓度点的数据，且数据分布不均匀。试用三次样条插值法对以下插值点插值并画出乙醇水二元溶液的气次样条插值法对以下插值点插值并画出乙醇水二元溶液的气液平衡曲线。液平衡曲线。 表表 乙醇水二元溶液的气液平衡曲线实验值乙醇水二元溶液的气液平衡曲线实验值x=

19、0, 0.018,0.054,0.078,0.111,0.147,0.172,0.203,0.210,0.241,. 0.297,0.320,0.381,0.391,0.423,0.425,0.497,0.507,0.567,0.685,. 0.767,0.838,0.882,0.894;y=0, 0.163,0.328,0.386,0.452,0.491,0.512,0.524,0.520,0.546,. 0.573,0.577,0.600,0.608,0.621,0.623,0.653,0.656,0.684,0.743,. 0.796,0.846,0.884,0.894;xk=0.010

20、,0.018,0.020,0.050,0.078,0.100,0.118,0.200,0.300,0.400,. 0.425,0.500,0.600,0.700,0.800,0.838,0.866,0.899;yk=spline(x,y,xk); %yk=csapi(x,y,xk);%pp=csape(x,y,second);%yk=ppval(pp,xk);plot(x,y,*,xk,yk,r);grid;legend(样本点,三次样条插值,Location,NorthWest)Matlab程序一程序一：这里：这里：x和和y分别表示样本点的液相和气相中乙醇摩尔分别表示样本点的液相和气相中乙醇

21、摩尔分数，分数，xk和和yk分别表示插值点的液相和气相中乙醇摩分别表示插值点的液相和气相中乙醇摩尔分数。尔分数。00.10.20.30.40.50.60.70.80.900.10.20.30.40.50.60.70.80.9 样 本 点三 次 样 条 插 值 已知二元函数已知二元函数g(x, y)在某矩形区域在某矩形区域R = a,b c,d上上互异节点互异节点(x i , y j) 的函数值的函数值 z ij ，如何求出在，如何求出在R 上任一点上任一点(x, y)处的函数值处的函数值g(x, y)的近似值。的近似值。 根据已知数据点根据已知数据点(x i , y j , z ij) ，构造

22、一个具有一定光滑性、，构造一个具有一定光滑性、简单易于计算的函数简单易于计算的函数z = f (x, y)（已知曲面）作为（已知曲面）作为z = g(x, y)（未知曲面）的近似，使曲面（未知曲面）的近似，使曲面z = f (x, y)通过已知数据点，即通过已知数据点，即 f (x i , y j ) = z ij 。然后计算点。然后计算点(x, y)的函数值的函数值f (x, y)作为作为g(x, y)的近似值。的近似值。二二维插值维插值数学描述数学描述2基本思想基本思想3网格节点插值网格节点插值 已知已知 m n 个节点个节点 ), 2, 1;,.,2, 1(),(njmizyxijji

23、其中其中jiyx ,互不相同，不妨设互不相同，不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点, ,即即再用再用),(yxf计算插值，即计算插值，即).,(*yxfz ), 1 ,0;, 1 ,0(),(njmizyxfijji q 最邻近点插值最邻近点插值网格节点插值法主要有以下几种形式：网格节点插值法主要有以下几种形式：q 分片线性插值分片线性插值q 分片双线性插值分片双线性插值q 分片双三次样条插值分片双三次样条插值4散乱节点插值散乱节点插值已知已知n n个节点个节点),.,2, 1(),(nizyxiii 其中

24、其中),(iiyx互不相同，互不相同， 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点, ,即即), 1 ,0(),(nizyxfiii 再用再用),(yxf计算插值，即计算插值，即).,(*yxfz 5二维插值的二维插值的Matlab实现实现q 网格节点插值网格节点插值要求要求x0,y0单调；单调；x，y可取可取为矩阵，或为矩阵，或x 取行向量，取行向量，y 取为列向取为列向量，量，x,y 的值分别不能超出的值分别不能超出x0,y0的范围。的范围。(1) z = interp2 ( x0, y0, z0, x, y, method)被插值点被插值点插值方法插

25、值方法插值插值节点节点被插值点的被插值点的函数值函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值splinespline 三次样条三次样条cubiccubic 双三次插值双三次插值缺省时缺省时 双线性插值双线性插值例例 假设由二元函数假设由二元函数222( , )(2 )xyxyzf x yxx e可以计算出一些较稀疏的网格数据，试根据这些数据对整个函可以计算出一些较稀疏的网格数据，试根据这些数据对整个函数曲面进行各种插值拟合，比较拟合结果，并进行误差分析。数曲面进行各种插值拟合，比较拟合结果，并进行误差分析。原始数据与原始面图原始数据与原

26、始面图x,y=meshgrid(-3:.6:3, -2:.4:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); %原始数据原始数据figure(1)surf(x,y,z) %原始数据面图原始数据面图axis(-3,3,-2,2,-0.7,1.5) %设置坐标轴范围设置坐标轴范围box ontitle(原始数据面图原始数据面图)xdat,ydat=meshgrid(-3:.2:3, -2:.2:2); %插值点坐标插值点坐标线性插值面图线性插值面图zdat1=interp2(x,y,z,xdat,ydat); %线性插值结果线性插值结果 figure(2)surf(xda

27、t,ydat,zdat1) %线性插值面图线性插值面图axis(-3,3,-2,2,-0.7,1.5) %设置坐标轴范围设置坐标轴范围box ontitle(线性插值面图线性插值面图)立方插值面图立方插值面图zdat2=interp2(x,y,z,xdat,ydat,cubic); %立方插值算法结果立方插值算法结果figure(3)surf(xdat,ydat,zdat2) %立方插值算法面图立方插值算法面图axis(-3,3,-2,2,-0.7,1.5) %设置坐标轴范围设置坐标轴范围box ontitle(立方插值算法面图立方插值算法面图)样条插值面图样条插值面图zdat3=interp

28、2(x,y,z,xdat,ydat,spline); %样条插值结果样条插值结果figure(4)surf(xdat,ydat,zdat3) %样条插值面图样条插值面图axis(-3,3,-2,2,-0.7,1.5) %设置坐标轴范围设置坐标轴范围box ontitle(样条插值面图样条插值面图)样条插值误差面图样条插值误差面图zdat=(xdat.2-2*xdat).*exp(-xdat.2-ydat.2-xdat.*ydat); figure(6)surf(xdat,ydat,abs(zdat-zdat3) %样条插值对应的误差面图样条插值对应的误差面图axis(-3,3,-2,2,0,0

29、.025) %设置坐标轴范围设置坐标轴范围box ontitle(样条插值算法对应的误差面图样条插值算法对应的误差面图)插值点处的函数值插值点处的函数值 z = ppval ( pp, x, y ) (2) pp = csape ( x0, y0 , z0, conds, valconds) 其中其中x0,y0,z0为已知的数据，为已知的数据，x0为为m维向量，维向量， y0为为n维向量，维向量，z0为为mn维矩阵。维矩阵。 x,y为一维数组表示插值点（为一维数组表示插值点（x为行向量，为行向量，y为为列向量。列向量。若不考虑边界条件，以上两条语句可用下面命令代替若不考虑边界条件，以上两条语句

30、可用下面命令代替z=csapi(x0,y0,z0,x,y);被插值点被插值点插值边界条件及取值，插值边界条件及取值，同一维情形同一维情形插值节点插值节点被插值点的被插值点的函数值函数值系数矩阵系数矩阵例例3-8 在一丘陵地带测量高程，在一丘陵地带测量高程，x和和y方向每隔方向每隔100米测一个点，米测一个点，得高程如下表，试拟合一曲面，确定合适的模型，并由此找出得高程如下表，试拟合一曲面，确定合适的模型，并由此找出最高点和该点的高程。最高点和该点的高程。高程x100200300400500y1006366976244784502006987126304784203006806745984124

31、00400662626552334310原始数据面图原始数据面图x=100:100:500;y=100:100:400;z=636 697 624 478 450698 712 630 478 420680 674 598 412 400662 626 552 334 310; %原始数据原始数据figure(1)surf(x,y,z)title(原始数据面图原始数据面图)100200300400500100200300400200400600800Y原 始 数 据 面 图XZ插值数据面图插值数据面图xdat,ydat=ndgrid(100:10:500,100:10:400); %插值点网格

32、数据插值点网格数据%pp=csape(x,y,z); %zdat=fnval(pp,100:10:500,100:10:400); %经插值运算得到的插值点处经插值运算得到的插值点处函数值函数值zdat=csapi(x,y,z,100:10:500,100:10:400); %经插值运算得到的插值经插值运算得到的插值点处函数值点处函数值figure(2)surf(xdat,ydat,zdat) %插值面图插值面图title(插值数据面图插值数据面图)100150200250300350400450500100150200250300350400200300400500600700800插 值

33、数 据 面 图最高点的高程求解最高点的高程求解x_temp=100:500; %以以1为步长产生数据为步长产生数据y_temp=100:400;z_temp=csapi(x,y,z,100:500,100:400);z_max=max(z_temp(:);x_max,y_max=find(z_temp=z_max);xyz_max= x_temp(x_max),y_temp(y_max), z_max %最高点和该点的高程最高点和该点的高程z = griddata ( x0, y0, z0, x, y, method )q 散乱节点插值散乱节点插值被插值被插值点点插值方插值方法法插值节插值节点

34、点被插值点的被插值点的函数值函数值nearest 最邻近插值最邻近插值linear 双线性插值双线性插值cubic 双三次插值双三次插值v4 Matlab提供的插值方法提供的插值方法缺省时缺省时 双线性插值双线性插值其中其中x0,y0,z0为已知的数据，这里并不要求是网格型的，可以为已知的数据，这里并不要求是网格型的，可以是任意分布的，均由向量给出。是任意分布的，均由向量给出。 x,y为插值点，可以是单个为插值点，可以是单个点，可以是向量或网格型矩阵，得出的应该维数和点，可以是向量或网格型矩阵，得出的应该维数和x，y一致，一致，表示插值的结果。表示插值的结果。 method 表示插值算法，可以

35、为表示插值算法，可以为linear，cubic，nearest，v4其中其中v4是是Matlab4.0版本中提供的插值算法，公认效果较好，版本中提供的插值算法，公认效果较好，但没有一个正式的名称，所以用但没有一个正式的名称，所以用v4表明该算法。表明该算法。例例 已经有了已经有了 x，y，z 的坐标值，如的坐标值，如 x=1 2 3 4; y=2 3 4 5; z=5 6 7 4。由这些坐标点通过插值方法拟合出三维曲面由这些坐标点通过插值方法拟合出三维曲面 。Matlab程序程序x=1 2 3 4; y=2 3 4 5; z=5 6 7 4;x1,y1=meshgrid(1:0.1:4,2:0

36、.1:5);z1=griddata(x,y,z,x1,y1,v4);surf(x1,y1,z1)12342345051015XYZ题例题例 粗糙山顶曲面的平滑处理粗糙山顶曲面的平滑处理( (散乱情形散乱情形) )rand(seed,0)x = rand(100,1)*4-2; y = rand(100,1)*4-2;z = x.*exp(-x.2-y.2);plot3(x,y,z,o);hold onti = -2:.25:2; XI,YI = meshgrid(ti,ti);ZI = griddata(x,y,z,XI,YI);mesh(XI,YI,ZI);题例题例 墨西哥草帽的平滑处理墨西

37、哥草帽的平滑处理( (散乱情形散乱情形) ) x = rand(100,1)*16 - 8;y = rand(100,1)*16 - 8;r = sqrt(x.2 + y.2) + eps;z = sin(r)./r;plot3(x,y,z,.,MarkerSize,15)hold onxlin = linspace(min(x),max(x),33);ylin = linspace(min(y),max(y),33);X,Y = meshgrid(xlin,ylin);Z = griddata(x,y,z,X,Y,cubic);mesh(X,Y,Z); axis tight; 南半球气旋变化

38、的可视图形山区地貌的可视化图形 Matlab命令命令：1. z=interp3(x0,y0,z0,x1,y1, method) 2. z=interpn(x1,x2,xn,z0,xd1,xd2,xdn, method) 3. z=csapi(x1,x2,xn,z0,xd1,xd2,xdn);4. v=griddata3(x0,y0,z0,v0,x,y,z, method);5. v=griddata3(X,Y,XI method); 具体用法请参考具体用法请参考Matlab帮助帮助高高维插值的维插值的Matlab实现实现 内容：拟合是离散函数逼近的重要方法，利用它可通过内容：拟合是离散函数逼近

39、的重要方法，利用它可通过函数在有限个点处的取值状况，拟合出近似替代函数，函数在有限个点处的取值状况，拟合出近似替代函数，进而估算出函数在其他点处的近似值。进而估算出函数在其他点处的近似值。 了解基于最小二乘法则拟合的基本思想了解基于最小二乘法则拟合的基本思想 掌握拟合函数掌握拟合函数 polyfit lsqcurvefit curvefit 掌握掌握cftool曲线拟合工具箱曲线拟合工具箱(多目标函数多法则多目标函数多法则)拟 合第九章第九章 第十第十节节最小二乘法最小二乘法 多项式曲线拟合 对给定的试验数据点(xi,yi)(i=1,2,N)，可以构造m次多项式： 由曲线拟合的定义，应该使得下

40、式取极小值 通过简单运算可以得出系数是下面线性方程组的解 )()(10NmxaxaaxPmm20)(iNimjjijyxa mmmmmmmbbbaaaccccccccc10102112110 其中NikiikNikikmkxybmkxc11)2 , 1 , 0( ,)2 , 1 , 0( ,最小二乘拟合 最小二乘法拟合在科学实验的统计方法中经常使用。它的具体操作过程是从一组实验数据(xi,yi)中拟合出函数关系y=f(x)，拟合的标准是使(f(xi)-yi)的平方取极小值。 在工程实践和科学计算中，用某种在工程实践和科学计算中，用某种经验函数经验函数解析式解析式y=f(x)来近似来近似刻画采集

41、数据刻画采集数据(x,y) 之间的关系的方法就叫拟合，所谓之间的关系的方法就叫拟合，所谓“拟合拟合”有有 “最贴近最贴近”之意之意 。 与插值不同，拟合的主要目与插值不同，拟合的主要目标是要离散点尽量靠近拟合函标是要离散点尽量靠近拟合函数。一般过程是，我们首先根数。一般过程是，我们首先根据采样点的散点分布图，大致据采样点的散点分布图，大致推测推测x与与y之间的经验函数形式之间的经验函数形式(比如多项式、指数函数等比如多项式、指数函数等)， 然后依据某种法则然后依据某种法则(比如最常用的最小二乘法则比如最常用的最小二乘法则)，确定出的经验，确定出的经验函数解析式中的待定参数。其中函数解析式中的待

42、定参数。其中经验函数经验函数和和拟合法则拟合法则是拟合的两个是拟合的两个关键要素！关键要素！ 经验函数形式：经验函数形式：已经拟定为多项式函数：已经拟定为多项式函数：y= at2 +bt+ c剩下的工作是确定拟合原则：剩下的工作是确定拟合原则：可选的法则很多，其中最常用的是最小二乘法则可选的法则很多，其中最常用的是最小二乘法则(method of Least Squares)，即，即各点残差平方和最小各点残差平方和最小用用Matlab进行进行数据拟合数据拟合1. 多项式曲线拟合多项式曲线拟合: polyfit.y0=polyval(p,x0)p=polyfit(x,y,m)其中其中, x, y

43、为已知数据点向量为已知数据点向量, 分别表示横分别表示横,纵坐纵坐标标, m为拟合多项式的次数为拟合多项式的次数, 结果返回结果返回m次拟合次拟合多项式系数多项式系数, 从高次到低次存放在向量从高次到低次存放在向量p中中.可求得多项式在可求得多项式在x0处的值处的值y0.用四次多项式拟合下表中所列的数据点。解：x=1:0.2:1.8; y=1 1.09541.18321.26491.3416; p=polyfit(x,y,4)p = -0.0104 0.0854 -0.3121 0.9086 0.3285即拟合的多项式为p= -0.0104 + 0.0854x -0.3121x2+ 0.908

44、6 x3+ 0.3285x4。x11.21.41.61.8 y11.09541.18321.26491.3416例例1 已知观测数据点如表所示已知观测数据点如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分别用分别用3次和次和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.

45、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: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(

46、x,y,6)例例2 用切削机床进行金属品加工时用切削机床进行金属品加工时, 为了适当地调整为了适当地调整机床机床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量刀在一定的时间测量刀具的厚度具的厚度, 得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm解解: 描出散点图描出散点图, 在命令窗口输入在命令窗口输入:t=

47、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 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,*)a = -0.3012 29.3804hold onplot(t, y1), hold offa=polyfit(t,y,

48、1)y1=-0.3012*t+29.3804例例2 用切削机床进行金属品加工时用切削机床进行金属品加工时, 为了适当地调整为了适当地调整机床机床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量刀在一定的时间测量刀具的厚度具的厚度, 得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm拟合曲线为拟合曲线为: y=-0.

49、3012t+29.3804例例3 一个一个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.953 10 3. 0362000 10 31.517 10 已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 2.896

50、2375 10 31.219 10 已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 解解 选取选取指数函数指数函数作拟合时作拟合时, 在拟合前需作变量代换在拟合前需作变量代换, 化为化为 k1, k2 的线性函数的线性函数.于是于是,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.

51、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*x)plot(x,w,*,x,y1,r-)已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 拟合曲线为拟合曲线为:3 -4

52、94.52094.0275 10 e0211ln,ln, zakak01 zaa0211-494.5209,ln8.3009, akak3124.0275 10 ,494.5209 kk令令则则求得求得于是于是在实际应用中常见的拟合曲线有在实际应用中常见的拟合曲线有:01ya xa直线直线101 nnnya xa xa多项式多项式一般一般 n=2, 3, 不宜过高不宜过高.01ayax双曲线双曲线(一支一支) bxyae指数曲线指数曲线2. 非线性曲线拟合非线性曲线拟合: nlinfit.功能功能:x=lsqcurvefit(fun, x0, xdata, ydata)x, resnorm=l

53、sqcurvefit(fun, x0, xdata, ydata)根据给定的数据根据给定的数据 xdata, ydata (对应点的横对应点的横, 纵坐标纵坐标), 按函数文件按函数文件 fun 给定的函数给定的函数, 以以x0为初值作为初值作最小二乘拟合最小二乘拟合, 返回函数返回函数 fun中的中的系数向量系数向量x和残差的平方和和残差的平方和resnorm.CFTool曲线拟合工具箱简介 在在command window中键入指令中键入指令cftool即可启动即可启动曲线拟合工具箱。在该集成环境里面，可以实现多种曲线拟合工具箱。在该集成环境里面，可以实现多种经验函数，多种法则的曲线拟合，

54、实时绘制图像并进经验函数，多种法则的曲线拟合，实时绘制图像并进行误差分析行误差分析。 需要注意的是：在进入需要注意的是：在进入Curve Fitting Toolbox环境进行曲线拟合之前，需要预先在环境进行曲线拟合之前，需要预先在workspace输入输入或载入供拟合的数据源或载入供拟合的数据源CFTool曲线拟合工具箱简介下面还是以引例的采样数据为例，进行演示：下面还是以引例的采样数据为例，进行演示：t=1:16;y=4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6; cftool导入数据导入数据

55、绘制散点图绘制散点图进行拟合进行拟合拟合方法拟合方法结果和误差分析结果和误差分析 这里可供选择的这里可供选择的拟合类型拟合类型和和可选参可选参数数比较多比较多,包括多项式函数包括多项式函数,指数函数指数函数,幂函数等幂函数等,如何确定如何确定最优的方案最优的方案？CFTool-拟合效果评价指标SSE - The sum of squares due to error. This statistic measures the deviation of the responses from the fitted values of the responses. A value closer to

56、0 indicates a better fit. R-square - The coefficient of multiple determination. This statistic measures how successful the fit is in explaining the variation of the data. A value closer to 1 indicates a better fit. Adjusted R-square - The degree of freedom adjusted R-square. It is generally the best

57、 indicator of the fit quality when you add additional coefficients to your model. A value closer to 1 indicates a better fit. RMSE - The root mean squared error. A value closer to 0 indicates a better fit. CFTool-选择Analysis分析数据分析数据分析数据 这里可以计算拟合函数在采样点的这里可以计算拟合函数在采样点的一阶二阶导函数值，并绘制相应的函一阶二阶导函数值，并绘制相应的函数图

58、像数图像CFTool-拟合薄膜渗透数学模型第一步：导入数据第一步：导入数据load bmstl.matcftool第二步：拟合数据第二步：拟合数据syms a b c d tf = a*exp(b*t) + c*exp(d*t)cbvst=subs(f,a,b,c,d,0.008839,-0.0001336,-0.004818,-0.001422);f = inline(cbvst)figuret=100:1:1000;plot(t,f(t),b-)CFTool-拟合薄膜渗透数学模型第三步：分析数据第三步：分析数据延伸阅读 样条插值实验.doc Matlab进行一般插值和三次样条插值比较.do

59、c 曲线拟合cftool工具箱.ppt Matlab的曲线拟合工具箱CFtool使用简介.doc曲 线 拟 合 和 插 值 函 数polyfit(x, y, n)对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合interp1(x, y, xo)1维线性插值interp1(x, y, xo, spline )1维3次样条插值interp1(x, y, xo, cubic )1维3次插值interp2(x, y, Z, xi, yi)2维线性插值interp2(x, y, Z, xi, yi, cubic )2维3次插值interp2(x, y, Z, xi, yi, nearest )2维

60、最近邻插值神经网络潘多拉的魔盒支持向量机 “XX青年”的分类 神经网络 支持向量机q 分类分类q 回归回归q 分类分类q 回归回归 一般而言, ANN与经典计算方法相比并非优越, 只有当常规方法解决不了或效果不佳时ANN方法才能显示出其优越性。尤其对问题的机理不甚了解或不能用数学模型表示的系统,如故障诊断、特征提取和预测等问题,ANN往往是最有利的工具。另一方面, ANN对处理大量原始数据而不能用规则或公式描述的问题, 表现出极大的灵活性和自适应性。人工神经网络人工神经网络 (Artificial Neural Netwroks -ANN)潘多拉的魔盒训练训练样本样本训练训练目标目标待测待测样