样条插值函数及应用资料

2

摘要

样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。

在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。

插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。

关键词 三次样条函数；插值法

3

目 录

引 言 .................................................. 1 第一章 三次样条插值 .................................... 2 1.1 样条插值函数简介 .................................. 2 1.2 三次样条函数应用 .................................. 3 第二章 AMCM91A 估计水塔水流量 .......................... 5 2.1 理论分析及计算 .................................... 6 2.2运用MATLAB软件计算 ............................... 9 参考文献 ............................................... 14

4

引 言

样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。

在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。

在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。

现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。

1

第一章 三次样条插值

1.1 样条插值函数简介

在实际问题中所遇到许多函数f(x)往往很复杂，有些甚至是很难找到解析表达式的。有时通过实验或者数值计算所得到的也只是一些离散的点xii0，1,2...n上的函数值，即yif(xi),i0,1,2...n。

根据函数f(x)已有的数据来计算函数f(x)在一些新的点x处的函数值，就是插值法所需要解决的问题。

插值法的基本思想就是，首先根据已有的函数值来构造一个简单的函数y(x)作为f(x)的近似表达式，然后用y(x)来计算新的点上的函数值作为

f(x)的近似值。通常可以选多项式函数作为近似函数y(x)，因为多项式具有各阶导数，

求值也比见方便。常用的有Lagrange插值、Newton插值、Hermite插值和样条插值。线性插值在分段点上仅连续而不可导，三次埃尔米特插值有连续的一阶导数，这样的光滑程度常不能满足物理问题的需要，样条函数可以同时解决这两个问题,使插值函数既是低阶分段函数,又是光滑的函数，并且只需在区间端点提供某些导数信息。

三次样条函数定义：

设在区间a,b上取n1个节点ax0x1x2...xnb，函数yf(x)在各个节点处的函数值为yif(xi),i0,1,...,n,若S(x)满足：

(1) S(xi)yi,i0,1,...n;

(2) 在区间a,b上，S(x)具有连续的二阶导数； (3) 在区间xi,xi1上,S(x)是xi0,1,...,n1）（则称S(x)是函数yf(x)在区间a,b上的三次样条插值函数。

由以上定义可以看出，虽然每个子区间上的多项式可以各不相同，但在相邻子区间的连接处却是光滑的。因此，样条插值也称为分段光滑插值。

从定义知要求出S(x)，在每一个小区间xi,xi1上确定4个待定系i0,1,...,n1）（数，共有n个小区间，故应有4n个参数。

根据S(x)在a,b上二阶导数连续，在节点xii1,2...n-1出满足连续性条件

2

S(xi0)S(xi0)

S'(xi0)S'(xi0) S''(xi0)S''(xi0)

共有3n-3个条件，再加上S(x)满足插值条件S(xi)yi,i0,1,...n;共有4n-2个条件，因此还需要2个条件才能确定S(x)。

通常可在区间端点上各加一个条件（称为边界条件），可根据实际问题的要求给定，通常有以下三种：

（1）已知端点的一阶导数值，即

S'(x0)f0,S'(xn)fn'

（2）俩端点的二阶导数已知，即

'S''(x0)f0,S''(xn)fn''

其特殊情况S''(x0)S''(xn)0，称为自然边界条件。

（3）当f(x)是以xnx0为周期的函数时，则要求S(x)也是周期函数。这时边界条件应满足

''S(xi0)S(xi0)

S'(xi0)S'(xi0) S''(xi0)S''(xi0)

而此时y0yn。这样确定的样条函数S(x)，称为周期函数。

1.2 三次样条函数应用

作函数y(x23x7)e4xsin(2x)在0,1取间隔为0.1的点图，用插值进行实验。 使用MATLAB软件 程序代码如下： %产生原始数据 x=0:0.1:1;

y=(x.^2-3*x+7).*exp(-4*x).*sin(2*x);

3

%作图 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%以内,更为重要的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时,水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。

表1 白某小镇某天的水塔水位

时间 0 3316 6635 10619 13937 17921 21240 25223 283 水位 3175 3110 30 2994 2947 22 2850 2797 2752 时间 32284 35935 39332 39435 43318 46636 49953 53936 572 水位 2697 水泵 水泵 3550 3445 3350 3260 3167 3087 工作 工作 时间 60574 5 68535 718 75021 791 829 85968 953 水位 3012 2927 2842 2767 2697 时间 93270 水位 3340 水泵 水泵 3475 3397 工作 工作 单位：时间/秒；水位/0.01英尺 试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量f(t)，并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。

表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位达到31.10英尺。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺，通常当水塔的水位降至27.00英尺时水泵开始向水塔充水，而当水塔的水位升至35.50英尺时水泵停止工作。

5

2.1 理论分析及计算

1. 水塔充水时间的确定 (1) 第一次充水时间的确定

当时间t=32284秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第一次开始充水时刻。

当t=39435秒时，水塔水位35.5英尺，恰为最高水位，因此可作为第一次充水的结束时刻。充水时间为dt=(39435-32284)/3600=1.98小时，也接近充水时间2小时。

(2) 第二次充水时间的确定

当时间t=75021秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第二次开始充水时刻。

当t=829秒时，水泵在工作，但充水时间达到dt=(82-7502)/3600=2.11小时；但下一时刻t=85968时，水塔水位34.75英尺，低于最高水位35.50 英尺。

因此可将t=829秒作为第二次充水的结束时刻,且该时刻水位为最大充水高度35.50 英尺。

2.计算各时刻塔内水的体积

单位转换为1英尺=0.3048米, 1升=1/3.7811加仑 体积计算公式为v.d2h/4

