数学建模算法与应用：预测算法（6）预测习题练习_某居民区的自来水是由一个园柱形的水塔提供。水塔高12.2米,直径17.4米。水塔由水

目录

 一，水塔总水量以及流速预测问题

        1.1、题目

        1.2、建立模型

        1.3、用MATLAB计算，将“-”替换为-1。

        1.4、拟合法

         二、预测产值问题

        2.1、题目

        2.2、建立模型

---

 一，水塔总水量以及流速预测问题

        1.1、题目

        某地区用水管理机构需要对居民的用水速度（单位时间的用水量）和日总用水量进行估计。现有一居民区，其自来水是由一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m。水塔是由水泵根据水塔中的水位自动加水。按照设计，当水塔中的水位降至最低水位，约8.2m时，水泵自动启动加水；当水泵升高到最高水位，约10.8m时，水泵停止工作。

        表1给出的是28个时刻的数据，但由于水泵正向水塔供水，有四个时刻无法测到水位（-）。

        试建立数学模型，来估计居民的用水速度和日总用水量。

        1.2、建立模型

        解：（1）插值法。要估计在任意时刻（包括水泵灌水期间）t 居民的用水速度和日总用水量，分如下三步。

        ①水塔中水的体积的计算。计算水的流量，首先需要计算出水塔中水的体积，即：

             ，式中：D为水塔的直径；h为水塔中水位高度。

        ②水塔中水流速度的估计。居民的用水速度就是水塔中的水流速度，水流速度应该是水塔中水的体积对时间的导数，但由于没有每一时刻水体积的具体数学表达式，只能用差商近似导数。

        由于在两个时段，水泵向水塔供水，无法确定水位的高度，因此在计算水塔中水流速度时要分三段计算。第一时段为0~8.97h，第二时段为10.95~20.84h，第三时段为23.88~25.91h。

        上面计算仅给出流速的离散值，如果需要得到流速的连续型曲线，需要做插值处理，这里可以使用三次样条插值。

        ③日总用水量的计算。日用水量是对水流速度做积分，其积分区间是[0,24]，可以采用数值积分的方法计算。

        1.3、用MATLAB计算，将“-”替换为-1。

计算的MATLAB程序如下： 

clc, clear, close all
a=load('bitter_tea_seeds.txt');
t0=a([1:2:end],:); t0=t0'; t0=t0(:); %提出时间数据，并展开成列向量
h0=a([2:2:end],:); h0=h0'; h0=h0(:); %提出高度数据，并展开成列向量
D=17.4;
V=pi/4*D^2*h0; %计算各时刻的体积
dv=gradient(V,t0); %计算各时刻的数值导数（导数近似值）
no1=find(h0==-1) %找出原始无效数据的地址
no2=[no1(1)-1:no1(2)+1,no1(3)-1:no1(4)+1] %找出导数数据的无效地址
t=t0; t(no2)=[]; %删除导数数据无效地址对应的时间
dv2=-dv; dv2(no2)=[]; %给出各时刻的流速
plot(t,dv2,'*'),title('流速图') %画出流速的散点图
pp=csape(t,dv2); %对流速进行插值
tt=0:0.1:t(end); %给出插值点
fdv=fnval(pp,tt); %计算各插值点的流速值
hold on, plot(tt,fdv) %画出插值曲线
I=trapz(tt(1:241),fdv(1:241)) %计算24小时内总流量的数值积分

结果展示 

 
no1 =
 
      11       
      12       
      24       
      25       
 
 
no2 =
 
  1 至 5 列
 
      10             11             12             13             23       
 
  6 至 8 列
 
      24             25             26       
 
 
I =
 
   21221/17    

        1.4、拟合法

        解：（2）拟合法。要估计在任意时刻（包括水泵灌水期间）t 居民的用水速度和日总用水量，分如下三步。

        ①水塔中水的体积的计算。计算水的流量，首先需要计算出水塔中水的体积，即：

             ，式中：D为水塔的直径；h为水塔中水位高度。

        ②水塔中水流速度的估计。居民的用水速度就是水塔中的水流速度，水流速度应该是水塔中水的体积对时间的导数，但由于没有每一时刻水体积的具体数学表达式，只能用差商近似导数。

        由于在两个时段，水泵向水塔供水，无法确定水位的高度，因此在计算水塔中水流速度时要分三段计算。第一时段为0~8.97h，第二时段为10.95~20.84h，第三时段为23.88~25.91h。  

 

        上面的计算仅给出流速的离散值，流速的散点图如图所示中的 “ * ” 点。如果需要得到流速的连续型曲线，可以拟合多项式曲线，原始数据总共有28个观测值，其中4个无效数据。

        图中总共有20个数据点，这里我们分三段进行三次多项式拟合，应用前6个数据点拟合三次多项式，即在时间区间[0，4.98]上拟合三次多项式；应用第6个数据点到第10个数据点，即在时间区间[4.98，12.03]，拟合第二个三次多项式；应用第10个数据点到第20个数据点，共11个数据点，即在时间区间[12.03，25.91]，拟合第三个三次多项式。

        拟合得到的分段三次多项式曲线如上图所示。

        ③日总用水量的计算。日用水量是对水流速度做积分，其积分区间是[0,24]，可以采用数值积分的方法计算。

