插值与数据拟合建模

插值与数据拟合就是通过一些已知数据去确定某类函数的参数或寻找某个近似函数，使所得的函数与已知数据具有较高的精度，并且能够使用数学分析的工具分析数据所反映的对象的性质．几种常用的方法：

1、一般插值法

2、样条插值法

3、最小二乘曲线

4、曲面的拟合数据拟合与插值建模

上大学二年级的小华正在做作业，“爸爸，计算这道题要用到sin，可是我的计算器坏了，怎么办。”当工程师的老张从厚厚的一摞旧书底下抽出一本数学用表来，“给你，这是我念大学时用的，那时候啊，计算器听都没听说过。”小华拿着表翻了一会儿，无奈的说：“表上每才有一个函数值，这里只sin和sin““表中没有的，都可以用插值方法计算”“插值！我们的数学实验课就要学了，不过，今天我要先自己想个办法，用这个算出sin”这本四位数学用表给出sin＝0.576，sin=0.5783。小华认为在sin到sin这样小的范围内，正弦可以近似为线性函数，于是很容易地得到Sin=0.576+(0.5783-0.5760)×0.6=0.5774

聪明的小华用的这个办法是一种插值方法——分段线性插值。实际上，插值可以理解为，要根据一个用表格表示的函数，计算表中没有的函数值。表中有的，如（sin，0.5760）(sin,0.5783)称为节点；要计算的，如sin，称为插值点，结果（0.5774）即为插值。小华作的线性函数为插值函数，插值函数所表示的直线当然要通过节点。

插值最初来源于天体计算——由若干观测值（即节点）计算任意时刻星球的位置（即插值点和插值）——的需要。现在，虽然人们已很少需要用它从函数表计算函数值了，但是插值仍然在诸如机械加工等工程技术和数据处理等科学研究中有着许多直接的应用，另一方面，插值又是数值微分、数值积分、常微分方程数值等数值计算的基础。

几天后，小华在物理实验里又碰到一个看起来非常类似的问题：有一只对温度敏感的电阻，已经测得了一组温度T和电阻R数据。

现在想知道600C时的电阻多大。温度t(0C)20.532.751.073.095.7电阻R()7658268739421032

小华征求老师的意见，老师给了他两点提示：一是在直角坐标系中把5个点（T，R）画一下，看看电阻R和温度T之间大致有什么关系；二是测量数据总有相当大的误差，这与用函数表作插值计算应该有不同之处吧（虽然函数表也存在舍入误差，但很小，可以认为表中数值是精确的）通过图形小华看到，R与T大致呈直线关系，于是用手画了一条靠近5个点的直线，又想起中学物理学过，金属材料的电阻率与温度成正比，从而确定R与T的关系应该是

R=at+b

其中a，b为待定常数。

正是由于测量误差的存在，由R=at+b表示的直线不可能通过全部5个点，所以，与插值曲线要通过全部节点不同，小华打算作一条尽量靠近所有的点的直线，求出a，b待定常数，由此计算t=600C

的R就十分简单了。

根据一组（二组）数据，即平面上的若干点，确定一个一元函数，即曲线，使这些节点与曲线总体来说尽量接近，这就是曲线拟合。函数值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法是完全不同的。

数据拟合1.拟合的基本原理；2.最小二乘拟合；3.用Matlab作最小二乘拟合；4.如何用拟合解决实际问题。

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01

对某人用快速静脉注射方式一次性注射某种药物300mg后，经过时间t采集血样，测得血药浓度c如下表：求血药浓度随时间的变化规律c(t).半对数坐标系(semilogy)下的图形Log10c(t)=at+b引例1：血药浓度的变化规律曲线拟合问题的提法

已知一组（二维）数据，即平面上n个点（xi,yi)i=1,…n,

寻求一个函数（曲线）y=f(x),

使f(x)

在某种准则下与所有数据点最为接近，即曲线拟合得最好。

+++++++++xyy=f(x)(xi,yi)

i

i为点（xi,yi)与曲线y=f(x)的距离最小二乘拟合

第一步：先选定一类函数f(x，a1，a2，

…，am)

