城市供水量预测

答卷编号：

论文题目：城市供水量预测

指导教师：teacher

参赛学校：南昌大学

报名序号：1657

证书邮寄地址（学校组织请填写组队负责人）：

邮政编码330031 南昌大学前湖180405室

卢雄 （收）

答卷编号：

城市供水量预测

摘要

随着科学技术的发展，人们生活水平的提高，水资源需求量越来越大，导致如今各大城市水资源严重匮乏，为了节约水资源，某城市对7年内的总供水量以及两个供水厂各自的供水量做了调查，本文利用这7年内供水量数据，建立了价格与供水量之间的函数图像，为了使模型简单化我们只分离出每一年同一个月的供水量，用matlab软件对函数图像进行多项式拟合，并用拟合程度最好与变化率方差最大找出函数关系式。发现理论图像与实际图像拟合得很好，能够很好的说明供水量受价格影响的变化趋势，并能准确预测出未来一年该月受价格影响的情况，总供水量以及两厂各自供水量的具体变化值。根据用上述方法求出的函数，能做到对水价进行合理的调整来控制未来的供水量。

一 问题提出

根据题目给出的数据附表，建立适当的数学模型来预测2007年1月份城市的计划供水量，同理可以推出2007年1月份城市中每个水厂的计划供水量，并由该城市8月份的供水量曲线得出在限定的供水量下所对应的水价，从而确定合理的水价调整方案。

二 问题相关的概念及原理

在(美) Frank R. Giordano, Maurice D. Weir, William P. Fox著的数学建模一书中的概念

1 误差平方和SSE：SSE=(yi-f(xi))2

i=1nn

2 修正平方和SST：SST=(yi-y)2

i=1

3 拟合程度：R=1-

SSE

SST

其中：yi为某一个水价对应得实际供水量， f(xi)为某水价对应的理论供水量，

y为实际供水量的平均值。

由于题目需要我们自己给出一下定义： 4 变化率方差：S=(f(xi)-f(x))2

i=1n

其中f(xi)为某水价对应导函数的理论供水量，f(xi)为某水价对应导函数的理论供水量的平均值。

变化率方差S是是用来说明模拟出来的函数在定义域变化的缓度，波动太大，

和实际的情况也不符，因此我们给出这个定义以便找出最符合实际的函数。后面的结果也验证我们定义的正确性。

5 拟合程度和变化率方差之比：T=R/S。

由于在实际操作中，由于符合上述R，和S的函数有多个，而对于拟合程度R越大越好（0

三 模型假设

根据实际情况，我们做出了如下假设：

忽略2000年到2007年在同一个月的温度变化;

某月水厂的供水量就是居民在该月所用的水的总量。

四 关键词

多项式拟合 水价 误差平方和 修正平方和 拟合程度 变化率方差 拟合程度/变化率方差之比

五 问题分析

模型一：

首先，预测2007年一月份城市的计划供水量。利用分离法，分离出从2000年到2006年一月份的供水量，利用多项式拟合找出函数关系式，由函数关系式计算出2007年一月份的计划供水量。

其次，预测2007年1月份城市中每个水厂的计划供水量。利用分离法，分离出从2000年到2006年一月份各厂的供水量，利用上述方法预测出一月份各厂的计划供水量。

模型二：

怎样合理调整水价使2007年8月份的供水量不超过5045万吨。利用分离法，分离出从2000年到2006年八月份的供水量，利用多项式拟合找出函数关系式，再把所要求限制条件代入函数关系式，求出合理的价格。

模型求解

现在来进行模型求解，首先说明一下用到的matlab 中的函数

使用polyfit()函数用来求解多项式拟合问题，使用polyval()函数用来计算在点x上的对应的函数值。

相关参数：

Y：每年一月的城市供水量几何平均数 T：拟合程度与变化率方差之比 X：水价

（另外一些符号已经在“问题相关的概念及原理”中给出，就不在这里重复。） 下面是求解过程：

