数据插值与拟合讲课文档
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第九页，共25页。
（3）三次样条插值
第十页，共25页。
2、曲线拟合的最小二乘法
给定平面上的点
进行曲线拟合有多种方法，最小二乘法是解决曲线拟合最常 用的一种方法 最小二乘法的原理是求f(x),使
达到最小 简单地说，最小二乘法准则就是使所有散点到曲线的距离平 方和最小
第十一页，共25页。
线性最小二乘法
第十六页，共25页。
例2 气旋变化情况的可视化
下表是气象学家测量得到的气象资料，它们分别表示在南半球地时 按不同纬度。不同月份的平均气旋数字．根据这些数据，绘制出气 旋分布曲面图形
第十七页，共25页。
y=5:10:85;x=1:12;
[x,y]=meshgrid(x,y); plot(x,y,'*'); pause
第十三页，共25页。
y iinetr1(p x,y,x,i'met')hod
表示采用的插值方法 MATLAB提供的插值方法有几种
：分段线性插值
' p c h ip ' ：三次Hermite插值（立方插值）
：三次分段样条插值
'nearest ' ：最近点等值方式
缺省时表示线性插值
第十四页，共25页。
第三页，共25页。
一、实例及其模型
1、船在该海域会搁浅吗 在某海域测得一些点(x，y)处的水深z（单位：英尺）由下表给 出，水深数据是在低潮时测得的．船的吃水深度为5英尺，问 在矩形区域(75,200)*(-50，150）里的哪些地方船要避免进人．
第四页，共25页。
分析 由于测量点是散乱分布的，先在平面上作出测量点的分布 图，再利用二维插值方法补充一些点的水深，然后作出海 底曲面图和等高线图，并求出水深小于5的海域范围．
p=lsqcurvefit(‘Fun’,p0,xdata,ydata)
其中Fun表示函数Fun(p,data)的M函数文件，p0表示函数的 初值.。
注：若要求解点x处的函数值可用程序f=Fun(p,x)计算
第二十四页，共25页。
(1)函数原型m文件 function y=fname(a,t) y=a(1)*exp(-a(2)*t);
数据插值与拟合
第一页，共25页。
（优选）第九讲数据插值 与拟合
第二页，共25页。
插值则要求函数在每个观测点处一定要满足
插值函数一般是已知函数的线性组合或者称为加权平 均．插值在工程实践和科学实验中有着非常广泛而又十分 重要的应用，例如，信息技术中的图像重建、图像放大中 为避免图像的扭曲失真的插值补点、建筑工程的外观设计。 化学工程实验数据与模型的分析、天文观测数据、地理信 息数据的处理如(天气预报）以及社会经济现象的统计分 析等等．
a=polyfit(xdata,ydata,n) 其中n表示多项式的最高阶数，xdata，ydata为将要拟合的数 据，它是用数组的方式输入
输出参数a为拟合多项式 的系数
注：多项式在x处的值y可用下面程序计算．
y=polyval(a,x)
第二十页，共25页。
解：Matlab程序
T=[19.1 25.0 30.1 36.0 40.0 45.1 50.0];
x,y,z是已知样本点的坐标，可以是任意分布的。
X0,y0是期望的插值位置，即被插值节点, 可以是单点，向量或者 网格型矩阵
插值方法，除了上面的 方法外，还有一个是4.0版本提供的 一个插值方法，选项为’v4’
第十九页，共25页。
四、曲线拟合的matlab实现
1、已知函数原型的
（1）多项式拟合
假设已知函数原型为 Matlab提供的拟合函数为
拟合函数可由一些简单的“基函数”（例如幂函数，三角
函数等等）
来线性表示
现在要确定系数
使
达到极小为此
第十二页，共25页。
三、插值的matlab实现
1、一维插值
MATLAB中的插值函数为interp1，其调用格式为
其中x，y为插值点，yi为在被插值点xi处的插值结果，x,
y为向量。
注意：所有的插值方法都要求x是单调的，并且xi不能 够超过x的范围。
第七页，共25页。
二、 插值与拟合
1、插值方法
（1）分段线性插值
分段线性插值的提法如下：
第八页，共25页。
（2）分段三次埃尔米特插值 在插值问题中，如果除了插值节点的函数值给定外，还要 求在节点的导数值为给定值，即插值问题变为
相当于在每一小段上应满足四个条件（方程），可以确定四 个待定参数．三次多项式正好有四个系数，所以可以考虑用 三次多项式函数作为插值函数，这就是分段三次埃尔米特插 值，它与分段线性插值一起都称为分段多项式插值
surf(x,y,z) pause
[xi,yi]=meshgrid(1:12,5:1:85); zi=interp2(x,y,z,xi,yi,'spline' ); figure mesh(xi,yi,zi)
xlabel('月份'), ylabel('纬度'), zlabel('气旋'), axis([0 12 0 90 0 50]) title('南半球气旋可视化图形')
通过求解线性方程可得待定系数，一般方法： X=[…] %已知数据x的列向量
Y=[…] %已知数据y的列向量
A=[f1(X),f2(X),…,fm(X)] %系数矩阵，fm（）为基函数 c=A\y
第二十二页，共25页。
解：matlab程序
X=[0 0.1 0.2 0.3 0.4 0.5 0.6]'
Y=[2 2.20254 2.40715 2.61592 2.83096 3.05448 3.28876]' A=[ones(size(X)),exp(X),exp(-X)]; c1=A\Y;
R=[76.30 77.80 79.25 80.80 82.35 83.90 85.10];
PR=polyfit(T,R,1);
t=10:60; r=polyval(PR,t);
plot(T,R,'*',t,r)
第二十一页，共25页。
（2）一般函数线性组合的曲线拟合
假设已知函数原型为 f( x ) c 00 ( x ) c 11 ( x ) c m m ( x )
问题分析
（1）首先将这些离散数据分布在直角坐标系下，由此可发现 浓度与时间之间呈现什么规律．这种数据分布在直角坐标 系下的图形被称为散点图；
（2）根据散点图，判段它接近于哪类函数曲线，
即确定函数形式 （3）函数形式确定以后，关键是要确定函数中含有的
待定参数。
最常用的确定待定系数的方法是，曲线拟合的最小二乘法
N维插值函数interpN（） 其中N可以为2，3，……，如N＝2为二维插值，调用格式为
其中 x，y，z为插值节点，zi为被插值点(xi，yi)处的插值结果 且， xi， yi为被插值节点构成的新的网格数据 ‘methods’代表的意思和可选择的插值方法和前面一样
注意：所有的插值方法都要求x和y是单调的网格，x和 y可以是等距的 也可以是不等距的
z=[2.4,1.6,2.4,3.2,1.0,0.5,0.4,0.2,0.5,0.8,2.4,3.6; 18.7 21.4 16.2 9.2 2.8 1.7 1.4 2.4 5.8 9.2 10.3 16;
20.8 18.5 18.2 16.6 12.9 10.1 8.3 11.2 12.5 21.1 23.9 25.5; 22.1 20.1 20.5 25.1 29.2 32.6 33.0 31.0 28.6 32.0 28.1 25.6; 37.3 28.8 27.8 37.2 40.3 41.7 46.2 39.9 35.9 40.3 38.2 43.4; 48.2 36.6 35.5 40 37.6 35.4 35 34.7 35.7 39.5 4在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一 组数据如表
表中的数据反映了浓度随时间变化的函数关系，它是一种离散关 系若需要推断20，40分钟时的浓度值，能否用一个显函数y=f(t) 来拟合表中的离散数据，然后再计算浓度值f(20), f（40）?
第六页，共25页。
注：因为后面可能要用到计算函数在一些点上的值，因此写 函数原型表达式时，记得用点运算
第二十五页，共25页。
C=c1' x=0:0.05:1;
y=C(1)+C(2)*exp(x)+C(3)*exp(-x);
plot(X,Y,'*',x,y)
第二十三页，共25页。
（3）一般的曲线拟合
假设已知函数原型是一般的函数，可以是多项式，可以是线 性，也可以是非线性的，一般情况下用这个来求解非线性情 况
Matlab在优化工具箱中提供的求解一般的曲线拟合函数 lsqcurvefit(),其调用格式如下
25.6 24.2 25.5 24.6 21.1 22.2 20.2 21.2 22.6 28.5 25.3 24.3;
5.3 5.3 5.4 4.9 4.9 7.1 5.3 7.3 7 8.6 6.3 6.6;
0.3,0,0,0.3,0,0,0.1,0.2,0.3,0,0.1,0.3]; figure
第十八页，共25页。
（2）、一般二维分布的数据插值
在实际应用问题中，大部分的数据以实测的多组（xi,yi,zi） 给出，所以不能直接使用interp2（）函数。
Matlab中提供了另一个函数griddata( )，用来专门解决 这类问题。其调用格式如下
Z=griddata(x,y,z,x0,y0,’method’)
例1 在一 天24小时内，从零点开始每间隔2小时测得的环境温度数据
分别为
12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，
推测中午（即13点）时的温度．
x=0：2：24； y=[12 9 9 10 18 24 28 27 25 20 18 15 13]；
x1＝13 ；
y1＝interp1（x，y，x1，‘spline’）
若要得到一天24小时的温度曲线
x=0：2：24； y=[12 9 9 10 18 24 28 27 25 20 18 15 13] xi＝0：1／3600：24； yi=interp1(x,y,xi,’spline’ ); plot(x, y, ’o’, xi, yi)
第十五页，共25页。
2、高维插值 (1) 网格数据插值问题
（3）三次样条插值
第十页，共25页。
2、曲线拟合的最小二乘法
给定平面上的点
进行曲线拟合有多种方法，最小二乘法是解决曲线拟合最常 用的一种方法 最小二乘法的原理是求f(x),使
达到最小 简单地说，最小二乘法准则就是使所有散点到曲线的距离平 方和最小
第十一页，共25页。
线性最小二乘法
第十六页，共25页。
例2 气旋变化情况的可视化
下表是气象学家测量得到的气象资料，它们分别表示在南半球地时 按不同纬度。不同月份的平均气旋数字．根据这些数据，绘制出气 旋分布曲面图形
第十七页，共25页。
y=5:10:85;x=1:12;
[x,y]=meshgrid(x,y); plot(x,y,'*'); pause
第十三页，共25页。
y iinetr1(p x,y,x,i'met')hod
表示采用的插值方法 MATLAB提供的插值方法有几种
：分段线性插值
' p c h ip ' ：三次Hermite插值（立方插值）
：三次分段样条插值
'nearest ' ：最近点等值方式
缺省时表示线性插值
第十四页，共25页。
第三页，共25页。
一、实例及其模型
1、船在该海域会搁浅吗 在某海域测得一些点(x，y)处的水深z（单位：英尺）由下表给 出，水深数据是在低潮时测得的．船的吃水深度为5英尺，问 在矩形区域(75,200)*(-50，150）里的哪些地方船要避免进人．
第四页，共25页。
分析 由于测量点是散乱分布的，先在平面上作出测量点的分布 图，再利用二维插值方法补充一些点的水深，然后作出海 底曲面图和等高线图，并求出水深小于5的海域范围．
p=lsqcurvefit(‘Fun’,p0,xdata,ydata)
其中Fun表示函数Fun(p,data)的M函数文件，p0表示函数的 初值.。
注：若要求解点x处的函数值可用程序f=Fun(p,x)计算
第二十四页，共25页。
(1)函数原型m文件 function y=fname(a,t) y=a(1)*exp(-a(2)*t);
数据插值与拟合
第一页，共25页。
（优选）第九讲数据插值 与拟合
第二页，共25页。
插值则要求函数在每个观测点处一定要满足
插值函数一般是已知函数的线性组合或者称为加权平 均．插值在工程实践和科学实验中有着非常广泛而又十分 重要的应用，例如，信息技术中的图像重建、图像放大中 为避免图像的扭曲失真的插值补点、建筑工程的外观设计。 化学工程实验数据与模型的分析、天文观测数据、地理信 息数据的处理如(天气预报）以及社会经济现象的统计分 析等等．
a=polyfit(xdata,ydata,n) 其中n表示多项式的最高阶数，xdata，ydata为将要拟合的数 据，它是用数组的方式输入
输出参数a为拟合多项式 的系数
注：多项式在x处的值y可用下面程序计算．
y=polyval(a,x)
第二十页，共25页。
解：Matlab程序
T=[19.1 25.0 30.1 36.0 40.0 45.1 50.0];
x,y,z是已知样本点的坐标，可以是任意分布的。
X0,y0是期望的插值位置，即被插值节点, 可以是单点，向量或者 网格型矩阵
插值方法，除了上面的 方法外，还有一个是4.0版本提供的 一个插值方法，选项为’v4’
第十九页，共25页。
四、曲线拟合的matlab实现
1、已知函数原型的
（1）多项式拟合
假设已知函数原型为 Matlab提供的拟合函数为
拟合函数可由一些简单的“基函数”（例如幂函数，三角
函数等等）
来线性表示
现在要确定系数
使
达到极小为此
第十二页，共25页。
三、插值的matlab实现
1、一维插值
MATLAB中的插值函数为interp1，其调用格式为
其中x，y为插值点，yi为在被插值点xi处的插值结果，x,
y为向量。
注意：所有的插值方法都要求x是单调的，并且xi不能 够超过x的范围。
第七页，共25页。
二、 插值与拟合
1、插值方法
（1）分段线性插值
分段线性插值的提法如下：
第八页，共25页。
（2）分段三次埃尔米特插值 在插值问题中，如果除了插值节点的函数值给定外，还要 求在节点的导数值为给定值，即插值问题变为
相当于在每一小段上应满足四个条件（方程），可以确定四 个待定参数．三次多项式正好有四个系数，所以可以考虑用 三次多项式函数作为插值函数，这就是分段三次埃尔米特插 值，它与分段线性插值一起都称为分段多项式插值
surf(x,y,z) pause
[xi,yi]=meshgrid(1:12,5:1:85); zi=interp2(x,y,z,xi,yi,'spline' ); figure mesh(xi,yi,zi)
xlabel('月份'), ylabel('纬度'), zlabel('气旋'), axis([0 12 0 90 0 50]) title('南半球气旋可视化图形')
通过求解线性方程可得待定系数，一般方法： X=[…] %已知数据x的列向量
Y=[…] %已知数据y的列向量
A=[f1(X),f2(X),…,fm(X)] %系数矩阵，fm（）为基函数 c=A\y
第二十二页，共25页。
解：matlab程序
X=[0 0.1 0.2 0.3 0.4 0.5 0.6]'
Y=[2 2.20254 2.40715 2.61592 2.83096 3.05448 3.28876]' A=[ones(size(X)),exp(X),exp(-X)]; c1=A\Y;
R=[76.30 77.80 79.25 80.80 82.35 83.90 85.10];
PR=polyfit(T,R,1);
t=10:60; r=polyval(PR,t);
plot(T,R,'*',t,r)
第二十一页，共25页。
（2）一般函数线性组合的曲线拟合
假设已知函数原型为 f( x ) c 00 ( x ) c 11 ( x ) c m m ( x )
问题分析
（1）首先将这些离散数据分布在直角坐标系下，由此可发现 浓度与时间之间呈现什么规律．这种数据分布在直角坐标 系下的图形被称为散点图；
（2）根据散点图，判段它接近于哪类函数曲线，
即确定函数形式 （3）函数形式确定以后，关键是要确定函数中含有的
待定参数。
最常用的确定待定系数的方法是，曲线拟合的最小二乘法
N维插值函数interpN（） 其中N可以为2，3，……，如N＝2为二维插值，调用格式为
其中 x，y，z为插值节点，zi为被插值点(xi，yi)处的插值结果 且， xi， yi为被插值节点构成的新的网格数据 ‘methods’代表的意思和可选择的插值方法和前面一样
注意：所有的插值方法都要求x和y是单调的网格，x和 y可以是等距的 也可以是不等距的
z=[2.4,1.6,2.4,3.2,1.0,0.5,0.4,0.2,0.5,0.8,2.4,3.6; 18.7 21.4 16.2 9.2 2.8 1.7 1.4 2.4 5.8 9.2 10.3 16;
20.8 18.5 18.2 16.6 12.9 10.1 8.3 11.2 12.5 21.1 23.9 25.5; 22.1 20.1 20.5 25.1 29.2 32.6 33.0 31.0 28.6 32.0 28.1 25.6; 37.3 28.8 27.8 37.2 40.3 41.7 46.2 39.9 35.9 40.3 38.2 43.4; 48.2 36.6 35.5 40 37.6 35.4 35 34.7 35.7 39.5 4在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一 组数据如表
表中的数据反映了浓度随时间变化的函数关系，它是一种离散关 系若需要推断20，40分钟时的浓度值，能否用一个显函数y=f(t) 来拟合表中的离散数据，然后再计算浓度值f(20), f（40）?
第六页，共25页。
注：因为后面可能要用到计算函数在一些点上的值，因此写 函数原型表达式时，记得用点运算
第二十五页，共25页。
C=c1' x=0:0.05:1;
y=C(1)+C(2)*exp(x)+C(3)*exp(-x);
plot(X,Y,'*',x,y)
第二十三页，共25页。
（3）一般的曲线拟合
假设已知函数原型是一般的函数，可以是多项式，可以是线 性，也可以是非线性的，一般情况下用这个来求解非线性情 况
Matlab在优化工具箱中提供的求解一般的曲线拟合函数 lsqcurvefit(),其调用格式如下
25.6 24.2 25.5 24.6 21.1 22.2 20.2 21.2 22.6 28.5 25.3 24.3;
5.3 5.3 5.4 4.9 4.9 7.1 5.3 7.3 7 8.6 6.3 6.6;
0.3,0,0,0.3,0,0,0.1,0.2,0.3,0,0.1,0.3]; figure
第十八页，共25页。
（2）、一般二维分布的数据插值
在实际应用问题中，大部分的数据以实测的多组（xi,yi,zi） 给出，所以不能直接使用interp2（）函数。
Matlab中提供了另一个函数griddata( )，用来专门解决 这类问题。其调用格式如下
Z=griddata(x,y,z,x0,y0,’method’)
例1 在一 天24小时内，从零点开始每间隔2小时测得的环境温度数据
分别为
12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，
推测中午（即13点）时的温度．
x=0：2：24； y=[12 9 9 10 18 24 28 27 25 20 18 15 13]；
x1＝13 ；
y1＝interp1（x，y，x1，‘spline’）
若要得到一天24小时的温度曲线
x=0：2：24； y=[12 9 9 10 18 24 28 27 25 20 18 15 13] xi＝0：1／3600：24； yi=interp1(x,y,xi,’spline’ ); plot(x, y, ’o’, xi, yi)
第十五页，共25页。
2、高维插值 (1) 网格数据插值问题