数值计算3插值和曲线拟合.docx
《数值计算3插值和曲线拟合.docx》由会员分享，可在线阅读，更多相关《数值计算3插值和曲线拟合.docx（15页珍藏版）》请在冰豆网上搜索。
数值计算3插值和曲线拟合
数值计算3－插值和曲线拟合
插值法是实用的数值方法，是函数逼近的重要方法。
在生产和科学实验中，自变量x与因变量y的函数y=f（x）的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。
当要求知道观测点之外的函数值时，需要估计函数值在该点的值。
如何根据观测点的值，构造一个比较简单的函数y=φ（x），使函数在观测点的值等于已知的数值或导数值。
用简单函数y=φ（x）在点x处的值来估计未知函数y=f（x）在x点的值。
寻找这样的函数φ（x），办法是很多的。
φ（x）可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ（x）可以是任意光滑（任意阶导数连续）的函数或是分段函数。
函数类的不同，自然地有不同的逼近效果。
在许多应用中，通常要用一个解析函数（一、二元函数）来描述观测数据。
根据测量数据的类型：
1．测量值是准确的，没有误差。
2．测量值与真实值有误差。
这时对应地有两种处理观测数据方法：
1．插值或曲线拟合。
2．回归分析（假定数据测量是精确时，一般用插值法，否则用曲线拟合）。
MATLAB中提供了众多的数据处理命令。
有插值命令，有拟合命令，有查表命令。
一维插值
插值定义为对数据点之间函数的估值方法，这些数据点是由某些集合给定。
当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。
例如，当数据点是某些实验测量的结果或是过长的计算过程时，就有这种情况。
interp1（x,y,xi,method）
x和y为既有数据的向量，其长度必须相同。
xi为要插值的数据点向量。
method插值方法,‘nearest’/‘linear’/‘cubic’/‘spline’之一，
分别为最近点插值/线性插值/分段三次Hermite插值/三次样条插值。
例
x=[1.02.03.04.05.0];%输入变量数据x
y=[11.216.520.426.330.5];%输入变量数据y
x1=2.55;%输入待插值点x
y11=interp1（x,y,x1,'nearest'）%最近点插值方法的插值结果
y12=interp1（x,y,x1,'linear'）%线性插值方法的插值结果
y13=interp1（x,y,x1,'cubic'）%三次Hermite插值方法的插值结果
y14=interp1（x,y,x1,'spline'）%样条插值方法的插值结果
y11=
20.4000
y12=
18.6450
y13=
18.6028
y14=
18.4874
plot（x,y）
或许最简单插值的例子是MATLAB的作图。
按缺省，MATLAB用直线连接所用的数据点以作图。
这个线性插值猜测中间值落在数据点之间的直线上。
当然，当数据点个数的增加和它们之间距离的减小时，线性插值就更精确。
例如：
x1=linspace（0,2*pi,60）;
x2=linspace（0,2*pi,6）;
plot（x1,sin（x1）,x2,sin（x2）,'-'）
xlabel（'x'）,ylabel（'sin（x）'）,title（'LinearInterpolation'）
上图是sin函数的两个图，一个在数据点之间用60个点，它比另一个只用6个点更光滑和更精确。
根据所作的假设，有多种插值。
而且，可以在一维以上空间中进行插值。
即如果有反映两个变量函数的插值，z=f（x,y），那么就可在x之间和在y之间，找出z的中间值进行插值。
MATLAB在一维函数interp1和在二维函数interp2中，提供了许多的插值选择。
其中的每个函数将在下面阐述。
为了说明一维插值，考虑下列问题，12小时内，一小时测量一次室外温度。
数据存储在两个MATLAB变量中。
hours=1:
12;%indexforhourdatawasrecorded
temps=[589152529313022252724];
%recordedtemperatures
plot（hours,temps,hours,temps,'+'）%viewtemperatures
title（'Temperature'）
xlabel（'Hour'）,ylabel（'DegreesCelsius'）
在线性插值下室外温度曲线
正如上图看到的，MATLAB画出了数据点线性插值的直线。
为了计算在任意给定时间的温度，人们可试着对可视的图作解释。
另外一种方法，可用函数interp1。
t=interp1（hours,temps,9.3）%estimatetemperatureathour=9.3
t=
22.9000
t=interp1（hours,temps,4.7）%estimatetemperatureathour=4.7
t=
22
t=interp1（hours,temps,[3.26.57.111.7]）%findtempatmanypoints!
t=
10.200030.000030.900024.9000
interp1的缺省用法是由interp1（x,y,xo）来描述，这里x是独立变量（横坐标），y是应变量（纵坐标），xo是进行插值的一个数值数组。
另外，该缺省的使用假定为线性插值。
若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。
最常用的方法是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致。
这种类型的插值被称为3次样条或简称为样条。
函数interp1也能执行3次样条插值。
hours=1:
12;
temps=[589152529313022252724];
plot（hours,temps,hours,temps,'+'）
title（'Temperature'）
xlabel（'时间轴'）,ylabel（'温度变化轴'）
x1=[9.32.5];
t1=interp1（hours,temps,x1,'spline'）
t1=
21.85778.2364
x2=[3.26.57.111.7];
t3=interp1（hours,temps,x2,'spline'）
t3=
9.673430.042731.175525.3820
t3=
9.673430.042731.175525.3820
注意，样条插值得到的结果，与上面所示的线性插值的结果不同。
因为插值是一个估计或猜测的过程，其意义在于，应用不同的估计规则导致不同的结果。
一个最常用的样条插值是对数据平滑。
也就是，给定一组数据，使用样条插值在更细的间隔求值。
例如，
hours=1:
12;
temps=[589152529313022252724];
h=1:
0.1:
12;%estimatetemperatureevery1/10hour
t=interp1（hours,temps,h,'spline'）;
plot（hours,temps,'-',hours,temps,'+',h,t）
%plotcomparativeresults
title（'SpringfieldTemperature'）
xlabel（'Hour'）,ylabel（'DegreesCelsius'）
在不同插值下室外温度曲线
在上图中，虚线是线性插值，实线是平滑的样条插值，标有'+'的是原始数据。
如要求在时间轴上有更细的分辨率，并使用样条插值，我们有一个更平滑、但不一定更精确地对温度的估计。
尤其应注意，在数据点，样条解的斜率不突然改变。
作为这个平滑插值的回报，3次样条插值要求更大量的计算，因为必须找到3次多项式以描述给定数据之间的特征。
在讨论二维插值之前，了解interp1所强制的二个强约束是很重要的。
首先，人们不能要求有独立变量范围以外的结果，例如，interp1（hours,temps,13.5）导致一个错误，因为hours在1到12之间变化。
其次，独立变量必须是单调的。
即独立变量在值上必须总是增加的或总是减小的。
在我们的例子里，hours是单调的。
然而，如果我们已经定义独立变量为一天的实际时间，
»time_of_day=[7:
121:
6]%startat7AM,endat6PM
time_of_day=
789101112123456
则独立变量将不是单调的，因为time_of_day增加到12，然后跌到1，再然后增加。
如果用time_of_day代替interp1中的hours，将会返回一个错误。
同样的理由，人们不能对temps插值来找出产生某温度的时间（小时），因为temps不是单调的。
二维插值
interp2（x,y,z,xi,yi,method）
x,y,z为既有数据的向量，其长度必须相同。
xi,yi为要插值的数据点。
method为插值方法，可为'nearest'、'linear'、'cubic'之一，
分别为最近点插值、双线性插值、双三次插值
%例：
x=[0.00.51.01.52.02.53.0];%输入变量数据x
y=[0.00.51.01.52.02.53.0];%输入变量数据y
%输入变量数据z
z=[100991009910099100
10099999910099100
9999989810099100
10098979799100103
1011009898100102104
102103101100102106106
99102100100103108103];
x1=2.3;y1=2.8;%输入待插值点[x1,y1]
z11=interp2（x,y,z,x1,y1,'nearest'）
%最近点插值方法的插值结果
z12=interp2（x,y,z,x1,y1,'bilinear'）
%双线性插值插值方法的插值结果
z13=interp2（x,y,z,x1,y1,'bicubic'）
%双三次插值插值方法的插值结果
z11=
108
z12=
105.3600
z13=
105.9744
二维插值是基于与一维插值同样的基本思想。
然而，正如名字所隐含的，二维插值是对两变量的函数z=f（x,y）进行插值。
为了说明这个附加的维数，考虑一个问题。
设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。
采集了下列的数据：
width=1:
5;%indexforwidthofplate（i.e.,thex-dimension）
depth=1:
3;%indexfordepthofplate（i,e,,they-dimension）
temps=[8281808284;7963616581;8484828586]
%temperaturedata
temps=
8281808284
7963616581
8484828586
如同在标引点上测量一样，矩阵temps表示整个平板的温度分布。
temps的列与下标depth或y-维相联系，行与下标width或x-维相联系（见图）。
为了估计在中间点的温度，我们必须对它们进行辨识。
wi=1:
0.2:
5;%estimateacrosswidthofplate
d=2;%atadepthof2
zlinear=interp2（width,depth,temps,wi,d）;%linearinterpolation
zcubic=interp2（width,depth,temps,wi,d,'cubic'）;
%cubicinterpolation
plot（wi,zlinear,'-',wi,zcubic,'r'）%plotresults
xlabel（'WidthofPlate'）,ylabel（'DegreesCelsius'）
title（['TemperatureatDepth='num2str（d）]）
在深度d=2处的平板温度
另一种方法，我们可以在两个方向插值。
先在三维坐标画出原始数据，看一下该数据的粗糙程度（见图11.7）。
mesh（width,depth,temps）%usemeshplot
xlabel（'WidthofPlate'）,ylabel（'DepthofPlate'）
zlabel（'DegreesCelsius'）,axis（'ij'）,grid
平板温度
然后在两个方向上插值，以平滑数据。
width=1:
5;
depth=1:
3;
temps=[8281808284;7963616581;8484828586];
di=1:
0.2:
3;%choosehigherresolutionfordepth
wi=1:
0.2:
5;%choosehigherresolutionforwidth
zcubic=interp2（width,depth,temps,wi,di,'cubic'）;%cubic
mesh（wi,di,zcubic）
xlabel（'WidthofPlate'）,ylabel（'DepthofPlate'）
zlabel（'DegreesCelsius'）,axis（'ij'）,grid
上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪。
interp2的基本形式是interp2（x,y,z,xi,yi,method）。
这里x和y是两个独立变量，z是一个应变量矩阵。
x和y对z的关系是
z（i,:
）=f（x,y（i））和z（:
j）=f（x（j）,y）.
也就是，当x变化时，z的第i行与y的第i个元素y（i）相关，当y变化时，z的第j列与x的第j个元素x（j）相关，。
xi是沿x-轴插值的一个数值数组；yi是沿y-轴插值的一个数值数组。
图11.8二维插值后的平板温度
可选的参数method可以是'linear'，'cubic'或'nearest'。
在这种情况下，cubic不意味着3次样条，而是使用3次多项式的另一种算法。
linear方法是线性插值，仅用作连接图上数据点。
nearest方法只选择最接近各估计点的粗略数据点。
在所有的情况下，假定独立变量x和y是线性间隔和单调的。
关于这些方法的更多的信息，可请求在线帮助。
二．曲线拟合
使用多项式曲线拟合
p=polyfit（x,y,n）
p为最小二乘意义上拟合多项式的相关系数。
x,y为数据点向量。
n为多项式阶次。
y=polyval（p,x）
y为自变量为x、系数为p的多项式值.
x为矩阵或向量时对其每一点进行操作.
p为一向量，各元素为多项式系数（降幂排列）.
曲线拟合涉及回答两个基本问题：
最佳拟合意味着什么？
应该用什么样的曲线？
可用许多不同的方法定义最佳拟合，并存在无穷数目的曲线。
所以，从这里开始，我们走向何方？
正如它证实的那样，当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。
数学上，称为多项式的最小二乘曲线拟合。
如果这种描述使你混淆，再研究图。
虚线和标志的数据点之间的垂直距离是在该点的误差。
对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。
这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。
最小二乘这个术语仅仅是使误差平方和最小的省略说法。
在MATLAB中，函数polyfit求解最小二乘曲线拟合问题。
为了阐述这个函数的用法，让我们以上图中的数据开始。
x=[0.1.2.3.4.5.6.7.8.91];
y=[-.4471.9783.286.167.087.347.669.569.489.3011.2];
为了用polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。
如果我们选择n=1作为阶次，得到最简单的线性近似。
通常称为线性回归。
相反，如果我们选择n=2作为阶次，得到一个2阶多项式。
现在，我们选择一个2阶多项式。
n=2;%polynomialorder
p=polyfit（x,y,n）
p=
-9.810820.1293-0.0317
polyfit的输出是一个多项式系数的行向量。
其解是y=－9.8108x2＋20.1293x－0.0317。
为了将曲线拟合解与数据点比较，让我们把二者都绘成图。
xi=linspace（0,1,100）;%x-axisdataforplotting
z=polyval（p,xi）;
为了计算在xi数据点的多项式值，调用MATLAB的函数polyval。
plot（x,y,'o',x,y,xi,z,':
'）
画出了原始数据x和y，用'o'标出该数据点，在数据点之间，再用直线重画原始数据，并用点':
'线，画出多项式数据xi和z。
»xlabel（'x'）,ylabel（'y=f（x）'）,title（'SecondOrderCurveFitting'）
将图作标志。
这些步骤的结果表示于前面的图中。
例１．设ｘ＝[0.0 0.3 0.8 1.1 1.6 2.3]'，y=[0.50 0.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.91.2]时对应的y1。
x=[0.00.30.81.11.62.3]';%输入变量数据x
y=[0.500.821.141.251.351.40]';%输入变量数据y
p=polyfit（x,y,2）%对x,y用二次多项式拟合,得到系数p
x1=[0.91.2];%输入点x1
y1=polyval（p,x1）%估计x1处对应的y1
p=
-0.23870.91910.5318
y1=
1.16561.2909
使用指定函数进行曲线拟合
例2．使用试根据土的压缩实验数据（ｅ－ｐ数据），用双曲线模型
确定模型参数ａ，ｂ（压力ｐ＝0时的孔隙比，可用天然孔隙比e0代替），已知p=0,0.05,0.1,0.2,0.4,0.6,0.8,1.2MPa对应的孔隙比分别为e=1.335,1.253,1.180,1.058,0.887,0.803,0.752,0.685。
Data=[0.00001.335;0.05001.253;0.10001.180;0.20001.058;0.40000.887;
0.60000.803;0.80000.752;1.20000.685];
p=Data（:
1）;
e=Data（:
2）;
plot（p,e,'g-'）;holdon;
x0=[10]';
x=lsqnonlin（'f0405',x0）
e1=e
（1）;
e=e1-p./（x
（1）+x
（2）*p）;
plot（p,e,'bo-'）;
Optimizationterminated:
relativefunctionvalue
changingbylessthanOPTIONS.TolFun.
x=
0.4833
1.1104
functionf=f0405（x,Data）
Data=[0.00001.335;0.05001.253;0.10001.180;0.20001.058;0.40000.887;
0.60000.803;0.80000.752;1.20000.685];
p=Data（:
1）;
e=Data（:
2）;
e1=e
（1）;
z=e1-p./（x
（1）+x
（2）*p）;
f=z-e;
作业：
1．已知ｘ＝〔0.0 0.3 0.8 1.1 1.6 2.3〕，ｙ＝〔0.40 0.56 1.14 1.31 1.56 2.10〕， 试用MatLab对数据点对ｘ－ｙ进行二次多项式拟合，并使用不同插值方法计算ｘ＝２处对应的ｙ。
2．实测数据如下：
ｘ＝〔0.0 0.5 1.0 1.5〕；ｙ＝〔0.0 0.5 1.0 1.5〕；
ｚ＝〔50 49 50 49；49 49 49 49；49 48 48 47；48 48 47 46〕
试用MatLab求ｘ1＝1.2、ｙ1＝1.4处的插值结果。
３．求解方程组
升级会员 