用matlab求得每年一月的几何平均数

Y= 1.3001 1.3498 1.3854 1.4110 1.4302 1.4527 1.4572 代入拟合曲线程序。（Y的单位是107吨） 得:T1=479.4976 T2=454.5882

T1>T2所以曲线一的拟合度更好，即：曲线一是最佳拟合曲线。于是: Y=-0.0133X2+0.1598X+0.9919

根据上述函数关系式，代入2007年1月份的所对应的价格X=5.2得出

Y=4.53593×10，所以我们预测2007年1月份的计划供水量为4.53593×10吨。

下面是在T=479.4976，R=0.99646,S=0.0020747情况下的最佳图

7

7

同理可得，1月份一号水厂最佳供水曲线为：

2

Y=-0.0821X+0.9799X+6.1300

根据上述函数关系式我们不难预测出2007年1月份一号水厂的计划供水量为2.791696×10吨。

下面是一号水厂最佳拟合曲线图：

7

同上可求出2000年71月份二号水厂最佳供水曲线为： Y=-0.0507X2+0.6150X+3.7458

从而预测出2007年1月份二号水厂的计划供水量为1.727520×107吨。

下面是其对应的最佳拟合曲线图：

先用上面给出的方法求出2007年8月份城市的最佳拟合程度供水量曲线为 Y=0.0148X2-0.0395X+1.4016

下面是其对应的图：

把Y=50.45/31=1.6273代入供水量曲线, 由matlab求得水价X=5.4614元。

考虑到实际生活中的情况，我们可以把价格调整为5.5元。

六 结论分析、评价及推广

结论分析：

题目给出的数据太多太杂，为了把问题简单化，我们先是忽略温度因素，之后用分离法，分离出每年的同一个月份，把时间函数变成一个只与价格有关的函数。这样问题就简单的变成了一个一元函数问题，对于一元函数的处理，我们采用了《数值计算方法》一书中的线性最小二乘拟合法来拟合找出最佳曲线，来预测将来的用水量。先由多项式模拟曲线入手，分析了供水量关于价格的函数，找出最佳拟合曲线，然后由求得的最佳拟合曲线出发，可求得某个价格所对应供水量以及限制一定量供水量所对应的价格。 评价：

我们采用的这种分离法是值得鼓励的，它把一个复杂问题简单化了好多，做起来，也轻松多了，在运算过程中，我们采用数型结合求解，使整个论文看上去，直观，易懂，这也是值得推崇的。但是在处理问题，把问题想得太过简单，忽略了太多的因素，如温度等。而处理所用的点都过少，使计算精确度也下降了不少。 推广：

对于给定的一组数据{(Xi,Yi)}im假设拟合函1(X1,X2,„,Xm∈[a,b]互不相同)，数的形式为:

(x)a00(x)a11(x)ann(x)