表2 不同时刻水体积表

时间 0 (1) 0.9211 1.8431 2.9497 3.8714 4.9781 5.9000 7.00 7.9286

其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始； 单位：时间(小时)，水流量：加仑/小时；

6

水体积 606125 593716 583026 571571 562599 552099 4081 533963 525372 时间 8.9678 10.92 (2) 12.0328 12.94 13.8758 14.9822 15.9039 16.8261 17.9317 水体积 时间 514872 19.0375 677715 19.9594 657670 20.8392 639534 22.9581 (3) 622352 23.8800 604598 24.9869 5325 25.9083 575008 558781 水体积 25 528236 514872 677715 663397 8506 637625 3.计算各时刻点的水流量(加仑/小时) 水流量公式为：

f(t)dv(t) dt以上25个时刻处的水流量采用差分的方法得到，共分三段分别处理。 差分公式为：

(1) 对每段前两点采用向前差分公式

f(ti)3Vi4Vi1Vi2

2(ti1ti)(2) 对每段最后两点采用向后差分公式

f(ti)3Vi4Vi1Vi2

2(titi1)(3) 对每段中间点采用中心差分公式

f(ti)得到各点水流量表

Vi28Vi18Vi1Vi2

12(ti1ti)表3 不同时刻水流量表

时间 0 (1) 0.9211 1.8431 2.9497 3.8714 4.9781 5.9000 7.00 7.9286 水流量 14404 11182 10063 11012 8798 9991 8124 10161 8487 时间 8.9678 10.9(2) 12.0328 12.94 13.8758 14.9822 15.9039 16.8261 17.9317 水流量 11023 19469 20195 141 15903 18055 156 13742 14962 时间 19.0375 19.9594 20.8392 22.958(3) 23.8800 24.9869 25.9083 水流量 16652 14495 148 15220 15263 13711 9634 其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始： 单位：时间(小时)，水流量：加仑/小时； 4．用三次样条拟合流量数据

对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图2.1。

7

图2.1 水塔流量图

其中*表示数据点，实线为样条曲线 5.一天总用水量计算 一天流水总量计算: 方法1:直接积分法:

24S1f(t)dt332986(加仑)

0方法2:分段计算法

第一次充水前用水V160612551487291253(加仑)

第一次充水后第二次充水前用水V2677715(加仑) 514872162843[22.9581，23.88]期间用水V367771566339714318(加仑)

第一次充水期间用水: V410.928.967822.9581f(t)dt30326(加仑) f(t)dt31605(加仑) f(t)dt1524(加仑)

8

第二次充水期间用水: V520.839224[23.88，24]期间用水: V6

23.88总共用水S2Vi331869(加仑)

i16两种方法结果相差err6.水泵水流量计算

S1S20.34% S1第一次充水期间水塔体积增加

V1677715514872162843(加仑)

充水时间:

t110.928.96781.98(小时)

第一次充水期间水泵平均流量

p1V110.928.9678f(t)dtt197246(加仑/小时)

第二次充水期间水塔体积增加

V2677715514872162843(加仑)

充水时间:

t222.958120.83922.11(小时)

第二次充水期间水泵平均流量

p2V222.958120.8392f(t)dtt291769(加仑/小时)

则整个充水期间水泵平均流量

p(p1p2)/294507(加仑/小时)

2.2运用MATLAB软件计算

1.依据理论分析，编写程序代码。 MATLAB程序代码

c=0.3048; %1英尺等于0.3048米 p=1.0/3.785; %1升=1/3.7811加仑 d=57*c; h=31.75*c;

9

v=pi*d*d*h/4*1000*p;

data=[0,3175;3316,3110;6635,30;10619,2994;13937,2947;17921,22; 21240,2850;25223,2797;283,2752;32284,2697;39435,3550; 43318,3445;46636,3350;49953,3260;53936,3167;572,3087; 60574,3012;5,2927;68535,2842;718,2767;75021,2697; 829,3550;85968,3475; 953,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-2 %采用三点中心差分公式

f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(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 end

n2=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-2

f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));

10

elseif i>=n2-1

f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1))); end end

n3=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-2

f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(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 end

plot(t,f,'r*'); %画原始点图 tmin=min(t); tmax=max(t);

tt=tmin:0.1:tmax; %获得离散的时间点,用于作样条曲线 ff=spline(t,f,tt); %计算三次样条插值 hold on

plot(tt,ff,'b'); %画样条曲线 xlabel('时间(小时)'); ylabel('流量(加仑/小时)'); title('水塔流量图'); hold off dt=0.05;

t2=0.5:dt:24.5; %获得离散的时间点,用于积分 nn=length(t2); f2=spline(t,f,t2);

11

s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)))*dt/2 ;%计算24小时用水量,采用复化梯形公式 fprintf('(全部积分法)1天总水流量s= %8.2f\\n',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.2f\\n',p);

vv1=v(1)-v(10); %第一次充水前总流量 vv2=v(11)-v(21); %两次充水间总流量 vv3=v(22)-v(23); %t=839到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到800流量

tc=t(23):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.2f\\n',ss);

%误差分析 err=abs((s-ss)/s);

fprintf('两种计算法总量相对误差%6.2f%%\\n',err*100);

2.运行代码，得出结果。 运行结果如下

图2.2 运行结果

13

参考文献

[1] 施吉林,刘淑珍.计算机数值方法[M].北京:高等教育出版社，1999. [2] 常庚哲.关于三次样条函数的两点注记[J].数学的实践与认识，1979. [3] 林成森编著.数值计算方法[M].科学出版社,2005. [4] 李庆扬著.数值分析[M].清华大学出版社,2010.

14