目录
一，水塔总水量以及流速预测问题
1.1、题目
某地区用水管理机构需要对居民的用水速度（单位时间的用水量）和日总用水量进行估计。现有一居民区，其自来水是由一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m。水塔是由水泵根据水塔中的水位自动加水。按照设计，当水塔中的水位降至最低水位，约8.2m时，水泵自动启动加水；当水泵升高到最高水位，约10.8m时，水泵停止工作。
表1给出的是28个时刻的数据，但由于水泵正向水塔供水，有四个时刻无法测到水位（-）。
| 时刻（t）/h | 0 | 0.92 | 1.84 | 2.95 | 3.87 | 4.98 | 5.90 |
|---|---|---|---|---|---|---|---|
| 水位/m | 9.68 | 9.48 | 9.31 | 9.13 | 8.98 | 8.81 | 8.69 |
| 时刻（t）/h | 7.01 | 7.93 | 8.97 | 9.98 | 10.92 | 10.95 | 12.03 |
| 水位/m | 8.52 | 8.39 | 8.22 | - | - | 10.82 | 10.5 |
| 时刻（t）/h | 12.95 | 13.88 | 14.98 | 15.9 | 16.83 | 17.93 | 19.04 |
| 水位/m | 10.21 | 9.94 | 9.65 | 9.41 | 9.18 | 8.92 | 8.66 |
| 时刻（t）/h | 19.96 | 20.84 | 22.01 | 22.96 | 23.88 | 24.99 | 25.91 |
| 水位/m | 8.43 | 8.22 | - | - | 10.59 | 10.35 | 10.18 |
试建立数学模型，来估计居民的用水速度和日总用水量。
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年的产值。
| 年份 | 1997 | 1998 | 1999 | 2000 |
|---|---|---|---|---|
| 产值/万元 | 27260 | 29547 | 32411 | 35388 |
2.2、建立模型（笔记打印，有曲面）
解：1.级比检验
建立原始序列数据
2.GM（1，1）建模
3.模型检验
模型的各种检验指标值的计算结果见表2。经检验，该模型的精度较高，可进行预测和预报。
| 年份 | 原始值 | 预测值 | 残差 | 相对误差 | 级比偏差 |
|---|---|---|---|---|---|
| 1997 | 27260 | 27260 | 0 | 0 | |
| 1998 | 29547 | 29553.4421 | -6.4421 | 0.0002 | -0.0095 |
| 1999 | 32411 | 32336.4602 | 74.5398 | 0.0023 | 0.0025 |
| 2000 | 35388 | 35381.5524 | 6.4476 | 0.0002 | -0.0022 |
4.预测值
2001---2005年预测值见表2。
| 年份 | 2001 | 2002 | 2003 | 2004 | 2005 |
|---|---|---|---|---|---|
| 产值/万元 | 38713.3978 | 42358.9998 | 46347.9045 | 50712.4404 | 55487.9803 |
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
8404
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