MATLAB程序如下：

clc, clear, close all
a=load('bitter_tea_seeds.txt');
t0=a([1:2:end],:); t0=t0'; t0=t0(:); %提出时间数据，并展开成列向量
h0=a([2:2:end],:); h0=h0'; h0=h0(:); %提出高度数据，并展开成列向量
D=17.4;
V=pi/4*D^2*h0; %计算各时刻的体积
dv=gradient(V,t0); %计算各时刻的数值导数（导数近似值）
no1=find(h0==-1) %找出原始无效数据的地址
no2=[no1(1)-1:no1(2)+1,no1(3)-1:no1(4)+1] %找出导数数据的无效地址
t=t0; t(no2)=[]; %删除导数数据无效地址对应的时间
dv2=-dv; dv2(no2)=[]; %给出各时刻的流速
hold on, plot(t,dv2,'*') %画出流速的散点图
a1=polyfit(t(1:6),dv2(1:6),3); %拟合第一个多项式的系数
a2=polyfit(t(6:10),dv2(6:10),3); %拟合第二个多项式的系数
a3=polyfit(t(10:20),dv2(10:20),3); %拟合第三个多项式的系数
dvf1=polyval(a1,[t(1):0.1:t(6)]); %计算第一个多项式的函数值
dvf2=polyval(a2,[t(6):0.1:t(10)]); %计算第二个多项式的函数值
dvf3=polyval(a3,[t(10):0.1:t(end)]); %计算第三个多项式的函数值
tt=t(1):0.1:t(end); dvf=[dvf1,dvf2,dvf3];
plot(tt,dvf) %画出拟合的三个分段多项式曲线
I=trapz(tt(1:241),dvf(1:241)) %计算24小时内总流量的数值积分

二、预测产值问题

        2.1、题目

         某大型企业1997---2000年的产值资料如表所示，试建立GM（1，1）预测模型，预测该企业2001---2005年的产值。

        2.2、建立模型（笔记打印，有曲面）

        解：1.级比检验

        建立原始序列数据

        2.GM（1，1）建模

         3.模型检验

        模型的各种检验指标值的计算结果见表2。经检验，该模型的精度较高，可进行预测和预报。

          4.预测值

        2001---2005年预测值见表2。

MATLAB程序设计代码，如下：

clc,clear, format long g
x0=[27260	29547	32411	35388]'; %注意这里为列向量
n=length(x0);
lamda=x0(1:n-1)./x0(2:n)  %计算级比
range=minmax(lamda')  %计算级比的范围
theta=[exp(-2/(n+1)),exp(2/(n+1))]  % 计算级比的容许区间
x1=cumsum(x0)  %累加运算
B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];
Y=x0(2:n);
u=B\Y  %拟合参数u(1)=a,u(2)=b
syms x(t)
x=dsolve(diff(x)+u(1)*x==u(2),x(0)==x0(1)); %求微分方程的符号解
xt=vpa(x,6) %以小数格式显示微分方程的解
yuce1=subs(x,t,[0:n+4]); %求已知数据和未来5期的预测值
yuce1=double(yuce1); %符号数转换成数值类型，否则无法作差分运算
yuce=[x0(1),diff(yuce1)]  %差分运算，还原数据
epsilon=x0'-yuce(1:n)    %计算已知数据预测的残差
delta=abs(epsilon./x0')  %计算相对误差
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'  %计算级比偏差值，u(1)=a
yhat=yuce(n+1:end)  %提取未来5期的预测值

 结果展示：

 
lamda =
 
         0.922597894879345
         0.911634938755361
         0.915875438001582
 
 
range =
 
         0.911634938755361         0.922597894879345
 
 
theta =
 
         0.670320046035639          1.49182469764127
 
 
x1 =
 
       27260
       56807
       89218
      124606
 
 
u =
 
       -0.0899951723868726
          25790.2838424515
 
 
xt =
 
313834.0*exp(0.0899952*t) - 286574.0
 
 
yuce =
 
  1 至 3 列
 
                     27260          29553.4420565748          32336.4601849797
 
  4 至 6 列
 
          35381.5523516033          38713.3978069354          42358.9998218412
 
  7 至 9 列
 
          46347.9045382397          50712.4404287318          55487.9803059021
 
 
epsilon =
 
  1 至 3 列
 
                         0         -6.44205657481507          74.5398150203109
 
  4 列
 
          6.44764839668642
 
 
delta =
 
  1 至 3 列
 
                         0      0.000218027433404917       0.00229983076795875
 
  4 列
 
      0.000182198722637233
 
 
rho =
 
      -0.00953940978346313       0.00245664647930999      -0.00218345852198909
 
 
yhat =
 
  1 至 3 列
 
          38713.3978069354          42358.9998218412          46347.9045382397
 
  4 至 5 列
 
          50712.4404287318          55487.9803059021