样条插值实验与Matlab仿真.doc

样条插值实验
Matlab的一维插值函数为interp1()，调用格式为：
yy=interp1(x,y,xx,方法)
其中x=[x1,x2,…,xn]’, y=[y1,y2,…,yn]’,两个向量分别为给定的一组自变量和函数值，用来表示已知样本点数据.xx为待求插值点处横坐标，yy返回的对应纵坐标。
插值方法可以选用默认的’linear’(线性插值)，’nearest’(最近邻等值方式)，’cubic’（三次Hermite插值），’spline’(三次样条插值)。一般可采用样条插值。

实验1
作函数在[0,1]取间隔为0.1的点图，用插值进行实验
Matlab程序chazhi1.m：
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'); %此处可用nearest,cubic,spline分别试验

subplot(1,2,2)
plot(x,y,'ro',xx,yy,'b') %作图

图1(1) 原始点图 图1(2) 样条插值图

实验2. AMCM91A 估计水塔水流量
美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出当地水塔的水量的设备,他们只能代之以每小时测量水塔中的水位,精度在0.5%以内,更为重要的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时,水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。

表1 白某小镇某天的水塔水位
时间(秒)
0
3316
6635
10619
13937
17921
21240
25223
28543
水位(0.01
英尺)
3175
3110
3054
2994
2947
2892
2850
2797
2752
时间(秒)
32284
35935
39332
39435
43318
46636
49953
53936
57254
水位
2697
水泵
工作
水泵
工作
3550
3445
3350
3260
3167
3087
时间
60574
64554
68535
71854
75021
79154
82649
85968
89953
水位
3012
2927
2842
2767
2697
水泵
工作
水泵
工作
3475
3397
时间
93270

水位
3340

试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量，并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。
表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位达到31.10英尺。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺，通常当水塔的水位降至27.00英尺时水泵开始向水塔充水，而当水塔的水位升至35.50英尺时水泵停止工作。

解答:
1. 水塔充水时间的确定
(1) 第一次充水时间的确定
当时间t=32284秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第一次开始充水时刻。
当t=39435秒时，水塔水位35.5英尺，恰为最高水位，因此可作为第一次充水的结束时刻。充水时间为dt=(39435-32284)/3600=1.9864小时，也接近充水时间2小时。
(2) 第二次充水时间的确定
当时间t=75021秒时，水位26.97英尺,约低于最低水位27英尺，因此可作为第二次开始充水时刻。
当t=82649秒时，水泵在工作，但充水时间达到dt=(82649-75021)/3600=2.1189小时；但下一时刻t=85968时，水塔水位34.75英尺，低于最高水位35.50 英尺。
因此可将t=82649秒作为第二次充水的结束时刻,且该时刻水位为最大充水高度35.50 英尺。
2.计算各时刻塔内水的体积
单位转换为1英尺=0.3048米, 1升=1/3.785411加仑
体积计算公式为

表2 不同时刻水体积表
时间(小时)
水体积(加仑)
时间(小时)
水体积(加仑)
时间(小时)
水体积(加仑)
0 (1)
606125
10.9542 (2)
677715
20.8392
514872
0.9211
593716
12.0328
657670
22.9581 (3)
677715
1.8431
583026
12.9544
639534
23.8800
663397
2.9497
571571
13.8758
622352
24.9869
648506
3.8714
562599
14.9822
604598
25.9083
637625
4.9781
552099
15.9039
589325

5.9000
544081
16.8261
575008

7.0064
533963
17.9317
558781

7.9286
525372
19.0375
542554

8.9678
514872
19.9594
528236

其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始

3.计算各时刻点的水流量(加仑/小时)
水流量公式为：
以上25个时刻处的水流量采用差分的方法得到，共分三段分别处理。
差分公式为：
(1) 对每段前两点采用向前差分公式

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

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

得到各点水流量表

表3 不同时刻水流量表
时间(小时)
水流量
(加仑/小时)
时间(小时)
水流量
(加仑/小时)
时间(小时)
水流量
(加仑/小时)
0 (1)
14404
10.9542 (2)
19469
20.8392
14648
0.9211
11182
12.0328
20195
22.9581 (3)
15220
1.8431
10063
12.9544
18941
23.8800
15263
2.9497
11012
13.8758
15903
24.9869
13711
3.8714
8798
14.9822
18055
25.9083
9634
4.9781
9991
15.9039
15646

5.9000
8124
16.8261
13742

7.0064
10161
17.9317
14962

7.9286
8487
19.0375
16652

8.9678
11023
19.9594
14495

其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始
4．用三次样条拟合流量数据
对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图1。

图1 水塔流量图
其中*表示数据点，实线为样条曲线

5.一天总用水量计算
一天流水总量计算:
方法1:直接积分法:

方法2:分段计算法
第一次充水前用水(加仑)
第一次充水后第二次充水前用水(加仑)
[22.9581，23.88]期间用水(加仑)
第一次充水期间用水:
第二次充水期间用水:
[23.88，24]期间用水:
总共用水
两种方法结果相差
6.水泵水流量计算
第一次充水期间水塔体积增加

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

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

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

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

附Matlab程序
%AMCM91A
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-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)));
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);
%计算24小时用水量,采用复化梯形公式
s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)))*dt/2 ;
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);
% t=83649到85968其间流量
vv3=v(22)-v(23);
%第一次充水期间流量
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):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);