其准则为（最小二乘准则）：使n个点（xi,yi)与曲线y=f(x，a1，a2，

…，am)的距离

i的平方和最小

。其中

a1,a2,…am

为待定常数。f可以为一些简单的“基函数”（如幂函数，三角函数等等）的线性组合：第二步：确定参数a1,a2,…am，记问题归结为，求

a1,a2,…am

使

J(a1,a2,…am)

最小。这样的拟合称为最小二乘拟合。

除了最小二乘准则（即各点误差的平方和最小），你认为还可以用怎样的拟合准则？比较起来，最小二乘准则有什么优点？思考最小二乘拟合函数f(x，a1，…am)的选取++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1exp(a2x)+++++f=a1exp(a2x)1.通过机理分析建立数学模型来确定f；2.将数据(xi,yi)i=1,…n作图，通过直观判断确定f：2.作一般的最小二乘曲线拟合，可利用已有程序curvefit,其调用格式为：

a=curvefit(‘f’,a0,x,y)

1.作多项式f(x)=a1xm+…+amx+am+1函数拟合,可利用已有程序polyfit,其调用格式为:a=polyfit(x,y,m)用MATLAB作最小二乘拟合数据点拟合多项式次数系数注：f为拟合函数y=f(a,x)的函数M—文件，f(a,x)为拟合函数。数据点待定常数a的初值函数M文件用MATLAB作多项式最小二乘拟合2.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.49181.选取函数R=

a1t+a2温度t(0C)20.532.751.073.095.7电阻R()7658268739421032例.由数据拟合R=f(t)用MATLAB作最小二乘曲线拟合例：用函数f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)拟合下列数据点：xdata=[0:.1:2]ydata=[5.89553.56392.51731.97901.89901.39381.13591.00961.03430.84350.68560.61000.53920.39460.39030.54740.34590.13700.22110.17040.2636]用命令curvefit(‘f’,a0,x,y)

拟合的应用——参数辨识

数学建模的方法：机理分析和测试分析。机理分析是根据对客观事物特性的认识，找出反映内部机理的数量规律，建立的模型常有明确的物理意义。测试分析将研究的对象看作一个“黑箱”，通过对实验数据的统计分析，找出与数据拟合得最好的模型。机理分析——>模型结构实验数据——>未知参数范例：薄膜渗透率的测定

一、问题：某种医用薄膜，具有从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜被物质分子穿透的能力。测定方法：用面积为S的薄膜将容器分成体积分别为的两部份，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数K表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定K。VAVBS二、问题分析考察时段[t，t+Δt]薄膜两侧容器中该物质质量的变化。设，对容器的B部分溶液浓度的测试结果如下表：（浓度单位）

1）在容器的一侧，物质质量的增加是由于另一侧的物质向该侧渗透的结果，因此物质质量的增量应等于另一侧的该物质向这侧的渗透量。

以容器A侧为例，在时段[t，t+Δt]物质质量的增量为：分别表示在时刻t膜两侧溶液设的浓度，浓度单位:

由于平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数为K。

因此，在时段[t，t+Δt]，从B侧渗透至A侧的该物质的质量为：于是有：两边除以Δt，并令Δt→0取极限再稍加整理即得：分别表示在初始时刻两侧溶液的浓度其中（1）2)注意到整个容器的溶液中含有该物质的质量不变,与初始时刻该物质的含量相同，因此从而：加上初值条件：代入式（1）得：便可得出CB(t)的变化规律，从而根据实验数据进行拟合，估计出参数K,。三、数学模型假设：1）薄膜两侧的溶液始终是均匀的；2）平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。3）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。基于假设和前面的分析，B侧的浓度CB(t)应满足如下微分方程和初始条件：四、求解方法：1.函数拟合法前面得到的模型是一个带初值的一阶线性微分方程，解之得：问题归结为利用CB在时刻tj的测量数据Cj(j=1,2,...,N)来辨识K和。引入从而

