《样条插值实验与Matlab仿真.doc》由会员分享，可在线阅读，更多相关《样条插值实验与Matlab仿真.doc（9页珍藏版）》请在三一文库上搜索。
1、 样条插值实验 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).*e
2、xp(-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 估计水塔水流量 美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出当地水塔的水量的设备,他们只能代之以每小时测
3、量水塔中的水位,精度在0.5%以内,更为重要的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时,水泵就启动向水塔重新充水至某一最高水位H,但也无法得到水泵的供水量的测量数据。因此，在水泵工作时，人们容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约两小时。表1 白某小镇某天的水塔水位时间(秒)033166635106191393717921212402522328543水位(0.01 英尺)317531103054299429472892285027972752时间(秒)322843593539332394354331846636499535393
4、657254水位2697水泵工作水泵工作355034453350326031673087时间605746455468535718547502179154826498596889953水位30122927284227672697水泵工作水泵工作34753397时间93270水位3340试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量，并估计一天的总用水量，表1中给出了某个真实小镇某一天的真实数据。表1中给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，3316秒后，水塔中的水位达到31.10英尺。水塔是一个垂直圆形柱体，高为4
5、0英尺，直径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英尺，因此可作为第二次开始充水时刻。
6、当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)60612510.9542 (2)67771520.8392514872 0.9211 59
7、3716 12.0328657670 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
8、19.9594 528236其中(1)表示第一段开始，(2) 表示第二段开始，(3) 表示第三段开始3.计算各时刻点的水流量(加仑/小时) 水流量公式为： 以上25个时刻处的水流量采用差分的方法得到，共分三段分别处理。差分公式为：(1) 对每段前两点采用向前差分公式(2) 对每段最后两点采用向后差分公式(3) 对每段中间点采用中心差分公式得到各点水流量表表3 不同时刻水流量表时间(小时)水流量(加仑/小时)时间(小时)水流量(加仑/小时)时间(小时)水流量(加仑/小时)0 (1)1440410.9542 (2)1946920.839214648 0.9211 11182 12.0328 201
9、95 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
10、) 表示第三段开始4用三次样条拟合流量数据 对表3中25个时刻点的流量数据采用三次样条插值得到一条光滑曲线，作为任意时刻的流量曲线,见图1。 图1 水塔流量图其中*表示数据点，实线为样条曲线5.一天总用水量计算 一天流水总量计算:方法1:直接积分法: 方法2:分段计算法 第一次充水前用水(加仑) 第一次充水后第二次充水前用水(加仑) 22.9581，23.88期间用水(加仑)第一次充水期间用水: 第二次充水期间用水: 23.88，24期间用水: 总共用水两种方法结果相差6.水泵水流量计算 第一次充水期间水塔体积增加 充水时间:第一次充水期间水泵平均流量第二次充水期间水塔体积增加 充水时间:第二
11、次充水期间水泵平均流量则整个充水期间水泵平均流量附Matlab程序%AMCM91Ac=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; 539
12、36,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
13、(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 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-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
14、+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(流量(加仑/小时)
15、;title(水塔流量图);hold offdt=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.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); %计算三次
16、样条插值 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; %第二次充水的平
17、均水流量 p=(p1+p2)/2; %两次充水平均水流量 fprintf(两次充水平均水流量p= %8.2fn,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); %获得
18、离散的时间点,用于积分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.2fn,ss); err=abs(s-ss)/s); fprintf(两种计算法总量相对误差%6.2f%n,err*100);