样条插值函数及应用

1、样条插值函数及应用摘要样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值，求取未知函数的近似值。早在一千多年前，我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中，有许多插信函数的曲线要求具有较高的光滑性，在整个曲线中

2、，曲线不但不能有拐点，而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号，这就需要用到三次样条插值。关键词三次样条函数；插值法引言1第一章三次样条插值21.1 样条插值函数简介21.2 三次样条函数应用3第二章AMCM91粘计水塔水流量52.1 理论分析及计算62.2 运用MATLA歆件计算9参考文献14引言样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年

3、Schoenberg将样条引入数学，即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。现在国内外学者对这方面的

4、研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。第一章三次样条插值1.1样条插值函数简介在实际问题中所遇到许多函数f(x)往往很复杂，有些甚至是很难找到解析表达式的。有时通过实验或者数值计算所得到的也只是一些离散的点X(i=0,1,2.n)上的函数值,即yi=f(xi),i=0,1,2.n。根据函数f(x)已有的数据来计算函数f(x)在一些新的点x处的函数值，就是插值法所需要解决的问题。插值法的基本思想就是，首先根据已有的函数值来构造一个简单的函数y(x)作为f(x)的近似表达式，然后用y(x)来

5、计算新的点上的函数值作为f(x)的近似值。通常可以选多项式函数作为近似函数y(x),因为多项式具有各阶导数，求值也比见方便。常用的有Lagrange插值、Newton插值、Hermite插值和样条插值。线性插值在分段点上仅连续而不可导，三次埃尔米特插值有连续的一阶导数，这样的光滑程度常不能满足物理问题的需要，样条函数可以同时解决这两个问题，使插值函数既是低阶分段函数，又是光滑的函数，并且只需在区间端点提供某些导数信息。三次样条函数定义：设在区间a,b】上取n+1个节点a=x0<x1<x2<xn=b,函数y=f(x)在各个节点处的函数值为yi=f(xj,i=0,1,.,n,若S

6、(x)满足：(1)S(x)=yi,i=Q1,.n;(2)在区间Ja,b】上，S(x)具有连续的二阶导数；(3)在区间k/i=0,1,.,n-1)上，S(x)是x三次的多项式；则称S(x)是函数y=f(x)在区间a,b上的三次样条插值函数。由以上定义可以看出，虽然每个子区间上的多项式可以各不相同，但在相邻子区间的连接处却是光滑的。因此，样条插值也称为分段光滑插值。从定义知要求出S(x),在每一个小区间k,xi_Ji=0,1,.,n-1)上确定4个待定系数，共有n个小区间，故应有4n个参数。根据S(x)在口力】上二阶导数连续，在节点x(i=1,2.n-1)出满足连续性条件S(Xi-0)=S(Xi0

7、)S(Xi-0)=S(Xi0).一、，二iS(Xi-0)=S(xi0)共有3n-3个条件，再加上S(x)满足插值条件S(x)=丫4=0,1,.n;共有4n-2个条件,因此还需要2个条件才能确定S(x)o通常可在区间端点上各加一个条件(称为边界条件)，可根据实际问题的要求给定，通常有以下三种：(1)已知端点的一阶导数值，即S(x0)-f0,S(xn)=fn(2)俩端点的二阶导数已知，即S(x0)=f0,S(xn)=fn其特殊情况S"(%)=S"(Xn)=0,称为自然边界条件。(3)当f(x)是以X。为周期的函数时，则要求S(X)也是周期函数。这时边界条件应满足S(Xi-0)=

8、S(Xi0)S(Xi-0)=S(Xi0)S(xi-0)-S(xi0)而此时y°=yn。这样确定的样条函数S(x),称为周期函数。1.2三次样条函数应用作函数y=(x2-3x+7)eXsin(2x)ft10,1取间FB为0.1的点图，用插值进行实验。使用MATLA歆件程序代码如下：%产生原始数据x=0:0.1:1;y=(x.A2-3*x+7).*exp(-4*x).*sin(2*x);%乍图subplot(1,2,1);plot(x,y,x,y,'ro')%寺求插值点xx=0:0,02:1;yy=interp1(x,y,xx,'spline');%乍图s