用函数CB(t)来拟合所给的实验数据，从而估计出其中的参数a,b,K。将代入上式有：用MATLAB软件进行计算.1）编写函数M-文件nongdu.mfunctionf=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);其中x(1)=a;x(2)=b;x(3)=k;2)在工作空间中执行以下命令(test1.m)tdata=linspace(100,1000,10);cdata=[4.544.995.355.655.906.10...6.266.396.506.59];x0=[0.2,0.05,0.05];x=curvefit(‘nongdu’,x0,tdata,cdata)3)输出结果:x=0.007-0.0030.1012

即k=0.1012,a=0.007,b=-0.003,进一步求得：2.非线性规划法

利用CB在时刻tj的测量数据Cj(j=1,2,...,N)来辨识K和。问题可转化为求函数即求函数的最小值点（K，a，b）。3.导函数拟合法前面得到的微分方程为：令上式变为：这可以看作随CB的变化规律(j=1,2,...,N)若知道一组数据则可用最小二乘拟合的方法来求出函数中的未知参数K和h。即为求参数K,a使下列误差函数达到最小：

该问题等价于用函数f(K,a,CB)=K(0.01a-0.02CB)来拟合数据(j=1,2,...,N)用MATLAB软件进行计算.％求数据点(j=1,2,...,N)tdata=linspace(100,1000,10);cdata=1e-05.*[454499535565590...

610626639650659];[d,ifail]=e01bef(tdata,cdata);[cj,dcj]=e01bgf(tdata,cdata,d,tdata);1）编写函数M-文件baomof.mfunctionf=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中x(1)=K;x(2)=h2)编写命令M文件(baomo21.m)3)输出结果:x=0.10090.014

即k=0.1009,h=0.014％作函数拟合x0=[0.2,0.1];x=curvefit('baomof',x0,cdata,dcj')4.线性化迭代法前面带初始条件的一阶线性微分方程的解为其中：

如果得到了参数K的一个较好的近似值K*，则将CB(t)关于K在K*处展开，略去K的二次及以上的项得CB(t)的一个近似式通过极小化确定a,b,d,再由

K=d/0.02b得到K*的修正值

K。K*K*-K,得到K的一个新的近似值，用同样的方法再求新的修正值

K。这个过程可以不断重复，直到修正值足够小为止。1）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；2）当K的初值取为k=0.2时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。五、结果及误差分析

几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方和。注：导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。因此可将导函数拟合法和线性化迭代法结合起来使用，把前者得到的参数K的值作为迭代法中K的初值，这样可使迭代法收敛或收敛更快。3）取K的初值为k=0.1009,只一次迭代就得到2）中的最后结果。函数拟合法的拟合效果求解参数辨识模型的方法：函数拟合；非线性规划；导函数拟合；线性化迭代；其它方法。用Logistic模拟水稻叶伸长生长时间11.82.63.44.14.85.46.16.87.48.1重量0.30.50.91.42.53.24.37.610.114.418.5时间8.89.410.110.811.712.413.114.415.115.7

重量23.025.230.433.738.841.743.744.845.545.3

生长观测记录数据模型表达式：程序！关于polyfit命令命令：p=polyfit(x,y,n)（1）x与y为模拟数据（2）n为拟合多项式的次数（3）当n=1时为用最小二乘法进行直线拟合（4）得到的向量p为长度n+1向量，对应p的分量依次是次数从高到底各多项式系数用Richard模拟

水稻叶伸长生长关于inline函数例如：y=inline(‘sin(x)-cos(x)’,’x’)输入y(0)，可得：-1作图：x=0:0.1:2*pi;plot(x,y(x))1、插值问题：不知道某一函数f(x)在待定范围[a,b]上的具体表达式，而只能通过实验测量得到该函数在一系列点a≤x1,x2,...,xn≤b上的值

y0,y1,y2,...,yn，需要找一个简单的函数P(x)

来近似地代替f(x)，要求满足：

P(xi)=yi(i=1,2,...,n)2、插值函数：P(x)3、插值法：求插值函数P(x)的方法

插值二、常用插值函数1、多项式函数2、样条函数1、多项式插值方法（1）n次代数插值（2）拉格朗日插值几点说明：（1）拉格朗日插值基函数仅与节点有关，而与被插值函数f(x)无关。（2）拉格朗日插值多项式仅由数对(xi,yi)(i＝

1,2,…,n)确定，而与数对排列次序无关。（3）多项式插值除了上述插值法外还有其它插值法，如newton插值法、hermite插值法等。2、样条插值方法（1）样条函数——m次半截幂函数（2）k次B样条或k次基本样条函数的定义（一）广泛使用的样条函数（1）广泛采用：二次样条、三次样条及B样条（2）力学意义：

A:二次样条在力学上解释为集中力偶作用下的弹性细梁挠度曲线。

B:弹性细梁受集中载荷作用形成的挠度曲线，在小挠度的情况下，恰好表示为三次样条函数，集中载荷的作用点，恰好就是三次样条函数的节点。（1）二次样条的定义

设[a,b]的一个划分：a=x0<x1,x2,...,xn=b，函数f(x)各节点的值分别为：

f(xi)=yi(i=1,2,...,n)

如果二次样条函数：满足：S(xi)=yi(i=1,2,...,n)（2）三次样条函数的定义

设[a,b]的一个划分：a=x0<x1,x2,...,xn=b，函数f(x)各节点的值分别为：

f(xi)=yi(i=1,2,...,n)

如果三次样条函数：3满足：S(xi)=yi(i=1,2,...,n)例：某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点0:2.5:10(米)，用h表示测量时间0:30:60(秒)，用T表示测试所得各点的温度(℃)。试用线性插值求出在一分钟内每隔20秒、钢轨每隔1米处的温度TI。

命令如下：

x=0:2.5:10;

h=[0:30:60]';

T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];

