样条插值函数及应用.doc

样条插值函数及应用摘要样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。关键词 三次样条函数；插值法目 录引 言1第一章三次样条插值21.1 样条插值函数简介21.2 三次样条函数应用3第二章AMCM91A 估计水塔水流量52.1 理论分析及计算62.2运用MATLAB软件计算9参考文献142引 言样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。14第一章三次样条插值1.1 样条插值函数简介 在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。有时通过实验或者数值计算所得到的也只是一些离散的点上的函数值，即。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。插值法的基本思想就是，首先根据已有的函数值来构造一个简单的函数作为的近似表达式，然后用来计算新的点上的函数值作为的近似值。通常可以选多项式函数作为近似函数，因为多项式具有各阶导数，求值也比见方便。常用的有Lagrange插值、Newton插值、Hermite插值和样条插值。线性插值在分段点上仅连续而不可导，三次埃尔米特插值有连续的一阶导数，这样的光滑程度常不能满足物理问题的需要，样条函数可以同时解决这两个问题,使插值函数既是低阶分段函数,又是光滑的函数，并且只需在区间端点提供某些导数信息。三次样条函数定义：设在区间上取个节点，函数在各个节点处的函数值为若满足： (1) (2) 在区间上，具有连续的二阶导数； (3) 在区间上,是三次的多项式；则称是函数在区间上的三次样条插值函数。 由以上定义可以看出，虽然每个子区间上的多项式可以各不相同，但在相邻子区间的连接处却是光滑的。因此，样条插值也称为分段光滑插值。从定义知要求出，在每一个小区间上确定4个待定系数，共有n个小区间，故应有4n个参数。根据在上二阶导数连续，在节点出满足连续性条件共有3n-3个条件，再加上满足插值条件共有4n-2个条件，因此还需要2个条件才能确定。通常可在区间端点上各加一个条件（称为边界条件），可根据实际问题的要求给定，通常有以下三种：（1）已知端点的一阶导数值，即（2）俩端点的二阶导数已知，即其特殊情况称为自然边界条件。（3）当是以为周期的函数时，则要求也是周期函数。这时边界条件应满足而此时。这样确定的样条函数，称为周期函数。1.2 三次样条函数应用作函数在取间隔为0.1的点图，用插值进行实验。使用MATLAB软件程序代码如下：%产生原始数据x=0:0.1:1;y=(x.2-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);%作图subplot(1,2,2)plot(x,y,ro,xx,yy,b) 运行截图图1.1 运行结果第二章AMCM91A 估计水塔水流量美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出当地水塔的水量的设备,他们只能代之以每小时测量水塔中的水位,精度在0.5%以内,更为重要的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时,水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。表1 白某小镇某天的水塔水位时间033166635106191393717921212402522328543水位317531103054299429472892285027972752时间322843593539332394354331846636499535393657254水位2697水泵工作水泵工作355034453350326031673087时间605746455468535718547502179154826498596889953水位30122927284227672697水泵工作水泵工作34753397时间93270水位3340单位：时间/秒；水位/0.01英尺试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量，并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位达到31.10英尺。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺，通常当水塔的水位降至27.00英尺时水泵开始向水塔充水，而当水塔的水位升至35.50英尺时水泵停止工作。2.1 理论分析及计算1. 水塔充水时间的确定(1) 第一次充水时间的确定当时间t=32284秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第一次开始充水时刻。当t=39435秒时，水塔水位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加仑体积计算公式为 表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其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始：单位：时间(小时)，水流量：加仑/小时；4用三次样条拟合流量数据对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图2.1。图2.1 水塔流量图其中*表示数据点，实线为样条曲线5.一天总用水量计算一天流水总量计算:方法1:直接积分法: 方法2:分段计算法第一次充水前用水(加仑) 第一次充水后第二次充水前用水(加仑)22.9581，23.88期间用水(加仑)第一次充水期间用水: 第二次充水期间用水: 23.88，24期间用水: 总共用水两种方法结果相差6.水泵水流量计算第一次充水期间水塔体积增加充水时间:第一次充水期间水泵平均流量第二次充水期间水塔体积增加充水时间:第二次充水期间水泵平均流量则整个充水期间水泵平均流量2.2运用MATLAB软件计算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;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; %计算第一段for i=1:n1 if i=2 %前两点采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i=n1-1 f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endn2=21; %计算第二段for i=n1+1:n2 if i=n1+2 %前两点采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i=n2-1 f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endn3=25; %计算第三段for i=n2+1:n3 if i=n2+2 %前两点采用向前差分 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2)/(2*(t(i+1)-t(i); elseif i=n3-1 f(i)=abs(3*v(i)-4*v(i-1)+v(i-2)/(2*(t(i)-t(i-1); end endplot(t,f,r*); %画原始点图tmin=min(t); tmax=max(t);tt=tmin:0.1:tmax; %获得离散的时间点,用于作样条曲线ff=spline(t,f,tt); %计算三次样条插值hold onplot(tt,ff,b); %画样条曲线xlabel(时间(小时);ylabel(流量(加仑/小时);title(水塔流量图);hold offdt=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)*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,p); vv1=v(1)-v(10