样条插值函数及应用.docx

样条插值函数及应用摘要样条函数具有广泛的应用，是现代函数论的一个十分活泼的分支，是计算方法的主要根底和工具之一，由于生产和科学技术向前开展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。插值法是数值逼近的重耍方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中，有许多插值函数的曲线要求具有较高的光滑性，在整个曲线中，曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号，这就需要用到三次样条插值。关键词三次样条函数；插值法目录引言O第一章三次样条插值11.1 样条插值函数简介11.2 三次样条函数应用2第二章AMCM91A估计水塔水流量32.1理论分析及计算42.2运用MATLAB软件计算6参考文献10引言样条函数具有广泛的应用，是现代函数论的一个十分活泼的分支，是计算方法的主要根底和工具之一，由于生产和科学技术向前开展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代,在研究数据处理的问题中引出了样条函数，例如，在1946年SChoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改良，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。第一章三次样条插值11样条插值函数简介在实际问题中所遇到许多函数/。)往往很复杂，有些甚至是很难找到解析表达式的。有时通过实验或者数值计算所得到的也只是一些离散的点NG=(),12M上的函数值，即yf-=/(xl),/=0,1,2.wo根据函数/(x)已有的数据来计算函数八幻在一些新的点X处的函数值，就是插值法所需要解决的问题。插值法的根本思想就是，首先根据已有的函数值来构造一个简单的函数y(x)作为/(尢)的近似表达式，然后用y(x)来计算新的点上的函数值作为/a)的近似值。通常可以选多项式函数作为近似函数y(),因为多项式具有各阶导数，求值也比见方便。常用的有Lagrange插值、NeWtOn插值、HelTnite插值和样条插值。线性插值在分段点上仅连续而不可导，三次埃尔米特插值有连续的一阶导数，这样的光滑程度常不能满足物理问题的需要，样条函数可以同时解决这两个问题,使插值函数既是低阶分段函数，又是光滑的函数，并且只需在区间端点提供某些导数信息。三次样条函数定义：设在区间a,b上取n+1个节点。=<x1<x2.<xtl=bf函数y=/(x)在各个节点处的函数值为Yi=/(xf),z=0,1,.,小假设S(X)满足：(1)S(Xi)=y,i=0,1,.；(2)在区间a,“上，S(X)具有连续的二阶导数；(3)在区间IXjk=O,上，S(X)是X三次的多项式；那么称Sa)是函数y=/(x)在区间a,"上的三次样条插值函数。由以上定义可以看出，虽然每个子区间上的多项式可以各不相同，但在相邻子区间的连接处却是光滑的。因此，样条插值也称为分段光滑插值。从定义知要求出S(R),在每一个小区间x”XzH=0,l,.,-1)上确定4个待定系数，共有n个小区间，故应有4n个参数。根据S(X)在a,"上二阶导数连续，在节点Ji=l,2.-1)出满足连续性条件共有3n-3个条件，再加上S(X)满足插值条件S(Xj)=y/=0,1,.;共有4n-2个条件,因此还需要2个条件才能确定S(X)。通常可在区间端点上各加一个条件(称为边界条件，可根据实际问题的要求给定，通常有以下三种：(1)端点的一阶导数值，即(2)俩端点的二阶导数，即其特殊情况S"(x0)=S(x,r)=0,称为自然边界条件。(3)当/(幻是以居-%为周期的函数时，那么要求Sa)也是周期函数。这时边界条件应满足而此时光=必。这样确定的样条函数S(x),称为周期函数。1.2三次样条函数应用作函数y=(2-3x+7)xsin(2x)在()1取间隔为0.1的点图,用插值进行实验。使用MATLAB软件程序代码如下：%产生原始数据x=0:0.1:1;y=(x.2-3*x+7).*exp(-4*x).*sin(2*x);%作图subplot(1,2,1);PIot(X,y,x,y,To')%待求插值点xx=0:0.02:l;yy=interpl(x,y,xx,'spline');IM乍图subplot(1,2,2)plot(x,y,'ro',xx,yy,b,)运行截图图Ll运行结果第二章AMCM91A估计水塔水流量美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位,精度在0.5%以内，更为重耍的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时,水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。表1白某小镇某天的水塔水位时间033166635106191393717921212402522328543水位317531103054299429472892285027972752时间322843593539332394354331846636499535393657254水位2697水泵工作水泵工作355034453350326031673087时间605746455468535718547502179154826498596889953水位30122927284227672697水泵工作水泵工作34753397时间93270水位3340单位：时间/秒；水位001英尺试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量/，并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位到达31.10英尺。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺，通常当水塔的水位降至27.00英尺时水泵开始向水塔充水，而当水塔的水位升至35.50英尺时水泵停止工作。2.1理论分析及计算1.水塔充水时间确实定（1）第一次充水时间确实定当时间t=32284秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第一次开始充水时刻。当439435秒时，水塔水位35.5英尺，恰为最高水位，因此可作为第一次充水的结束时刻。充水时间为力=（39435-32284）/3600=1.9864小时，也接近充水时间2小时。（2）第二次充水时间确实定当时间t=75021秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第二次开始充水时刻。当t=82649秒时，水泵在工作，但充水时间到达力二（8264-7502）/3600=2.1189小时；但下一时刻t=85968时，水塔水位34.75英尺，低于最高水位35.50英尺。因此可将t=82649秒作为第二次充水的结束时刻,且该时刻水位为最大充水高度35.50英尺。2.计算各时刻塔内水的体积单位转换为1英尺二0.3048米，1升=1/3.785411加仑体积计算公式为U=加筋/4表2不同时刻水体积表时间水体积时间水体积时间水体积0(1)6061258.967851487219.03755425540.921159371610.9542(2)67771519.95945282361.843158302612.032865767020.83925148722.949757157112.954463953422.9581(3)6777153.871456259913.875862235223.88006633974.978155209914.982260459824.98696485065.900054408115.903958932525.90836376257.006453396316.82615750087.928652537217.9317558781其中（1）表示第一段开始，（2）表示第二段开始，（3）表示第三段开始;单位：时间（小时），水流量：加仑/小时；3.计算各时刻点的水流量（加仑/小时）水流量公式为：以上25个时刻处的水流量采用差分的方法得到，共分三段分别处理。差分公式为：（1）对每段前两点采用向前差分公式（2）对每段最后两点采用向后差分公式（3）对每段中间点采用中心差分公式得到各点水流量表表3不同时刻水流量表时间水流量时间水流量时间水流量0(1)144048.96781102319.0375166520.92111118210,954(2)1946919.9594144951.84311006312.03282019520.8392146482.94971101212.95441894122.958(3)152203.8714879813.87581590323.8800152634.9781999114.98221805524.9869137115.9000812415.90391564625.908396347.00641016116.8261137427.9286848717.931714962其中（I）表示第一段开始，（2）表示第二段开始，（3）表示第三段开始：单位：时间（小时），水流量：加仑/小时；4.用三次样条拟合流量数据对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图2.1。图2.1水塔流量图其中*表示数据点，实线为样条曲线5.一天总用水量计算一天流水总量计算：方法1:直接积分法：方法2:分段计算法第一次充水前用水K=606125-514872=91253（加仑）第一次充水后第二次充水前用水=677715-514872=162843（加仑）22.9581,23.88期间用水匕=677715-663397=14318（加仑）Clo9542第一次充水期间用水：V4=j9678/（r）Jr=3（）326（加仑）,2）9581第二次充水期间用水：V5=Q）M=31605（力口仑）'J20.839223.88,24期间用水：V6=/Q）d=1524（加仑）oJ23.886总共用水邑=WX=331869（加仑）i=l两种方法结果相差e"=生邑=0.34%Sl6.水泵水流量计算第一次充水期间水塔体积增加充水时间：第一次充水期间水泵平均流量第二次充水期间水塔体积增加充水时间：第二次充水期间水泵平均流量那么整个充水期间水泵平均流量2.2运用MATLAB软件计算1.依据理论分析，编写程序代码。MATLAB程序代码c=0.3048;%l英尺等于0.3048米p=L0/3.785;%l升=1/3.785411加仑d=57*c;h=31.75*c;v=pi*d*d*h4*1000*p;data=l,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;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)/l00*c4*1Oo0*p;%计算体积%计算差分n=length(v);f=zeros(nj);/存储差分值nl=10;%计算第一段fori=l:nlifi<=2%前两点采用向前差分f(i)=abs(-3*v(i)+4*v(i+l)-v(i+2)(2*(t(i+l)-t(i);elseifi<=n12%采用三点中心差分公式f(i)=abs(-v(i÷2)+8*v(i+l)-8*v(i-l)+v(i-2)(l2*(t(i+l)-t(i);elseifi>=nl-lf(i)=abs(3*v(i)-4*v(i-1)+v(i-2)(2*(t(i)-t(i-1);endendn2=21;%计算第二段fori=nl+kn2ifi<=nl+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)-8*v(i-1)+v(i-2)(12*(t(i+l)-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+kn3ifi<=n2+2胴前两点采用向前差分f(i)=abs(-3*v(i)÷4*v(i+l)-v(i+2)(2*(t(i+l)-t(i);elseifi<=n3-2f(i)=abs(-v(i+2)+8*v(i+l)-8*v(i-l)+v(i-2)(l2*(t(i+l)-t(i);elseifi>=n3-lf(i)=abs(3*v(i)-4*v(i-1)+v(i-2)(2*(t(i)-t(i-1);endendPIOt(Lf画原始点图tmin=min(t);tmax=max(t);tt=tmin:0.1:tmax;%获得离散的时间点,用于作样条曲线ff=spline(t,f,tt);%计算三次样条插值holdonplot(tt,ff,b)%画样条曲线XlabelC时间(小时),);ylabM流量(加仑/小时);MtleC水塔流量图上holdoffdt=0.05;t2=0.5:dt:24.5;%获得离散的时间点,用于积分nn=length(t2);f2=spline(t,f,t2);s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)*dt2;%计算24小时用水量,采用复化梯形公式fprintfC(全部积分法)1天总水流量s=%8.2fn,s);V10=v(ll)-v(l0);%第一次水塔增加的水dtl=t(ll)-t(10);tp=t(10):dt:t(ll);%第一次充水其间流出的水nn=length(tp);yp=spline(t,f,tp)也计算三次样条插值V11=(yp(1)+yp(nn)+2*sum(yp(2:nn-1)*dt2;Vl=VlO+vll;%第一次充水总量pl=vldtl制第一次充水的平均水流量v20=v(22)v);%第二次水塔增加的水dt2=t(22)-t(21);tpl=t(21):dt:t(22);%第二次充水其间流出的水nn=length(tpl);yl=sline(t,f,tl);%计算三次样条插值v21=(ypl(l)+ypl(nn)÷2*sum(yp1(2:nn-l)*dt/2;v2=v20+v21;*第二次充水总量p2=v2dt2;%第二次充水的平均水流量p=(pl+p2)2%两次充水平均水流量fprintfC两次充水平均水流量p=%8.2fn,p);VVl=V(I第一次充水前总流量vv2=v(ll)-v(21)*两次充水间总流量vv3=v(22)-v(23);%t=83649到85968其间流量%第一次充水期间流量ta=t(10):dt:t(ll);%获得离散的时间点,用于积分nn=length(ta);fa=spline(t,f,ta);si=(fa(1)+fa(nn)+2*sum(fa(2:nn-1)*dt2;(第二次充水期间流量tb=t(21):dt:t(22);%获得离散的时间点,用于积分nn=length(tb);fb=spline(t,f,tb);s2=(fb(1)+fb(nn)+2*sum(fb(2:nn-1)*dt2;%t=85968到86400流量tc=t(23):dt:24;国获得离散的时间点,用于积分nn=length(tc);fc=spline(t,f,tc);s3=(fc(1)+fc(nn)+2*sum(fc(2:nn-1)*dt2;ss=vvl+vv2+vv3+s1+s2+s3;fprintf('(局部积分法)1天总水流量ss=%8.2fn,ss);国误差分析err=abs(s-ss)s);fprintf(,两种计算法总量相对误差%6.2f%nerr*100);2.运行代码，得出结果。运行结果如下图2.2运行结果参考文献1施吉林,刘淑珍.计算机数值方法M.北京:高等教育出版社，1999.2常庚哲.关于三次样条函数的两点注记J.数学的实践与认识，1979.3林成森编著.数值计算方法M.科学出版社,2005.4李庆扬著.数值分析M.清华大学出版社,2010.