xi=[0:10];

hi=[0:20:60]';

TI=interp2(x,h,T,xi,hi)例：设z=x2+y2，对z函数在[0,1]×[0,2]区域内进行插值。为了说明这个维数的插值，再考虑一个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。采集了下列的数据：

»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-维相联系(见图6)。为了估计

在中间点的温度，我们必须对它们进行辨识。

»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)%plotresults

»xlabel('WidthofPlate'),ylabel('DegreesCelsius')

»title(['TemperatureatDepth='num2str(d)])图6在深度d=2处的平板温度

另一种方法，我们可以在两个方向插值。先在三维坐标画出原始数据，看一下该数据的粗糙程度(见图7)。

»mesh(width,depth,temps)%usemeshplot

»xlabel(‘WidthofPlate’),ylabel(‘DepthofPlate’)

»zlabel(‘DegreesCelsius’),axis(‘ij’),grid

[xi,yi]=meshgrid(width,depth);

zi=interp2(width,depth,temps,xi,yi,‘cubic’);

mesh(xi,yi,zi)

图7平板温度然后在两个方向上插值，以平滑数据(见图8)。

»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图8二维插值后的平板温度

上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持跟踪。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-轴插值的一个数值数组。

可选的参数method可以是‘linear’，‘cubic’或‘nearest’。在这种情况下，cubic不意味着3次样条，而是使用3次多项式的另一种算法。linear方法是线性插值，仅用作连接图上数据点。nearest方法只选择最接近各估计点的粗略数据点。

插值的优点：利用已知点确定未知点粗糙——精确集合大的——简化的例：某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度(℃)，用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度(℃)。

设时间变量h为一行向量，温度变量t为一个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下：

h=6:2:18;

t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';

XI=6.5:2:17.5

YI=interp1(h,t,XI,‘spline’)

%用3次样条插值计算

在讨论二维插值之前，强调interp1所强制的二个强约束是很重要的。首先，人们不能要求有独立变量范围以外的结果，例如，interp1(hours,temps,13.5)导致一个错误，因为hours在1到12之间变化。其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减小的。

在我们的例子里，hours是单调的。然而，如果我们已经定义独立变量为一天的实际时间，

»time_of_day=[7:121:6]%startat

7AM,endat6PM

time_of_day=

789101112123456

则独立变量将不是单调的，因为time_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours，将会返回一个错误。同样的理由，人们不能对temps插值来找出产生某温度的时间(小时)，因为temps不是单调的。案例3估计水箱的水流量模型长度单位：E(＝30．24cm)容积单位：G(=3.785L(升))