1、0样 条 插 值 函 数 及 应 用1摘要样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的
2、光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。关键词 三 次 样 条 函 数 ； 插 值 法2目 录引 言 .1第一章 三次样条插值 .21.1 样条插值函数简介 21.2 三次样条函数应用 3第二章 AMCM91A 估计水塔水流量 .52.1 理论分析及计算 62.2 运用 MATLAB 软件计算 9参考文献 140引 言样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深
3、刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在 1946 年 Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地
4、逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象” 。现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。1第一章 三次样条插值1.1 样条插值函数简介 在实际问题中所遇到许多函数 往往很复杂，有些甚至是很难找到解析表达)(xf式的。有时通过实验或者数值计算所得到的也只是一些离散的点 上的nix.2,10，函数值，即 。nixfi .2,10),(y根据函数 已有的数据来计算函数 在一些新的点 处的函数值，就是插)(xf值法所需要解决的问题。插 值 法 的
5、 基 本 思 想 就 是 ， 首 先 根 据 已 有 的 函 数 值 来 构 造 一 个 简 单 的 函 数作 为 的 近 似 表 达 式 ， 然 后用 来计算新的点上的函数值作为 的)(yx)(f )(yx )(xf近 似 值 。 通 常 可 以 选 多 项 式 函 数 作 为 近 似 函 数 ，因为多项式具有各阶导数，求值也比见方便。常用的有 Lagrange 插值、Newton 插值、Hermite 插 值 和 样 条 插 值 。 线性 插 值 在 分 段 点 上 仅 连 续 而 不 可 导 ， 三 次 埃 尔 米 特 插 值 有 连 续 的 一 阶 导 数 ， 这 样的 光 滑 程 度
6、 常 不 能 满 足 物 理 问 题 的 需 要 ， 样 条 函 数 可 以 同 时 解 决 这 两 个 问 题 ,使插 值 函 数 既 是 低 阶 分 段 函 数 ,又 是 光 滑 的 函 数 ， 并 且只需在区间端点提供某些导数信息。三次样条函数定义：设在区间 上取 个节点 ，函数 在 各 个 节b,a1nbxxan.210 )(yxf点 处 的 函 数 值 为 若 满足： ,.),(yiixf)(S(1) ;.0,Sxii(2) 在区间 上， 具 有 连 续 的 二 阶 导 数 ； ba)(3) 在区间 上, 是 三 次 的 多 项 式 ； ）（ 1,.,1nii )(Sx则称 是函数
7、在区间 上的 三 次 样 条 插 值 函 数 。Sx)(yxfba由以上定义可以看出，虽然每个子区间上的多项式可以各不相同，但在相邻子区间的连接处却是光滑的。因此，样条插值也称为分段光滑插值。从定义知要求出 ，在每一个小区间 上确定4个待定系)(Sx）（ 1,.0x,1nii数，共有 n 个小区间，故应有4n 个参数。根据 在 上二阶导数连续，在节点 出满足连续性条件)(xb,a-.2,i2)0()(Siix ii )()( iixx共有3n-3个条件，再加上 满足插值条件 共有4n-2个条件，)S;,.10,ySnii因此还需要2个条件才能确定 。(x通常可在区间端点上各加一个条件（称为边界
8、条件） ，可根据实际问题的要求给定，通常有以下三种：（1）已知端点的一阶导数值，即 0 )(S,)(nfxfx（2）俩端点的二阶导数已知，即 0 )(,)(nff其特殊情况 称为自然边界条件。，)(S0 nx（3）当 是以 为周期的函数时，则要求 也是周期函数。这时边)(xf0)(Sx界条件应满足 )0()(iixSii )()( iixx而此时 。这样确定的样条函数 ，称为周期函数。ny01.2 三次样条函数应用作函数 在 取间隔为 0.1 的点图，用插值进行实)2sin()73(42xexyx1,0验。使用 MATLAB 软件程序代码如下：%产生原始数据x=0:0.1:1;3y=(x.2-
9、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 运行结果4第二章 AMCM91A 估计水塔水流量美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出当地水塔的水量的设备,他们只能代之以每小时测量水塔中的水位,精度在 0.5%以内,更为重要的是,无论什么时候,只要水塔中的水
10、位下降到某一最低水位 L 时,水泵就启动向水塔重新充水至某一最高水位 H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。表 1 白某小镇某天的水塔水位时间 0 3316 6635 10619 13937 17921 21240 25223 28543水位 3175 3110 3054 2994 2947 2892 2850 2797 2752时间 32284 35935 39332 39435 43318 46636 49953 53936 57254水位 2697 水泵工作 水泵工作
11、 3550 3445 3350 3260 3167 3087时间 60574 64554 68535 71854 75021 79154 82649 85968 89953水位 3012 2927 2842 2767 2697 水泵工作 水泵工作 3475 3397时间 93270水位 3340单位：时间/秒；水位/0.01 英尺试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量 ，)(tf并估计一天的总用水量，表 1 中给出了某个真实小镇某一天的真实数据。表 1 中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316
12、秒后，水塔中的水位达到31.10 英尺。水塔是一个垂直圆形柱体，高为 40 英尺，直径 57 英尺，通常当水塔的水位降至 27.00 英尺时水泵开始向水塔充水，而当水塔的水位升至 35.50 英尺时水泵停止工作。52.1 理论分析及计算1. 水塔充水时间的确定(1) 第一次充水时间的确定当时间 t=32284 秒时，水位 26.97 英尺,约低于最低水位 27 英尺，因此可作为第一次开始充水时刻。当 t=39435 秒时，水塔水位 35.5 英尺，恰为最高水位，因此可作为第一次充水的结束时刻。充水时间为 =(39435-32284)/3600=1.9864 小时，也接近充水时间 2dt小时。(
13、2) 第二次充水时间的确定当时间 t=75021 秒时，水位 26.97 英尺,约低于最低水位 27 英尺，因此可作为第二次开始充水时刻。当 t=82649 秒时，水泵在工作，但充水时间达到 =(8264-7502)/3600=2.1189dt小时；但下一时刻 t=85968 时，水塔水位 34.75 英尺，低于最高水位 35.50 英尺。因此可将 t=82649 秒作为第二次充水的结束时刻,且该时刻水位为最大充水高度35.50 英尺。2.计算各时刻塔内水的体积单位转换为 1 英尺=0.3048 米, 1 升=1/3.785411 加仑体积计算公式为 4/.2hdv表 2 不同时刻水体积表其中
14、(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始；单位：时间(小时)，水流量：加仑/小时；时间 水体积 时间 水体积 时间 水体积0 (1) 606125 8.9678 514872 19.0375 5425540.9211 593716 10.9542 (2) 677715 19.9594 5282361.8431 583026 12.0328 657670 20.8392 5148722.9497 571571 12.9544 639534 22.9581 (3) 6777153.8714 562599 13.8758 622352 23.8800 6633974.97
15、81 552099 14.9822 604598 24.9869 6485065.9000 544081 15.9039 589325 25.9083 6376257.0064 533963 16.8261 5750087.9286 525372 17.9317 55878163.计算各时刻点的水流量(加仑/小时)水流量公式为：dtvtf)(以上 25 个时刻处的水流量采用差分的方法得到，共分三段分别处理。差分公式为：(1) 对每段前两点采用向前差分公式 )(243)(12iii tVtf(2) 对每段最后两点采用向后差分公式 1234()()iiiiftt(3) 对每段中间点采用中心差分公式
16、 )(128)( 21iiiii tVVtf得到各点水流量表表 3 不同时刻水流量表其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始：单位：时间(小时)，水流量：加仑/小时；4用三次样条拟合流量数据对表 3 中 25 个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图 2.1。时间 水流量 时间 水流量 时间 水流量0 (1) 14404 8.9678 11023 19.0375 166520.9211 11182 10.954(2) 19469 19.9594 144951.8431 10063 12.0328 20195 20.8392
17、 146482.9497 11012 12.9544 18941 22.958(3) 152203.8714 8798 13.8758 15903 23.8800 152634.9781 9991 14.9822 18055 24.9869 137115.9000 8124 15.9039 15646 25.9083 96347.0064 10161 16.8261 137427.9286 8487 17.9317 149627图 2.1 水塔流量图其中*表示数据点，实线为样条曲线5.一天总用水量计算一天流水总量计算:方法 1:直接积分法: 2401 )32986()(加 仑dtfS方法 2:
18、分段计算法第一次充水前用水 (加仑) 155761V第一次充水后第二次充水前用水 (加仑)16284376222.9581，23.88期间用水 (加仑)393第一次充水期间用水: 9542.106784 )026()(加 仑dtfV第二次充水期间用水: .3920515加 仑23.88，24期间用水: 48.6 )(4)(加 仑tf8总共用水 )(31869612加 仑iVS两种方法结果相差 %34.012Ser6.水泵水流量计算第一次充水期间水塔体积增加 )(162843576711 加 仑V充水时间: )(9942.01 小 时t第一次充水期间水泵平均流量 )/(746)(19542.06
19、781 小 时加 仑tdtfVp第二次充水期间水塔体积增加 )(128352 加 仑充水时间: )(90981.2 小 时t第二次充水期间水泵平均流量 )/(76)(25.83902 小 时加 仑tdtfVp则整个充水期间水泵平均流量 )/(450/)(21 小 时加 仑p2.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;9v=pi*d*d*h/4*1000*p;data=0,3175;3316,3110;66
20、35,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*dat
21、a(:,2)/100*c/4*1000*p; %计算体积%计算差分n=length(v);f=zeros(n,1); %存储差分值n1=10; %计算第一段for i=1:n1if i=n1-1f(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:n2if i=n2-1f(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:n3if i=n3-1f(i)=abs(3*v(i)-4*v(i-1
22、)+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);11s=(f2(1
23、)+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
24、(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); %第一次充水前总流量vv2=v(11)-v(21); %两次充水间总
25、流量vv3=v(22)-v(23); %t=83649 到 85968 其间流量12%第一次充水期间流量ta=t(10):dt:t(11); %获得离散的时间点,用于积分nn=length(ta);fa=spline(t,f,ta);s1=(fa(1)+fa(nn)+2*sum(fa(2:nn-1)*dt/2 ;%第二次充水期间流量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)*dt/2 ;%t=85968 到 86400 流量tc=t(23):
26、dt:24; %获得离散的时间点,用于积分nn=length(tc);fc=spline(t,f,tc);s3=(fc(1)+fc(nn)+2*sum(fc(2:nn-1)*dt/2 ;ss=vv1+vv2+vv3+s1+s2+s3;fprintf(部分积分法)1 天总水流量 ss= %8.2fn,ss);%误差分析err=abs(s-ss)/s);fprintf(两种计算法总量相对误差 %6.2f%n,err*100);2.运行代码，得出结果。运行结果如下13图2.2 运行结果14参考文献1 施吉林,刘淑珍.计算机数值方法M.北京:高等教育出版社，1999.2 常庚哲.关于三次样条函数的两点注记J.数学的实践与认识，1979.3 林成森编著.数值计算方法M.科学出版社,2005.4 李庆扬著.数值分析M.清华大学出版社,2010.