其中{k(x)}nk0(n

E2(a0,a1,,an)[(xi)yi][ajj(xi)yi]2

i1

i1

j0

m

2

mn

为极小的问题成为线性最小二乘问题。

七 心得体会

看到题目后，经过讨论，我们决定选接近我们生活的题目——供水量预测作为我们建模的内容，确定题目后，我们三人开始讨论，去图书馆找相关资料，通过一天的努力，终于建立了初步模型，虽然开始感到很盲目，没有一个明确的目标，甚至无从下手，可是通过互相鼓励，我们找回了属于我们的自信，这让我体会到团队精神在数学建模中的重要性，同时，除了相互鼓励，在很多时候，一个人的思考是不全面的，只有大家一起讨论才有可能把问题搞清楚，因此无论做任何板块，三个人都要一起齐心才行。第二天，为了防止盲目的思想重新出现，我们重新安排了时间，让我们有一种事半功倍的感觉。的确，做任何事情，合理的时间安排非常重要，建模也是一样，事先要做好一个规划，确定好哪一段时间要完成哪些内容，这样做才能游刃有余，终于在第二天，我们建立了完整的模型，进行了必要的数据处理，为第三天的论文做好了铺垫，在第三天，我们终于做好了所有的工作，在我们每个人看来这都是 比较难熬的三天，但没有人想过放弃，其实，数学建模锻炼和考察的是我们每个人的综合素质，从这次数学建模中我们不只学到的是数学知识，更重要的是建模思想，甚至是心理素质。

附：模型中计算拟合曲线以及预测点的函数值，以模型一中一月份城市的供水量的程序为例，其它模型中的只需在此基础上修改些数据即可 Matlab源程序：

1.求一月份城市供水量几何平均值 function [ y ] = GeometricMean(x)

x0=[1221790 1282410 1241980 1265880 1301360 1298670 1273770 1300620 1301750 1318300 1327550 1356910 1329360 1312580 1330460 1316710 1293410 1303150 1304690 1301800 1286760 1303180 1302230 1285860 1297670 1302460 1304330 1297080 1315840 1307490 1325940];

x1=[1224777.6 1270069.2 1280584.8 1310641.2 1282716.9 1319987.3 1342115.6 1353800.7 1359700.2 1378040.4 1372140.9 1353780 1367618 1360290.2 1378402.7 1377284.9 1401762.6 1376901.9 1370505.6 1387169.1 1362391.2 1366862.4 1387624.5 1361056.1 1351130.4 1347156 1376694.9 1369760.4 1356605.6 1357992.5 1354690.8 ];

x2=[1302579.2 1324042.2 1337937.9 1381181.4 1379235.8 1429176.3 1384325.1

1384888.3 1379553.3 1358940.2 1346775 1366282.2 1377157.1 1370828.8 1396572.2 1396152.3 1384069.1 1401753.6 1433538.6 1446799.4 1439303.7 1442232.3 1449257 1440921.6 1437132.8 1453230.1 1435105.3 1313812.5 1297152 1338306.6 1341624.3];

x3=[1331466.85 1392996.15 1351959.7 1376218.2 1412230.4 1409500.05

1384226.55 1411479.3 1412626.25 1429424.5 1438813.25 1468613.65 1440650.4 1423618.7 1441766.9 1427810.65 1404161.15 1414047.25 1415610.35 1412677 1397411.4 1414077.7 1413113.45 1396497.9 1408485.05 1413346.9 1415244.95 1407886.2 1426927.6 1418452.35 1437179.1];

x4=[1306987.2 1351622.4 1361985.6 1391606.4 1364086.8 1400817 1422624.6 1434140.4 1439954.4 1458028.8 1452214.8 1434120 1447757.4 1440535.8 1458385.8 1457284.2 1481407.2 1456906.8 1450603.2 1467025.2 1442606.4 1447012.8 1467474 1441290.6 1431508.8 1427592 1456702.8 1449868.8 1436904.6 1438271.4 1435017.6];

x5=[1391050.08 1404728.64 1447296.48 1445381.28 1494541.44 1450391.04 1450945.44 1445693.76 1425402.72 1413427.68 1432630.08 1443335.04 1437105.6 1462446.72 1462033.44 1450139.04 1467547.2 1498835.52 1511889.12 1504510.56 1507393.44 1514308.32 1506103.2 1502373.6 1518219.36 1500377.76 1380980.16 1364580 1405091.52 1408357.44 1407157.92];

x6=[1436961.729 1397505.082 1420829.689 1455455.508 1452830.27

1428529.737 1454733.324 1455836.119 1471987.678 1481014.984 1509668.142 1482781.408 1466405.387 1483854.926 1470435.957 1447696.904 1457202.414 1458705.338 1455884.915 1441207.003 1457231.692 1456304.563 1440328.671 1451854.345 1456529.026 1458354.005 1451278.549 1469586.902 1461437.928 1479443.745 1465117.166];

y=[1 1 1 1 1 1 1]

for k=1:length(x0) y(1)=y(1)*x0(k) end

for k=1:length(x1) y(2)=y(2)*x1(k) end

for k=1:length(x2) y(3)=y(3)*x2(k) end

for k=1:length(x3) y(4)=y(4)*x3(k) end

for k=1:length(x4) y(5)=y(5)*x4(k) end

for k=1:length(x5) y(6)=y(6)*x5(k) end

for k=1:length(x6) y(7)=y(7)*x6(k) end

y(1)=exp(1.0/length(x0)*log(y(1))) y(2)=exp(1.0/length(x1)*log(y(2))) y(3)=exp(1.0/length(x2)*log(y(3))) y(4)=exp(1.0/length(x3)*log(y(4))) y(5)=exp(1.0/length(x4)*log(y(5))) y(6)=exp(1.0/length(x5)*log(y(6))) y(7)=exp(1.0/length(x6)*log(y(7)))

2.求拟合曲线，并求出预测点的函数值 function [ y,x] =DateMultiRegress x=[3.0 3.4 3.9 4.3 4.7 5.0]

y=[1.3803 1.4353 1.4755 1.5064 1.5392 1.5516];

a=2:0.5:6

% 下面拟合二次和三次多项式 for n=2:3

f(:,n)=polyval(polyfit(x,y,n),a)'; z=polyfit(x,y,n); disp('z='); disp(z);

plot(a,f(:,n),x,y,'o')

% 下面求误差平方和sse,修正平方和sst sse=0 sst=0

g(:,n)=polyval(z,x) ym=mean(y)

for j=1:length(y)

psse=y(j)-g(j,n);

sse=sse+psse^2

psst=y(j)-ym;

sst=sst+psst^2

end

r=1-sse/sst %求拟合度r

%下面求变化率方差subsst

subsst=0

h(:,n)=polyval(polyder(z),x)

eh=h(:,n)

for j=1:length(x)

subpsst=h(j,n)-mean(eh);

subsst=subsst+subpsst^2

end

title(['Poly. regression, deg=',int2str(n),' R=',num2str(r),' S=',num2str(subsst),' R/S=',num2str(r/subsst)])

xlabel('Price'), ylabel('Total'), grid

z %显示所求函数的各项系数

o=polyval(z,5.2) %求出预测点的函数值

pause

end

3. 求拟合曲线，并根据拟合曲线求出在限定的供水量下所对应的水价 function [ y,x] =DateMultiRegress

x=[3.0 3.4 3.9 4.3 4.7 5.0]

y=[1.3803 1.4353 1.4755 1.5064 1.5392 1.5516];

a=2:0.5:6

% 下面拟合二次和三次多项式

for n=2:3

f(:,n)=polyval(polyfit(x,y,n),a)';

z=polyfit(x,y,n);

disp('z=');

disp(z);

plot(a,f(:,n),x,y,'o')

% 下面求误差平方和sse,修正平方和sst

sse=0

sst=0

g(:,n)=polyval(z,x)

ym=mean(y)

for j=1:length(y)

psse=y(j)-g(j,n);

sse=sse+psse^2

psst=y(j)-ym;

sst=sst+psst^2

end

r=1-sse/sst %求拟合度r

%下面求变化率方差subsst

subsst=0

h(:,n)=polyval(polyder(z),x)

eh=h(:,n)

for j=1:length(x)

subpsst=h(j,n)-mean(eh);

subsst=subsst+subpsst^2

end

title(['Poly. regression, deg=',int2str(n),' R=',num2str(r),' S=',num2str(subsst),' R/S=',num2str(r/subsst)])

xlabel('Price'), ylabel('Total'), grid

z %显示所求函数的各项系数

roots() %求出预测点的函数值

pause

end

参考文献

[1] 《数学建模》(美) Frank R. Giordano, Maurice D. Weir, William P. Fox著；

[2] 《数值计算方法》（下册） 林成森 编著

答卷编号：

论文题目：城市供水量预测

指导教师：teacher

参赛学校：南昌大学

报名序号：1657

证书邮寄地址（学校组织请填写组队负责人）：

邮政编码330031 南昌大学前湖180405室

卢雄 （收）

答卷编号：

城市供水量预测

摘要

随着科学技术的发展，人们生活水平的提高，水资源需求量越来越大，导致如今各大城市水资源严重匮乏，为了节约水资源，某城市对7年内的总供水量以及两个供水厂各自的供水量做了调查，本文利用这7年内供水量数据，建立了价格与供水量之间的函数图像，为了使模型简单化我们只分离出每一年同一个月的供水量，用matlab软件对函数图像进行多项式拟合，并用拟合程度最好与变化率方差最大找出函数关系式。发现理论图像与实际图像拟合得很好，能够很好的说明供水量受价格影响的变化趋势，并能准确预测出未来一年该月受价格影响的情况，总供水量以及两厂各自供水量的具体变化值。根据用上述方法求出的函数，能做到对水价进行合理的调整来控制未来的供水量。

一 问题提出

根据题目给出的数据附表，建立适当的数学模型来预测2007年1月份城市的计划供水量，同理可以推出2007年1月份城市中每个水厂的计划供水量，并由该城市8月份的供水量曲线得出在限定的供水量下所对应的水价，从而确定合理的水价调整方案。

二 问题相关的概念及原理

在(美) Frank R. Giordano, Maurice D. Weir, William P. Fox著的数学建模一书中的概念

1 误差平方和SSE：SSE=(yi-f(xi))2

i=1nn

2 修正平方和SST：SST=(yi-y)2

i=1

3 拟合程度：R=1-

SSE

SST

其中：yi为某一个水价对应得实际供水量， f(xi)为某水价对应的理论供水量，

y为实际供水量的平均值。

由于题目需要我们自己给出一下定义： 4 变化率方差：S=(f(xi)-f(x))2

i=1n

其中f(xi)为某水价对应导函数的理论供水量，f(xi)为某水价对应导函数的理论供水量的平均值。

变化率方差S是是用来说明模拟出来的函数在定义域变化的缓度，波动太大，

和实际的情况也不符，因此我们给出这个定义以便找出最符合实际的函数。后面的结果也验证我们定义的正确性。

5 拟合程度和变化率方差之比：T=R/S。

由于在实际操作中，由于符合上述R，和S的函数有多个，而对于拟合程度R越大越好（0

三 模型假设

根据实际情况，我们做出了如下假设：

忽略2000年到2007年在同一个月的温度变化;

某月水厂的供水量就是居民在该月所用的水的总量。

四 关键词

多项式拟合 水价 误差平方和 修正平方和 拟合程度 变化率方差 拟合程度/变化率方差之比

五 问题分析

模型一：

首先，预测2007年一月份城市的计划供水量。利用分离法，分离出从2000年到2006年一月份的供水量，利用多项式拟合找出函数关系式，由函数关系式计算出2007年一月份的计划供水量。

其次，预测2007年1月份城市中每个水厂的计划供水量。利用分离法，分离出从2000年到2006年一月份各厂的供水量，利用上述方法预测出一月份各厂的计划供水量。

模型二：

怎样合理调整水价使2007年8月份的供水量不超过5045万吨。利用分离法，分离出从2000年到2006年八月份的供水量，利用多项式拟合找出函数关系式，再把所要求限制条件代入函数关系式，求出合理的价格。

模型求解

现在来进行模型求解，首先说明一下用到的matlab 中的函数

使用polyfit()函数用来求解多项式拟合问题，使用polyval()函数用来计算在点x上的对应的函数值。

相关参数：

Y：每年一月的城市供水量几何平均数 T：拟合程度与变化率方差之比 X：水价

（另外一些符号已经在“问题相关的概念及原理”中给出，就不在这里重复。） 下面是求解过程：

用matlab求得每年一月的几何平均数

Y= 1.3001 1.3498 1.3854 1.4110 1.4302 1.4527 1.4572 代入拟合曲线程序。（Y的单位是107吨） 得:T1=479.4976 T2=454.5882

T1>T2所以曲线一的拟合度更好，即：曲线一是最佳拟合曲线。于是: Y=-0.0133X2+0.1598X+0.9919

根据上述函数关系式，代入2007年1月份的所对应的价格X=5.2得出

Y=4.53593×10，所以我们预测2007年1月份的计划供水量为4.53593×10吨。

下面是在T=479.4976，R=0.99646,S=0.0020747情况下的最佳图

7

7

同理可得，1月份一号水厂最佳供水曲线为：

2

Y=-0.0821X+0.9799X+6.1300

根据上述函数关系式我们不难预测出2007年1月份一号水厂的计划供水量为2.791696×10吨。

下面是一号水厂最佳拟合曲线图：

7

同上可求出2000年71月份二号水厂最佳供水曲线为： Y=-0.0507X2+0.6150X+3.7458

从而预测出2007年1月份二号水厂的计划供水量为1.727520×107吨。

下面是其对应的最佳拟合曲线图：

先用上面给出的方法求出2007年8月份城市的最佳拟合程度供水量曲线为 Y=0.0148X2-0.0395X+1.4016

下面是其对应的图：

把Y=50.45/31=1.6273代入供水量曲线, 由matlab求得水价X=5.4614元。

考虑到实际生活中的情况，我们可以把价格调整为5.5元。

六 结论分析、评价及推广

结论分析：

题目给出的数据太多太杂，为了把问题简单化，我们先是忽略温度因素，之后用分离法，分离出每年的同一个月份，把时间函数变成一个只与价格有关的函数。这样问题就简单的变成了一个一元函数问题，对于一元函数的处理，我们采用了《数值计算方法》一书中的线性最小二乘拟合法来拟合找出最佳曲线，来预测将来的用水量。先由多项式模拟曲线入手，分析了供水量关于价格的函数，找出最佳拟合曲线，然后由求得的最佳拟合曲线出发，可求得某个价格所对应供水量以及限制一定量供水量所对应的价格。 评价：

我们采用的这种分离法是值得鼓励的，它把一个复杂问题简单化了好多，做起来，也轻松多了，在运算过程中，我们采用数型结合求解，使整个论文看上去，直观，易懂，这也是值得推崇的。但是在处理问题，把问题想得太过简单，忽略了太多的因素，如温度等。而处理所用的点都过少，使计算精确度也下降了不少。 推广：

对于给定的一组数据{(Xi,Yi)}im假设拟合函1(X1,X2,„,Xm∈[a,b]互不相同)，数的形式为:

(x)a00(x)a11(x)ann(x)

其中{k(x)}nk0(n

E2(a0,a1,,an)[(xi)yi][ajj(xi)yi]2

i1

i1

j0

m

2

mn

为极小的问题成为线性最小二乘问题。

七 心得体会

看到题目后，经过讨论，我们决定选接近我们生活的题目——供水量预测作为我们建模的内容，确定题目后，我们三人开始讨论，去图书馆找相关资料，通过一天的努力，终于建立了初步模型，虽然开始感到很盲目，没有一个明确的目标，甚至无从下手，可是通过互相鼓励，我们找回了属于我们的自信，这让我体会到团队精神在数学建模中的重要性，同时，除了相互鼓励，在很多时候，一个人的思考是不全面的，只有大家一起讨论才有可能把问题搞清楚，因此无论做任何板块，三个人都要一起齐心才行。第二天，为了防止盲目的思想重新出现，我们重新安排了时间，让我们有一种事半功倍的感觉。的确，做任何事情，合理的时间安排非常重要，建模也是一样，事先要做好一个规划，确定好哪一段时间要完成哪些内容，这样做才能游刃有余，终于在第二天，我们建立了完整的模型，进行了必要的数据处理，为第三天的论文做好了铺垫，在第三天，我们终于做好了所有的工作，在我们每个人看来这都是 比较难熬的三天，但没有人想过放弃，其实，数学建模锻炼和考察的是我们每个人的综合素质，从这次数学建模中我们不只学到的是数学知识，更重要的是建模思想，甚至是心理素质。

附：模型中计算拟合曲线以及预测点的函数值，以模型一中一月份城市的供水量的程序为例，其它模型中的只需在此基础上修改些数据即可 Matlab源程序：

1.求一月份城市供水量几何平均值 function [ y ] = GeometricMean(x)

x0=[1221790 1282410 1241980 1265880 1301360 1298670 1273770 1300620 1301750 1318300 1327550 1356910 1329360 1312580 1330460 1316710 1293410 1303150 1304690 1301800 1286760 1303180 1302230 1285860 1297670 1302460 1304330 1297080 1315840 1307490 1325940];

x1=[1224777.6 1270069.2 1280584.8 1310641.2 1282716.9 1319987.3 1342115.6 1353800.7 1359700.2 1378040.4 1372140.9 1353780 1367618 1360290.2 1378402.7 1377284.9 1401762.6 1376901.9 1370505.6 1387169.1 1362391.2 1366862.4 1387624.5 1361056.1 1351130.4 1347156 1376694.9 1369760.4 1356605.6 1357992.5 1354690.8 ];

x2=[1302579.2 1324042.2 1337937.9 1381181.4 1379235.8 1429176.3 1384325.1

1384888.3 1379553.3 1358940.2 1346775 1366282.2 1377157.1 1370828.8 1396572.2 1396152.3 1384069.1 1401753.6 1433538.6 1446799.4 1439303.7 1442232.3 1449257 1440921.6 1437132.8 1453230.1 1435105.3 1313812.5 1297152 1338306.6 1341624.3];

x3=[1331466.85 1392996.15 1351959.7 1376218.2 1412230.4 1409500.05

1384226.55 1411479.3 1412626.25 1429424.5 1438813.25 1468613.65 1440650.4 1423618.7 1441766.9 1427810.65 1404161.15 1414047.25 1415610.35 1412677 1397411.4 1414077.7 1413113.45 1396497.9 1408485.05 1413346.9 1415244.95 1407886.2 1426927.6 1418452.35 1437179.1];

x4=[1306987.2 1351622.4 1361985.6 1391606.4 1364086.8 1400817 1422624.6 1434140.4 1439954.4 1458028.8 1452214.8 1434120 1447757.4 1440535.8 1458385.8 1457284.2 1481407.2 1456906.8 1450603.2 1467025.2 1442606.4 1447012.8 1467474 1441290.6 1431508.8 1427592 1456702.8 1449868.8 1436904.6 1438271.4 1435017.6];

x5=[1391050.08 1404728.64 1447296.48 1445381.28 1494541.44 1450391.04 1450945.44 1445693.76 1425402.72 1413427.68 1432630.08 1443335.04 1437105.6 1462446.72 1462033.44 1450139.04 1467547.2 1498835.52 1511889.12 1504510.56 1507393.44 1514308.32 1506103.2 1502373.6 1518219.36 1500377.76 1380980.16 1364580 1405091.52 1408357.44 1407157.92];

x6=[1436961.729 1397505.082 1420829.689 1455455.508 1452830.27

1428529.737 1454733.324 1455836.119 1471987.678 1481014.984 1509668.142 1482781.408 1466405.387 1483854.926 1470435.957 1447696.904 1457202.414 1458705.338 1455884.915 1441207.003 1457231.692 1456304.563 1440328.671 1451854.345 1456529.026 1458354.005 1451278.549 1469586.902 1461437.928 1479443.745 1465117.166];

y=[1 1 1 1 1 1 1]

for k=1:length(x0) y(1)=y(1)*x0(k) end

for k=1:length(x1) y(2)=y(2)*x1(k) end

for k=1:length(x2) y(3)=y(3)*x2(k) end

for k=1:length(x3) y(4)=y(4)*x3(k) end

for k=1:length(x4) y(5)=y(5)*x4(k) end

for k=1:length(x5) y(6)=y(6)*x5(k) end

for k=1:length(x6) y(7)=y(7)*x6(k) end

y(1)=exp(1.0/length(x0)*log(y(1))) y(2)=exp(1.0/length(x1)*log(y(2))) y(3)=exp(1.0/length(x2)*log(y(3))) y(4)=exp(1.0/length(x3)*log(y(4))) y(5)=exp(1.0/length(x4)*log(y(5))) y(6)=exp(1.0/length(x5)*log(y(6))) y(7)=exp(1.0/length(x6)*log(y(7)))

2.求拟合曲线，并求出预测点的函数值 function [ y,x] =DateMultiRegress x=[3.0 3.4 3.9 4.3 4.7 5.0]

y=[1.3803 1.4353 1.4755 1.5064 1.5392 1.5516];

a=2:0.5:6

% 下面拟合二次和三次多项式 for n=2:3

f(:,n)=polyval(polyfit(x,y,n),a)'; z=polyfit(x,y,n); disp('z='); disp(z);

plot(a,f(:,n),x,y,'o')

% 下面求误差平方和sse,修正平方和sst sse=0 sst=0

g(:,n)=polyval(z,x) ym=mean(y)

for j=1:length(y)

psse=y(j)-g(j,n);

sse=sse+psse^2

psst=y(j)-ym;

sst=sst+psst^2

end

r=1-sse/sst %求拟合度r

%下面求变化率方差subsst

subsst=0

h(:,n)=polyval(polyder(z),x)

eh=h(:,n)

for j=1:length(x)

subpsst=h(j,n)-mean(eh);

subsst=subsst+subpsst^2

end

title(['Poly. regression, deg=',int2str(n),' R=',num2str(r),' S=',num2str(subsst),' R/S=',num2str(r/subsst)])

xlabel('Price'), ylabel('Total'), grid

z %显示所求函数的各项系数

o=polyval(z,5.2) %求出预测点的函数值

pause

end

3. 求拟合曲线，并根据拟合曲线求出在限定的供水量下所对应的水价 function [ y,x] =DateMultiRegress

x=[3.0 3.4 3.9 4.3 4.7 5.0]

y=[1.3803 1.4353 1.4755 1.5064 1.5392 1.5516];

a=2:0.5:6

% 下面拟合二次和三次多项式

for n=2:3

f(:,n)=polyval(polyfit(x,y,n),a)';

z=polyfit(x,y,n);

disp('z=');

disp(z);

plot(a,f(:,n),x,y,'o')

% 下面求误差平方和sse,修正平方和sst

sse=0

sst=0

g(:,n)=polyval(z,x)

ym=mean(y)

for j=1:length(y)

psse=y(j)-g(j,n);

sse=sse+psse^2

psst=y(j)-ym;

sst=sst+psst^2

end

r=1-sse/sst %求拟合度r

%下面求变化率方差subsst

subsst=0

h(:,n)=polyval(polyder(z),x)

eh=h(:,n)

for j=1:length(x)

subpsst=h(j,n)-mean(eh);

subsst=subsst+subpsst^2

end

title(['Poly. regression, deg=',int2str(n),' R=',num2str(r),' S=',num2str(subsst),' R/S=',num2str(r/subsst)])

xlabel('Price'), ylabel('Total'), grid

z %显示所求函数的各项系数

roots() %求出预测点的函数值

pause

end

参考文献

[1] 《数学建模》(美) Frank R. Giordano, Maurice D. Weir, William P. Fox著；

[2] 《数值计算方法》（下册） 林成森 编著