9、ubplot(1,2,2)plot(x,y,'ro',xx,yy,'b')运行截图运行结果第二章AMCM91AT计水塔水流量美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量，但许多社区并没有测量水流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位，精度在0.5%以内，更为重要的是，无论什么时候,只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两

10、次，每次约两小时。表1白某小镇某天的水塔水位时间033166635106191393717921212402522328543水位317531103054299429472892285027972752时间322843593539332394354331846636499535393657254水位2697水泵工作水泵工作355034453350326031673087时间605746455468535718547502179154826498596889953水位30122927284227672697水泵工作水泵工作34753397时间93270水位3340单位：时间/秒；水位/0.01英

11、尺试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量f(t),并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位达到31.10英尺。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺，通常当水塔的水位降至27.00英尺时水泵开始向水塔充水，而当水塔的水位升至35.50英尺时水泵停止工作。2.1 理论分析及计算1 .水塔充水时间的确定(1)第一次充水时间的确定当时间t=32284秒时，水位26.97英尺，约低于最低水位27英尺，因此可

12、作为第一次开始充水时刻。当t=39435秒时，水塔水位35.5英尺，恰为最高水位，因此可作为第一次充水的结束时刻。充水时间为dt=(39435-32284)/3600=1.9864小时，也接近充水时间2小时。(2)第二次充水时间的确定当时间t=75021秒时，水位26.97英尺，约低于最低水位27英尺，因此可作为第二次开始充水时刻。当t=82649秒时，水泵在工作，但充水时间达到dt=(8264-7502)/3600=2.1189小时；但下一时刻t=85968时，水塔水位34.75英尺，低于最高水位35.50英尺。因此可将t=82649秒作为第二次充水的结束时刻，且该时刻水位为最大充水高度35

13、.50英尺。2 .计算各时刻塔内水的体积单位转换为1英尺=0.3048米，1升=1/3.785411加仑体积计算公式为v=：.d2h/4表2不同时刻水体积表时间水体积时间水体积时间水体积0(1)6061258.967851487219.03755425540.921159371610.954267771519.95945282361.843158302612.032865767020.83925148722.949757157112.954463953422.9581(3)6777153.871456259913.875862235223.88006633974.978155209914.98

14、2260459824.98696485065.900054408115.903958932525.90836376257.006453396316.82615750087.928652537217.9317558781其中(1)表示第一段开始，(2)表示第二段开始，(3)表示第三段开始;单位：时间(小时)，水流量：加仑/小时；3 .计算各时刻点的水流量（加仑/小时）水流量公式为：|dv(t)I出以上25个时刻处的水流量采用差分的方法得到，共分三段分别处理差分公式为：（1）对每段前两点采用向前差分公式f(ti)=-3Vi4Vi-Vi22(tii-1)（2）对每段最后两点采用向后差分公式f(ti)

15、=W-4VVc2(t-t)（3）对每段中间点采用中心差分公式f(ti)=V七+8Vt8/-+V/12(tkti)得到各点水流量表表3不同时刻水流量表时间水流量时间水流量时间水流量0(1)P144048.96781102319.03751166520.92111118210.954(2)1946919.9594144951.84311006312.03282019520.8392146482.9497:1101212.95441894122.958(3)152203.8714879813.87581590323.88001152634.9781999114.98221805524.9869137

16、115.9000:812415.90391564625.908396347.00641016116.8261137427.9286848717.931714962其中（1）表示第一段开始，（2）表示第二段开始，（3）表示第三段开始：单位：时间（小时），水流量：加仑/小时；4.用三次样条拟合流量数据对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线，见图2.1水塔流量图0.80510152025时间（小时）28.B421111苣47中县煤1图2.1水塔流量图其中*表示数据点，实线为样条曲线5.一天总用水量计算一大流水总量计算：方法1:直接积分法：24§

17、;=（f（t）dt=332986（加仑）方法2:分段计算法第一次充水前用水Vi=606125514872=91253（力口仑）第一次充水后第二次充水前用水V2=677715-514872=162843（加仑）22.9581,23.88期间用水V3=677715663397=14318（力口仑）10.9542第一次充水期间用水：V4=f（t）dt=30326（加仑）8.967822.9581第二次充水期间用水：V5=口392f（t）dt=31605（加仑）,2423.88,24期间用水:V6=f（t）dt=1524（加仑）23.886总共用水S2=£Vi=331869（加仑）i4两种方

