matlab插值曲线拟合
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
三、插值、曲线拟合
插值就是已知一组离散的数据点集，在集合内部某两 个点之间预测函数值的方法。
插值法是实用的数值方法，是函数逼近的重要方法。在 生产和科学实验中，自变量x与因变量y的函数y = f(x)的 关系式有时不能直接写出表达式，而只能得到函数在若干 个点的函数值或导数值。当要求知道观测点之外的函数值 时，需要估计函数值在该点的值。
x=-1:0.1:1; y=1./(1+9*x.^2); xi=-1:0.1:1; yi=interp1(x,y,xi); plot(x,y,'r-',xi,yi,'*')
例2. 在普通V带设计中，带轮的包角α与包角系数ka之间的关系如 表所示。求α=133.5°时的包角系数ka。
包角与包角系数
包角 ( °) 包角系 数 包角 ( °) 90 100 110 120 125 130 135 140
解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之 间, 所以对上述数据用三次样条插值加细为步长为1的数据: x=[18:2:30] y=[9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150] xi=[18:1:30] yi=interp1(x,y,xi,'spline') A=[xi;yi] A= 18.0000 25.0000 9.9618 9.9352
y01 = 0.639700000000000 y02 = 0.693100000000000 y03 = 0.641818996851421 y04 = 0.641831200000000
例4 对 y
性插值和三次样条插值, 用m=21个插值点作图,比较结果. 解 在命令窗口输入:
n=11, m=21 x=-5:10/(m-1):5 y=1./(1+x.^2) z=0*x x0=-5:10/(n-1):5 y0=1./(1+x0.^2) y1=interp1(x0,y0,x) y2=interp1(x0,y0,x,'spline') [x' y' y1' y2'] plot(x,z,'r',x,y,'k:',x,y1,'b',x,y2,'g') gtext('Piece.linear.'),gtext('Spline'),gtext('y=1/(1+x^2)')
Inteቤተ መጻሕፍቲ ባይዱp2
zi=interp2(z,xi,yi)
zi=interp2(x,y,z,xi,yi,method) 用指定的算法method计算二维插值。 Linear双线性插值(默认方式)， Nearest最邻近插值，spline三次样条 插值，cubic双三次插值。
函数
格式
功能 一维Fourier插值。适用于周期函数生 成数据的插值。X为抽样序列，n为 需要计算的等距点数，y为n点的插 值计算结果。 沿着指定的方向dim进行计算。 三次样条数据插值。用三次样条插 值计算出由向量x与y确定的一元函 数在点xx处的值yy。 返回由向量x与y确定的分段样条多 项式的系数矩阵pp。
y=interpft(x,n) Interpft y=interpft(x,n,dim)
yy=spline(x,y,xx)
spline
pp=spline(x,y)
注意： （1）只对已知数据点集内部的点进行的插值运算称为内插， 可比较准确的估测插值点上的函数值。
（2）当插值点落在已知数据集的外部时的插值称为外插， 要估计外插函数值很难。
例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图 纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘 轮廓线的部分数据如下:
x=[0 4.74 9.05 19 38 57 76 95 114 133 152 171 190] y=[0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 9.69 7.03 3.99 0] xi=[0:0.001:190] yi=interp1(x,y,xi,'spline') plot(xi,yi)
t=0:2:24 T=[12 9 9 10 18 24 28 27 25 20 18 15 13] plot(t,T,'*') ti=0:1/3600:24 T1i=interp1(t,T,ti) plot(t,T,'*',ti,T1i,'r-') T2i=interp1(t,T,ti,'spline') plot(t,T,'*',ti,T1i,'r-',ti,T2i,'g-')
例7 天文学家在1914年8月份的7次观测中, 测得地球与金 星之间距离(单位: m), 并取其常用对数值与日期的一组历 史数据如下所示, 试推断何时金星与地球的距离(单位: m)
的对数值为 9.9352.
日期 18 20 22 24 26 28 30
距离 9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150 对数 解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数 值之间, 所以对上述数据用三次样条插值加细为步长为1 的数据:
y2
1.0000 0.8205 0.5000 0.2973 0.2000 0.1401 0.1000 0.0745 0.0588 0.0484 0.0385
例 5 在一天24h内, 从零点开始每间隔2h测得的环境温度为
12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 (单位: C ) 推测在每1s时的温度. 并描绘温度曲线. 解 在命令窗口输入:
例3. 已知实验数据如表。
x y
0 0.25 0.50 0.75 1.00 0.9162 0.8109 0.6931 0.5596 0.4055
计算插值点x0=0.6处的函数值y0。
>> x=[0 0.25 0.50 0.75 1.00]; >> y=[0.9162 0.8109 0.6931 0.5596 0.4055]; >> x0=0.6; >> y01=interp1(x,y,x0); >> y02=interp1(x,y,x0,'nearest'); >> y03=interp1(x,y,x0,'pchip'); >> y04=interp1(x,y,x0,'spline'); >> y01,y02,y03,y04
y
1.0000 0.8000 0.5000 0.3077 0.2000 0.1379 0.1000 0.0755 0.0588 0.0471 0.0385
y1
1.0000 0.7500 0.5000 0.3500 0.2000 0.1500 0.1000 0.0794 0.0588 0.0486 0.0385
Lagrange插值
方法介绍 对给定的n个插值点 x1, x2 , , xn 及对应的函数 值 y1 , y2 , , y，利用 n次Lagrange插值多项式， n 则对插值区间内任意x的函数值y可通过下式求的：
y ( x) yk (
k 1 j 1 j k
n
n
x xj xk x j
四种插值方法比较
函数
格式
功能
zi=interp2(x,y,z,xi,yi)
二维插值。Z为由已知点的值组成的 矩阵，参量x与y是与z同维的已知点 的矩阵，且必须是单调的。xi与yi为 需要插值的点。若xi与yi中有在x与y 范围之外的点，则相应地返回NaN。
此格式默认x=1:n、y=1:m，其中 [m,n]=size(z)。再按 zi=interp2(x,y,z,xi,yi)情形进行计算。
1 在[-5, 5]上, 用n=11个等距分点作分段线 2 1 x
例4 对 y
性插值和三次样条插值, 用m=21个插值点作图,比较结果. x
0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 4.5000 5.0000
1 在[-5, 5]上, 用n=11个等距分点作分段线 2 1 x
例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图
纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘 轮廓线的部分数据如下: x 0 4.74 9.05 19 38 57 76 95 114 133 y 0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 6.69 x 152 171 190 y 7.03 3.99 0
Matlab常用数据插值函数及功能 函数 格式 yi=interp1(x,y,xi) yi=interp1(y,xi) Interp1 功能 一维插值。对一组点(x,y)进行插 值，计算插值点xi的函数值。 此格式默认x=1:n，其中n是向量 y的长度，y为矩阵时， n=size(y,1)
yi=interp1(x,y,xi,method) Method为指定的插值算法。 Nearest最近邻点插值，linear线 性插值(默认方式)，spline三次 样条函数插值，cubic三次插值。
曲线拟合
据人口统计年鉴，知我国从1949 年至1994年人口数据资料如下： (人口数单位为：百万)
y ( x) yk (
k 1 j 1 j k
n
n
x xj xk x j
)
例：给出f(x)=ln(x)的数值表，用Lagrange计算ln(0.54)的 近似值。
>> x=[0.4:0.1:0.8]; >> y=log(x); >> x0=0.54; >> y0=lagrange(x,y,x0); >> ys=log(x0); >> y0,ys y0 = -0.616142610505419 ys = -0.616186139423817
>> x=linspace(0,2*pi,11); >> y=sin(x).*exp(-x/5); >> xi=linspace(0,2*pi,21); >> yi=interpft(y,21); >> plot(x,y,'o',xi,yi); >> legend('Original','Curve by interpft')
19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 26.0000 27.0000 28.0000 29.0000 30.0000 9.9581 9.9544 9.9506 9.9468 9.9430 9.9391 9.9312 9.9272 9.9232 9.9191 9.9150
0.69
145 0.91
0.74
150 0.92
0.78
155 0.93
0.82
160 0.95
0.84
165 0.96
0.86
170 0.98
0.88
175 0.99
0.89
180 1
包角系 数
>>a1=[90,100,110,120,125,130,135,140,145,150,155,160,165,170,175,180]; >>a2=[0.69,0.74,0.78,0.82,0.84,0.86,0.88,0.89,0.91,0.92,0.93,0.95,0.96,0.98,0.99 ,1]; >>ka=interp1(a1,a2,133.5) >>ka=0.8740
MATLAB对已知数据集外部点上函数值的预测都返回NaN， 但可通过为interp1函数添加‘extrap’参数指明也用于外插。 MATLAB的外插结果偏差较大。
例1 对
f x
1 在[-1, 1 9 x2
1]上, 用n=20的等距分点
进行线性插值, 绘制 f(x)及插值函数的图形.
解 在命令窗口输入:
)
MATLAB实现
function yi=lagrange(x,y,xi) yi=zeros(size(xi)); np=length(y); for i=1:np z=ones(size(xi)); for j=1:np if i~=j z=z.*(xi-x(j))/(x(i)-x(j)); end end yi=yi+z*y(i); end
插值就是已知一组离散的数据点集，在集合内部某两 个点之间预测函数值的方法。
插值法是实用的数值方法，是函数逼近的重要方法。在 生产和科学实验中，自变量x与因变量y的函数y = f(x)的 关系式有时不能直接写出表达式，而只能得到函数在若干 个点的函数值或导数值。当要求知道观测点之外的函数值 时，需要估计函数值在该点的值。
x=-1:0.1:1; y=1./(1+9*x.^2); xi=-1:0.1:1; yi=interp1(x,y,xi); plot(x,y,'r-',xi,yi,'*')
例2. 在普通V带设计中，带轮的包角α与包角系数ka之间的关系如 表所示。求α=133.5°时的包角系数ka。
包角与包角系数
包角 ( °) 包角系 数 包角 ( °) 90 100 110 120 125 130 135 140
解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之 间, 所以对上述数据用三次样条插值加细为步长为1的数据: x=[18:2:30] y=[9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150] xi=[18:1:30] yi=interp1(x,y,xi,'spline') A=[xi;yi] A= 18.0000 25.0000 9.9618 9.9352
y01 = 0.639700000000000 y02 = 0.693100000000000 y03 = 0.641818996851421 y04 = 0.641831200000000
例4 对 y
性插值和三次样条插值, 用m=21个插值点作图,比较结果. 解 在命令窗口输入:
n=11, m=21 x=-5:10/(m-1):5 y=1./(1+x.^2) z=0*x x0=-5:10/(n-1):5 y0=1./(1+x0.^2) y1=interp1(x0,y0,x) y2=interp1(x0,y0,x,'spline') [x' y' y1' y2'] plot(x,z,'r',x,y,'k:',x,y1,'b',x,y2,'g') gtext('Piece.linear.'),gtext('Spline'),gtext('y=1/(1+x^2)')
Inteቤተ መጻሕፍቲ ባይዱp2
zi=interp2(z,xi,yi)
zi=interp2(x,y,z,xi,yi,method) 用指定的算法method计算二维插值。 Linear双线性插值(默认方式)， Nearest最邻近插值，spline三次样条 插值，cubic双三次插值。
函数
格式
功能 一维Fourier插值。适用于周期函数生 成数据的插值。X为抽样序列，n为 需要计算的等距点数，y为n点的插 值计算结果。 沿着指定的方向dim进行计算。 三次样条数据插值。用三次样条插 值计算出由向量x与y确定的一元函 数在点xx处的值yy。 返回由向量x与y确定的分段样条多 项式的系数矩阵pp。
y=interpft(x,n) Interpft y=interpft(x,n,dim)
yy=spline(x,y,xx)
spline
pp=spline(x,y)
注意： （1）只对已知数据点集内部的点进行的插值运算称为内插， 可比较准确的估测插值点上的函数值。
（2）当插值点落在已知数据集的外部时的插值称为外插， 要估计外插函数值很难。
例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图 纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘 轮廓线的部分数据如下:
x=[0 4.74 9.05 19 38 57 76 95 114 133 152 171 190] y=[0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 9.69 7.03 3.99 0] xi=[0:0.001:190] yi=interp1(x,y,xi,'spline') plot(xi,yi)
t=0:2:24 T=[12 9 9 10 18 24 28 27 25 20 18 15 13] plot(t,T,'*') ti=0:1/3600:24 T1i=interp1(t,T,ti) plot(t,T,'*',ti,T1i,'r-') T2i=interp1(t,T,ti,'spline') plot(t,T,'*',ti,T1i,'r-',ti,T2i,'g-')
例7 天文学家在1914年8月份的7次观测中, 测得地球与金 星之间距离(单位: m), 并取其常用对数值与日期的一组历 史数据如下所示, 试推断何时金星与地球的距离(单位: m)
的对数值为 9.9352.
日期 18 20 22 24 26 28 30
距离 9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150 对数 解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数 值之间, 所以对上述数据用三次样条插值加细为步长为1 的数据:
y2
1.0000 0.8205 0.5000 0.2973 0.2000 0.1401 0.1000 0.0745 0.0588 0.0484 0.0385
例 5 在一天24h内, 从零点开始每间隔2h测得的环境温度为
12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 (单位: C ) 推测在每1s时的温度. 并描绘温度曲线. 解 在命令窗口输入:
例3. 已知实验数据如表。
x y
0 0.25 0.50 0.75 1.00 0.9162 0.8109 0.6931 0.5596 0.4055
计算插值点x0=0.6处的函数值y0。
>> x=[0 0.25 0.50 0.75 1.00]; >> y=[0.9162 0.8109 0.6931 0.5596 0.4055]; >> x0=0.6; >> y01=interp1(x,y,x0); >> y02=interp1(x,y,x0,'nearest'); >> y03=interp1(x,y,x0,'pchip'); >> y04=interp1(x,y,x0,'spline'); >> y01,y02,y03,y04
y
1.0000 0.8000 0.5000 0.3077 0.2000 0.1379 0.1000 0.0755 0.0588 0.0471 0.0385
y1
1.0000 0.7500 0.5000 0.3500 0.2000 0.1500 0.1000 0.0794 0.0588 0.0486 0.0385
Lagrange插值
方法介绍 对给定的n个插值点 x1, x2 , , xn 及对应的函数 值 y1 , y2 , , y，利用 n次Lagrange插值多项式， n 则对插值区间内任意x的函数值y可通过下式求的：
y ( x) yk (
k 1 j 1 j k
n
n
x xj xk x j
四种插值方法比较
函数
格式
功能
zi=interp2(x,y,z,xi,yi)
二维插值。Z为由已知点的值组成的 矩阵，参量x与y是与z同维的已知点 的矩阵，且必须是单调的。xi与yi为 需要插值的点。若xi与yi中有在x与y 范围之外的点，则相应地返回NaN。
此格式默认x=1:n、y=1:m，其中 [m,n]=size(z)。再按 zi=interp2(x,y,z,xi,yi)情形进行计算。
1 在[-5, 5]上, 用n=11个等距分点作分段线 2 1 x
例4 对 y
性插值和三次样条插值, 用m=21个插值点作图,比较结果. x
0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000 3.5000 4.0000 4.5000 5.0000
1 在[-5, 5]上, 用n=11个等距分点作分段线 2 1 x
例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图
纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘 轮廓线的部分数据如下: x 0 4.74 9.05 19 38 57 76 95 114 133 y 0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 6.69 x 152 171 190 y 7.03 3.99 0
Matlab常用数据插值函数及功能 函数 格式 yi=interp1(x,y,xi) yi=interp1(y,xi) Interp1 功能 一维插值。对一组点(x,y)进行插 值，计算插值点xi的函数值。 此格式默认x=1:n，其中n是向量 y的长度，y为矩阵时， n=size(y,1)
yi=interp1(x,y,xi,method) Method为指定的插值算法。 Nearest最近邻点插值，linear线 性插值(默认方式)，spline三次 样条函数插值，cubic三次插值。
曲线拟合
据人口统计年鉴，知我国从1949 年至1994年人口数据资料如下： (人口数单位为：百万)
y ( x) yk (
k 1 j 1 j k
n
n
x xj xk x j
)
例：给出f(x)=ln(x)的数值表，用Lagrange计算ln(0.54)的 近似值。
>> x=[0.4:0.1:0.8]; >> y=log(x); >> x0=0.54; >> y0=lagrange(x,y,x0); >> ys=log(x0); >> y0,ys y0 = -0.616142610505419 ys = -0.616186139423817
>> x=linspace(0,2*pi,11); >> y=sin(x).*exp(-x/5); >> xi=linspace(0,2*pi,21); >> yi=interpft(y,21); >> plot(x,y,'o',xi,yi); >> legend('Original','Curve by interpft')
19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 26.0000 27.0000 28.0000 29.0000 30.0000 9.9581 9.9544 9.9506 9.9468 9.9430 9.9391 9.9312 9.9272 9.9232 9.9191 9.9150
0.69
145 0.91
0.74
150 0.92
0.78
155 0.93
0.82
160 0.95
0.84
165 0.96
0.86
170 0.98
0.88
175 0.99
0.89
180 1
包角系 数
>>a1=[90,100,110,120,125,130,135,140,145,150,155,160,165,170,175,180]; >>a2=[0.69,0.74,0.78,0.82,0.84,0.86,0.88,0.89,0.91,0.92,0.93,0.95,0.96,0.98,0.99 ,1]; >>ka=interp1(a1,a2,133.5) >>ka=0.8740
MATLAB对已知数据集外部点上函数值的预测都返回NaN， 但可通过为interp1函数添加‘extrap’参数指明也用于外插。 MATLAB的外插结果偏差较大。
例1 对
f x
1 在[-1, 1 9 x2
1]上, 用n=20的等距分点
进行线性插值, 绘制 f(x)及插值函数的图形.
解 在命令窗口输入:
)
MATLAB实现
function yi=lagrange(x,y,xi) yi=zeros(size(xi)); np=length(y); for i=1:np z=ones(size(xi)); for j=1:np if i~=j z=z.*(xi-x(j))/(x(i)-x(j)); end end yi=yi+z*y(i); end