18、法结果相差err=0.34%6,水泵水流量计算第一次充水期间水塔体积增加V1=677715-514872=162843（力口仑）充水时间：纯=10.9542-8.9678=1.9864（小时）第一次充水期间水泵平均流量10.9542Mf（t）dtP1=678=97246（力1仑/小时）:t1第二次充水期间水塔体积增加忸2=677715-514872=162843（力口仑）充水时间：t2=22.9581-20.8392=2.1189（小时）第二次充水期间水泵平均流量22.9581V2f（t）dtP2=208392-91769（力1仑/小时）上则整个充水期间水泵平均流量P=（P+P2）/2=945

19、07（加仑/小时）2.2运用MATLA歆件计算1,依据理论分析，编写程序代码。MATLAB程序代码c=0.3048;%1英尺等于0.3048米p=1.0/3.785;%1升=1/3,785411力口仑d=57*c;h=31,75*c;v=pi*d*d*h/4*1000*p;data=0,3175;3316,3110;6635,3054;10619,2994;13937,2947;17921,2892;21240,2850;25223,2797;28543,2752;32284,2697;39435,3550;43318,3445;46636,3350;49953,3260;53936,3167

20、;57254,3087;60574,3012;64554,2927;68535,2842;71854,2767;75021,2697;82649,3550;85968,3475;89953,3397;93270,3340;%原始数据t=data(:,1)/3600;%计算时间(小时为单位)v=pi*d*d*data(:,2)/100*c/4*1000*p;%计算体积%计算差分n=length(v);f=zeros(n,1);%存储差分值n1=10;%算第一段fori=1:n1ifi<=2%前两点采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+

21、1)-t(i);elseifi<=n1-2%采用三点中心差分公式f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2)/(12*(t(i+1)-t(i);elseifi>=n1-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1);endendn2=21;%算第二段fori=n1+1:n2ifi<=n1+2%前两点采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i);elseifi<=n2-2f(i)=abs(-v(i+2)+8*v(i+1)-

22、8*v(i-1)+v(i-2)/(12*(t(i+1)-t(i);elseifi>=n2-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1);endendn3=25;%计算第三段fori=n2+1:n3ifi<=n2+2%前两点采用向前差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i);elseifi<=n3-2f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2)/(12*(t(i+1)-t(i);elseifi>=n3-1f(i)=abs(3*

23、v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1);endendplot(t,f,'r*');刎原始点图tmin=min(t);tmax=max(t);tt=tmin:0.1:tmax;%获得离散的时间点，用于作样条曲线ff=spline(t,f,tt);%计算三次样条插值holdonplot(tt,ff,'b');%画样条曲线xlabel(,时间(小时),)；ylabel('流量(加仑/小时),)；title(,水塔流量图,);holdoffdt=0.05;t2=0.5:dt:24.5;%获得离散白时间点，用于积分nn=lengt

24、h(t2);f2=spline(t,f,t2);s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)*dt/2;%计算24小时用水量，采用复化梯形公式fprintf('(全部积分法)1天总水流量s=%8.2fn',s);v10=v(11)-v(10);%第一次水塔增加的水dt1=t(11)-t(10);tp=t(10):dt:t(11);%第一次充水其间流出的水nn=length(tp);yp=spline(t,f,tp);%计算三次样条插值v11=(yp(1)+yp(nn)+2*sum(yp(2:nn-1)*dt/2;v1=v10+v11;%第一次充水总量p1=v1/dt1;%第一次充水的平均水流量v20=v(22)-v(21);%第二次水塔增加的水dt2=t(22)-t(21);tp1=t(21):dt:t(22);%第二次充水其间流出的水nn=length(tp1);yp1=spline(t,f,tp1);%+算三次样条插值v21=(yp1(1)+yp1(nn)+2*sum(yp1(2:nn-1)*dt/2;v2=v20+v21;%!二次充水总量p2=v2/dt2;%第二次充水的平均水流量p=(p1+p2)/2;%两次充水平均水流量fprintf('两次充水平均水流量p=%8.2fn