matlab数学实验

《管理数学实验》实验报告
班级姓名
实验1：MATLAB的数值运算
【实验目的】
（1）掌握MATLAB变量的使用
（2）掌握MATLAB数组的创建，
（3）掌握MA TLAB数组和矩阵的运算。

（4）熟悉MATLAB多项式的运用
【实验原理】
矩阵运算和数组运算在MA TLAB中属于两种不同类型的运算，数组的运算是从数组元素出发，针对每个元素进行运算，矩阵的运算是从矩阵的整体出发，依照线性代数的运算规则进行。

【实验步骤】
（1）使用冒号生成法和定数线性采样法生成一维数组。

（2）使用MA TLAB提供的库函数reshape，将一维数组转换为二维和三维数组。

（3）使用逐个元素输入法生成给定变量，并对变量进行指定的算术运算、关系运算、逻辑运算。

（4）使用MA TLAB绘制指定函数的曲线图，将所有输入的指令保存为M文件。

【实验内容】
（1）在[0,2*pi]上产生50个等距采样数据的一维数组，用两种不同的指令实现。

0:(2*pi-0)/(50-1):2*pi 或linspace(0,2*pi,50)
（2）将一维数组A=1:18，转换为2×9数组和2×3×3数组。

reshape(A,2,9)
ans =
Columns 1 through 7
1 3 5 7 9 11 13
2 4 6 8 10 12 14
Columns 8 through 9
15 17
16 18
reshape(A,2,3,3)
ans(:,:,1) =
1 3 5
2 4 6
ans(:,:,2) =
7 9 11
8 10 12 ans(:,:,3) =
13 15 17 14 16 18
（3）A=[0 2 3 4 ;1 3 5 0],B=[1 0 5 3;1 5 0 5]，计算数组A 、B 乘积，计算A&B,A|B,~A,A= =B,A>B 。

A.*B
ans=
0 0 15 12
1 15 0 0 A&B
ans =
0 0 1 1
1 1 0 0 A|B
ans =
1 1 1 1
1 1 1 1
~A
ans =
1 0 0 0
0 0 0 1
A==B
ans =
0 0 0 0
1 0 0 0
A>=B
ans =
0 1 0 1
1 0 1 0
（4）绘制y= 0.53
t e -t*t*sin(t),t=[0,pi]并标注峰值和峰值时间,添加标题y= 0.53
t e -t*t*sint ，将所有输入的指令保存为M 文件。

a=0.5 b=1/3
t=0:0.001:pi
y=a*exp(b*t)-t.*t.*sin(t) [y_max,t_max]=max(y)
t_text=['t=',num2str(t(t_max))] y_text=['y=',num2str(y_max)]
max_text=char('maximum',t_text,y_text) tit=['y=a*exp(',num2str(b),'t)-t*t*sin(t)'] hold on plot(t,y,'y.')
plot(t(t_max),y_max,'r')
text(t(t_max)+0.3,y_max+0.1,max_text)
title(tit),xlabel('t'),ylabel('y'),hold off
【实验心得与总结】
通过这次试验让我了解常用简单函数的功能，学会利用函数解决一些；数值计算和符号计算的实际问题;利用Matlab的help命令查询一些函数的功能。

利用MA TLAB可以让繁琐的
计算问题变得更加简单化，如矩阵运算等。

\
实验2：MATLAB 绘图
【实验步目的】
利用MTALAB 画墨西哥帽子，及参数方程的图像 【实验原理】
(1)二维绘图命令：plot(x,y)函数
(2)三维绘图命令中三维曲线：plot3(x,y,z), (3)利用mesh 函数画三维的网格表面的。

【实验内容】
（含参考程序、实验结果及结果分析等）
画出函数图形π100)cos(23
≤≤⎪⎩
⎪⎨⎧===t t z t
y t x 。

方程：
π100)cos(23
≤≤⎪⎩
⎪⎨⎧===t t z t
y t x 【参考程序】 >> t=0:0.1:4*pi;
>> plot3(2*cos(t),t.^3,t) 【实验结果】
画出曲面
]5.7,5.7
[
,
sin
)
,
(
2
2
2
2
-
∈
+
+
=
=x
y
x
f
z
y
x
y
x
的图像。

方程：
(,)[7.5,7.5],[7.5,7.5]
z f x y x y
=∈-∈-
【参考程序】
x = -7.5:0.5:7.5;
y = x;
[xx, yy] = meshgrid(x, y);
R = sqrt(xx.^2 + yy.^2) + eps;
z = sin(R)./R;
surf(xx, yy, z)
【实验结果】
【实验心得与总结】
Matlab的常见错误：Inner matrix dimensions must agree
1、因为在Matlab的输入变量是矩阵，参与运算的矩阵维数必须对应，矩阵相应元素的运算必须全部加dot（点），例2中方程如果这样输入：x=2*(cos(t)+t*sin(t))，就会出现该错误.
2、mesh函数是用来画三维的网格表面的。

三维空间中的一个点是用(x,y,z)来表示的，mesh就是把这些点之间用网格连接起来。

实验3：MATLAB 微积分问题的计算
【实验目的】利用MTALAB 求解二重积分、勒展开式及级数求和。

【实验原理】
1．利用int(int(f,x,a,b),y,c,d)函数求二重积分计算累次积分
⎰⎰
d c
b
a
dxdy y x f ),(
2．利用泰勒函数taylor （f,n,x,a)来求f(x,y)的n-1阶泰勒展开式
k n k k a x k a f x f )(!
)
()(1
)(-∙=∑
-=； 3.利用函数symsum(f,k,n1,n2)来求级数的和函数
∑=2
1
)(n n k k f
【实验内容】（含参考程序、实验结果及结果分析等）
求
⎰⎰+1
01
22
x xydydx
x。

【参考程序】 >> syms x y >> z=x*y;
>> f=int(int(z,y,2*x,x^2+1),x,0,1) 【实验结果】 f =1/12
将f （x ）=ln x 展开为幂为（x-2)的5阶泰勒展开式。

【参考程序】 >> syms x n;
>> f=(-1)^n*x^(n+1)/(n+1); >> symsum(f,n,1,inf) 【实验结果】
ans =log(1+x)-x
级数求和)1,1(,1)1(1
1
-∈+-+∞
=∑x n x n n n。

【参考程序】
>> syms x n;
>> f=(-1)^n*x^(n+1)/(n+1); >> symsum(f,n,1,inf) 【实验结果】
ans = log(1+x)-x 【实验心得与总结】
1、在实验过程中，要是一句程序结束后加了分号，则说明，不要求执行程序时输出执行结果；
2、在matlab 中是区别大小写的，如果N 写成n 会出现Undefined function or variable
'n'.Undefined function or variable 'n'.的错误提示.
实验4： MATLAB 优化计算
【实验目的】
掌握应用matlab 求解无约束最优化问题的方法
【实验原理与方法】 1：标准形式：
元函数
为其中n R R f X f n
R x n
→∈:)
(min
2．无约束优化问题的基本算法一．最速下降法（共轭梯度法）算法步骤： ⑴ 给定初始点n E X ∈0，允许误差0>ε,令k=0； ⑵ 计算()
k X f ∇；
⑶ 检验是否满足收敛性的判别准则： ()
ε≤∇k X f ，
若满足，则停止迭代，得点k X X ≈*，否则进行⑷； ⑷ 令()
k k X f S -∇=，从k X 出发，沿k S 进行一维搜索， 即求k λ使得： ()()
k k k k k S X f S X f λλλ+=+≥0
min ；
⑸ 令k k k k S X X λ+=+1，k=k+1返回⑵.
最速下降法是一种最基本的算法，它在最优化方法中占有重要地位.最速下降法的优点是工作量小，存储变量较少,初始点要求不高；缺点是收敛慢，最速下降法适用于寻优过程的前期迭代或作为间插步骤，当接近极值点时，宜选用别种收敛快的算法.．牛顿法算法步骤：
(1) 选定初始点n E X ∈0，给定允许误差0>ε，令k=0； (2) 求()k X f ∇,()()
1
2-∇k
X f ,检验：若()
ε<∇k X f ,则
停止迭代,k X X ≈*.否则, 转向(3)； (3) 令 ()()
k k k X f X f S ∇∇-=-12][（牛顿方向）； (4) k k k S X X +=+1,1+=k k ,转回(2).
如果f 是对称正定矩阵A 的二次函数，则用牛顿法经过一次迭代 就可达到最优点,如不是二次函数，则牛顿法不能一步达到极值点， 但由于这种函数在极值点附近和二次函数很近似,因此牛顿法的收 敛速度还是很快的.
牛顿法的收敛速度虽然较快，但要求Hessian 矩阵要可逆，要计算二阶导数和逆矩阵，就加大了计算机计算量和存储量. 【实验内容】
1. 求 f = 2x e x
sin -在0<x<8中的最小值与最大值 主程序为wliti1.m: f='2*exp(-x).*sin(x)'; fplot(f,[0,8]); %作图语句 [xmin,ymin]=fminbnd (f, 0,8) f1='-2*exp(-x).*sin(x)';
[xmax,ymax]=fminbnd (f1, 0,8) 运行结果：
xmin = 3.9270 ymin = -0.0279 xmax = 0.7854 ymax = 0.6448
2. 对边长为3米的正方形铁板，在四个角剪去相等的正方形以制成方形无盖水槽，问如 何剪法使水槽的容积最大？
先编写M 文件fun0.m 如下: function f=fun0(x) f=-(3-2*x).^2*x;
主程序为wliti2.m:
[x,fval]=fminbnd('fun0',0,1.5); xmax=x
fmax=-fval 运算结果为: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边长为0.5米时水槽的容积最大,最大容积为2立方米.
实验5： MATLAB 图论问题计算
【实验目的】
了解用Matlab 软件求解图论模型及层次分析模型的方法。

【实验内容与原理】
内容：1.某城市要建立一个消防站，为该市所属的七个区服务，如图所示．问应设在那个区，才能使它至最远区的路径最短。

2.某矿区有七个矿点，如图所示．已知各矿点每天的产矿量
)
(j v q （标在图的各顶点上）．现
要从这七个矿点选一个来建造矿厂．问应选在哪个矿点，才能使各矿点所产的矿运到选矿厂
所在地的总运力（千吨公里）最小．
原理：利用层次分析法和图论方法模型的一般概念，理解建立层次分析法和图论方法模型的一般方法，初步学会建立层次分析法和图论方法模型以解决实际问题。

【操作方法与步骤】 步骤：1.(1)用Floyd 算法求出距离矩阵D=
ν
ν⨯)(ij d ．
设剪去的正方形的边长为x ，则水槽的容积为：x x )23(2- 建立无约束优化模型为：min y=-x x )23(2
-， 0<x<1.5
解：
(2)计算在各点i
v 设立服设 施的最大服务距离
)
(i v S
}
{max )(1ij j i d v S ν
≤≤= ν ,2,1=i
(3) 求出顶点
k
v ，使
)}
({min )(1i i k v S v S ν
≤≤=
建立M 文件
a=[0 3 inf inf inf inf inf; 3 0 2 inf 18 2.5 inf; inf 2 0 6 2 inf inf; inf inf 6 0 3 inf inf; inf 18 2 3 0 4 inf; inf 2.5 inf inf 4 0 1.5; inf inf inf inf inf 1.5 0];
[D,R]=floyd(a) 点击运行 则
k
v 就是要求的建立消防站的地点．此点称为图的中心点．
2. （1）求距离阵D=
ν
ν⨯)(ij d ．
（2） 计算各顶点作为选矿厂的总运力
)
(i v m
ij
j j i d v q v m ⨯=∑=)()(1
ν
ν ,2,1=i （3）求
k
v 使
)}
({min )(1i i k v m v m ν
≤≤=，
（4）建立M 文件 a=[0 3 inf inf inf inf inf; 3 0 2 inf inf 4 inf; inf 2 0 6 2 inf inf; inf inf 6 0 1 inf inf; inf inf 2 1 0 4 inf; inf 4 inf inf 4 0 1.5; inf inf inf inf inf 1.5 0]; [D,R]=floyd(a) q=[3,2,7,1,6,1,4]; m=0; for i=1:7 for j=1:7
m=m+q(i)*D(i,j); end m m=0;
end
点击运行 （5）则
k
v 就是选矿厂应设之矿点．此点称为图G 的重心或中位点．
【实验结果与分析】 实验结果与分析：1.
S(v1)=10, S(v2)=7, S(v3)=6, S(v4)=8.5, S(v5)=7, S(v6)=7, S(v7)=8.5 S(v3)=6,故应将消防站设在v3处。

2.
由上述公式可得：
m(v1)=38.5*3=115.5, m(v2)=23.5*2=47, m(v3)=23.5*7=164.5, m(v4)=28.5, m(v5)=23.5*6=141, m(v6)=27.5, m(v7)=35*4=140
再求其中的最小值，m(v6)=27.5，则6v 就是选矿厂应设之矿点．
实验6：MATLAB 计算机模拟计算
【实验目的】
学会用数学软件matlab 和蒙特卡洛方法估计积分值，并于其中应用概率论中的概率密度等知识点。

【实验问题】
估计积分值，并对误差进行估计。

【实验要求】
针对要估计的积分选择适当的概率分布设计蒙特卡洛估计算法； 利用计算机产生所选分布的随机样本的估计积分值； 通过计算平均误差对估计结果进行评价。

【实验过程分析】
（x 为运行结果平均值，y 为样本方差）
估计程序如下： clc; clear; m=10; n=10000; d=0; e=0;
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫
⎝⎛=05.15.55.65.75.55.85.10456475.54012475.65
103585.76230255.54452035.8778530D
for i=1:m
d=0;
a=rand(1,n);
for j=1:n
b=a(j)+2;
c=b^2;
e(i)=d+c/n;
d=e(i);
end
fprintf('e=%.8f\n',e(i))
end
p=sum(e)/m;
for j=1:m;
s(j)=(e(j)-p)^2;
end
q=sum(s);
fprintf('x=%.8f\ny=%.8f\n',p,q); 结果为：
e=6.34879520
e=6.34068140
e=6.35081124
e=6.31353632
e=6.35586630
e=6.33058791
e=6.32419121
e=6.33707454
e=6.30357011
e=6.35063255
x=6.33557468
y=0.00304243
2.估计程序变动如下：
b=a(j)*3;
c=b*sin(b)*3;
结果为：
e=3.12211717
e=3.11373037
e=3.07484948
e=3.08660758
e=3.10052243
e=3.10475698
e=3.13762746
e=3.16481618
e=3.11552000
e=3.09615989
x=3.11167076
y=0.00669223
3.估计程序变动如下：b=exp(-a(j)^2/2);
c=b/n*(2*pi)^0.5;
e(i)=d+c/2;
结果为：
e=0.88617655
e=0.88538972
e=0.88635209
e=0.88575809
e=0.88653705
e=0.88606366
e=0.88634011
e=0.88613926
e=0.88573325
e=0.88644154
x=0.88609313
y=0.00000132
4.估计程序变动如下：b=exp(a(j)^2);
c=b/n;
e(i)=d+c;
结果为：
e=1.46211146
e=1.46154792
e=1.46327379
e=1.46256348
e=1.46318297
e=1.46235828
e=1.46241378
e=1.46316145
e=1.46203052
e=1.46280489
x=1.46254485
y=0.00000320
5.估计程序变动如下：
b=a(i)*4;
c=1/((1+b^2)^0.5);
e(i)=d+c*4/n;
结果为：
e=1.98511173
e=1.02167881
e=1.26713031
e=0.98770837
e=1.14216662
e=1.75642022
e=1.97055988
e=1.96227794
e=1.83229787
e=1.06190231
x=1.49872541
y=1.90228258
【实验结果与分析】
通过对实验所得平均值与真实值的比较，可以看出实验结果与真实值相比非常接近，而且样本方差很小，从而说明概率分布的选取比较适当，计算机实验很准确。

实验7：MATLAB与马尔科夫预测模型
【实验目的】
基于matlab编程应用马尔可夫预测模型
【实验原理】
马尔可夫通过实践认为：世界上无论是社会领域还是自然领域，有一类事物的变化过程只与事物的近期状态有关，与事物的过去状态无关，这类事物的性质称为无后效性。

例如，事物π，从初始状态π(0)起，变动一次后为π(1)，变动n次后为π(n)，则π(n)仅与π(n-1)有关，与n-1以后的各次变动无关。

马尔可夫链：如果n个连续变动的事物，在变动的过程中，其中任一次变动的结果都具有无后效性，那么，这n个连续变动事物的集合，就叫做马尔可夫链，这类事物的演变过程就叫做马尔可夫过程。

【实验内容】
1.农业收成变化预测
考虑某地区农业收成变化的三个状态，即“丰收”、“平收”和“欠收”。

记E1为“丰收”状态，E2为“平收”状态，E3为“欠收”状态。

下表给出了该地区1965～2004年期间农业收成的状态变化情况。

试计算该地区农业收成变化的状态转移概率矩阵，并进行预测。

使用matlab实现如下：
P=[0.2000 0.4667 0.3333;0.5385 0.1538 0.3077;0.3636 0.4545 0.1818];
% 读入状态转移概率矩阵
x=[0,1,0];
% 读入初始状态概率向量(2004年的农业收成状态)
for i=1:11 % 预测今后11 年(2005－2015)的农业收成状态
y= x*P^i
end
运行结果如下：
y =0.5385 0.1538 0.3077
y =0.3024 0.4148 0.2827
y =0.3867 0.3334 0.2798
y =0.3586 0.3589 0.2823
y =0.3677 0.3509 0.2813
y =0.3648 0.3534 0.2817
y =0.3657 0.3526 0.2815
y =0.3654 0.3529 0.2816
y =0.3655 0.3528 0.2815
y =0.3654 0.3528 0.2815
y =0.3654 0.3528 0.2815
2.市场占有率预测
某厂对某产品的市场占有率和销售情况进行了调查：一月份共销售了50万件，其中普通、一级、特级品分别为35、10、5万件。

二月份中，一月份买普通品的顾客25%的顾客转买一级品，8%的顾客转买特级品；一月份买一级品的顾客10%转买特级品，3%转买普通品；一月份买特级品的顾客2%买普通品，15%转买一级品。

请预测以后月份各个等级产品的市
场占有率。

由所给的资料可知
35105
[,,][0.7,0.2,0.1]
505050
π==
使用matlab实现如下：
P=[0.67 0.25 0.08;0.03 0.87 0.1;0.02 0.15 0.83];
%读入状态转移矩阵
x=[0.7 0.2 0.1];
% 读入初始状态概率向量（一月份各等级产品的市场占有率）
for i=1:11 % 预测今年剩余11个月各产品等级的市场占有率
y= x*P^i
end
运行结果如下：
y =0.4770 0.3640 0.1590
y =0.3337 0.4598 0.2065
y =0.2415 0.5144 0.2441
y =0.1821 0.5445 0.2734
y =0.1438 0.5603 0.2959
y =0.1191 0.5678 0.3131
y =0.1031 0.5707 0.3262
y =0.0927 0.5712 0.3361
y =0.0860 0.5705 0.3435
y =0.0816 0.5694 0.3490
y =0.0787 0.5681 0.3532
结论：顾客对普通品的需求有减少的趋势，对一级品和特级品的需求有增加的趋势，因此，可以调整相应等级产品的产量。

实验8：基于MATLAB的灰色预测模型
【实验目的】
实验目的：掌握灰色预测模型及其应用
基本内容：灰色预测模型的提出，建模以及实现代码。

【实验内容】
灰色系统( Grey System)理论是我国著名学者邓聚龙教授20世纪80年代初创立的一种兼备软硬科学特性的新理论。

该理论将信息完全明确的系统定义为白色系统，将信息完全不明确的系统定义为黑色系统，将信息部分明确、部分不明确的系统定义为灰色系统。

由于客观世界中，诸如工程技术、社会、经济、农业、环境、军事等许多领域，大量存在着信息不完全的情况。

要么系统因素或参数不完全明确，因素关系不完全清楚；要么系统结构不完全知道，系统的作用原理不完全明了等，从而使得客观实际问题需要用灰色系统理论来解决。

灰色预测是应用灰色模型GM( 1，1) 对灰色系统进行分析、建模、求解、预测的过程。

由于灰色建模理论应用数据生成手段，弱化了系统的随机性，使紊乱的原始序列呈现某种规律，规律不明显的变得较为明显，建模后还能进行残差辨识，即使较少的历史数据，任意随机分布，也能得到较高的预测精度。

因此，灰色预测在社会经济、管理决策、农业规划、气象生态等各个部门和行业都得到了广泛的应用。

一、GM（1，1）模型建立
设有k个原始非负样本序列为揭示系统的客观规律，灰色系统理论采用了独特的数据预处理方式，对序列{}进行一阶累加生成，即AGO生成，
由此得生存数列：
GM(1,1)模型的原始形式为：
为的紧邻值生存序列
其中
GM(1,1)模型的基本形式：
若=为参数列，且
,
则GM(1,1)模型的最小二乘估计参数列满足
关于的白化方程也叫影子方程为：
定理：白化方程的解也称时间响应函数为：
GM(1,1)模型的时间响应序列为
实际预测值为
二、模型检验
为确保所建立的GM(1,1)模型有较高的预测精度，还需要进行以下检验
（1）求出及之残差e(k)、相对误差和平均相对误差：
,，.
（2）求出原始数据平均值x:
.
三、残差修正模型
记，当然我们知道并不一定全为负或者全为正，这时我们令
,同时令,则是一个非负序列，我们可以用方法来建立它的GM模型，求解，得出其预测值,而后还原残差预测值,最后用修正原来的预测值，得到修正后的预测值.
四、实现代码
clear all
X0=input('请输入序列矩阵 ');%输入数据请用如例所示形式：[48.7 57.17 68.76 92.15]或者 [43823050,44649620,45793750 46613050 47715440 48526270 49713050]，该向量为原始向量X0n=length(X0);
for i=2:n %开始进行建模可行性分析
Q(i)=X0(i-1)/X0(i);
end
Q(1)=[];
ma=max(Q);
mi=min(Q);
if ma>exp(2/(n+1))
disp(['序列无法进行灰色预测']);
return
elseif mi<exp(-2/(n+1))
disp(['序列无法进行灰色预测']);
return
else
disp(['序列可以进行灰色预测']);
end
clear Q ma mi %检验结束
X1=cumsum(X0);%累加生成算子向量X1
Z1=ones(n-1,2);
for i=1:(n-1)
Z1(i,1)=-(X1(i)+X1(i+1))/2;
Z1(i,2)=1;%均值生成算子Z1
end
Z1T=Z1';%均值生成算子矩阵Z1的转置Z1T
for j=1:n-1
Y(j)=X0(j+1);
end
YT=Y';
A=inv(Z1T*Z1)*Z1T*YT;%最小二乘估计计算参数a、u
a=A(1);%Z1参数a
u=A(2);%系统给定参数u
t=u/a;
t_test=input('请输入需要预测个数：');
i=1:t_test+n;
X1S(i+1)=(X0(1)-t).*exp(-a.*i)+t;%计算时间响应序列，得出估计累加向量X1S
X1S(1)=X0(1);
X0S(1)=X0(1);
for j=n+t_test:-1:2
X0S(j)=X1S(j)-X1S(j-1);%计算X1S的逆累加向量X0S,还原X0得到估计值
end
for i=1:n
Q(i)=X0S(i)-X0(i);%求残差
E(i)=abs(Q(i))/X0(i);%求相对误差
end
AVG=sum(E)/(n-1);%求平均相对误差
av=input('请输入允许的平均相对误差');%输入如0.1，或0.05等形式，不要用5%这类形式if AVG>=av;%如果平均相对误差大于av%，则进入残差GM模型
clear cn Q1 CZ1 CZ1T CY CYT CA ca cu ct Q1S
cn=length(Q);
Q1=cumsum(Q);%累加生成算子向量Q1
CZ1=ones(cn-1,2);
for i=1:(n-1)
CZ1(i,1)=-(Q1(i)+Q1(i+1))/2;
CZ1(i,2)=1;%均值生成算子CZ1
end
CZ1T=CZ1';%均值生成算子矩阵CZ1的转置CZ1T
for j=1:cn-1
CY(j)=Q(j+1);
end
CYT=Y';
CA=inv(CZ1T*CZ1)*CZ1T*CYT;%最小二乘估计计算参数ca、cu
ca=CA(1);%CZ1参数a
cu=CA(2);%系统给定参数cu
ct=cu/ca;
i=1:t_test+cn;
Q1S(i+1)=(Q(1)-ct).*exp(-ca.*i)+ct;%计算时间响应序列，得出估计累加向量Q1S X1S=X1S+Q1S;%将残差拟合值加入，提高精度
for j=n+t_test:-1:2
X0S(j)=X1S(j)-X1S(j-1);%计算X1S的逆累加向量X0S
end
clear av
for i=1:cn
Q(i)=X0S(i)-X0(i);%求残差
E(i)=abs(Q(i))/X0(i);%求相对误差
end
AVG=sum(E)/(n-1);%求平均相对误差
end
x=1:n;
xs=2:n+t_test;
yn=X0S(2:n+t_test);
plot(x,X0,'^r',xs,yn,'*-b');%作图
disp(['百分平均相对误差为：',num2str(AVG*100),'%']);
disp(['拟合值为： ',num2str(X0S(1:n+t_test))]);
disp(A);
【实验目总结】
1、灰色预测模型的应用范围有限，在使用时不能超过其能预测的范围
2、灰色预测模型的模型建立，以及检验非常重要，可以减少预测带来的误差
实验9：基于MATLAB的模糊聚类分析
【实验目的】
1. 加强对MA TLAB 软件使用的能力。

2. 加强对模糊数学学习的理解。

3. 基于MATLAB 对数据挖掘中使用模糊数学进行聚类分析
学会使用MATLAB 中的模糊工具箱。

【实验方案】
1. 实验原理
首先，我们要了解一般聚类分析分为三个步骤：（1）数据标准化，（2）标定，（3）聚类。

然后我们在来了解下什么是模糊数学：
（1） 模糊数学是1956年，美国加利福尼亚大学控制论专家扎德（L.A.Zadeh ）教授
提出的。

是一门研究和处理模糊性现象的数学方法。

（2） 对于有限论域 构造映射A(x)：U →[0,1]，确定U 上的模糊子集A ，映射A(x)称为A 的隶属函数，它表示x 对A 的隶属程度。

（3） 模糊集的基本运算：
相等：A = B ⇔ A (x ) = B (x )；
包含：A ⊆B ⇔ A(x )≤B(x )；
并：A ∪B 的隶属函数为
(A ∪B )(x )=A (x )∨B (x )=max{A (x )，B (x ）}；
交：A ∩B 的隶属函数为
(A ∩B )(x )=A (x )∧B (x )=min{A (x )，B (x ）}；
余：Ac 的隶属函数为
Ac (x ) = 1- A (x ).
2. 基于MATLAB 利用模糊数学实现聚类分析。

（1） 数据标准化。

设论域X = {x 1, x 2, …, xn }为被分类对象,每个对象又由m 个指标表示其性状: xi = { xi 1, xi 2, …, xim }, i = 1, 2, …, n 。

根据模糊矩阵的要求，将数据压缩到[0,1]区间上。

通常需要作如下几种变换： ①平移 • 标准差变换
②平移 • 极差变换
③对数变换
取对数以缩小变量间的数量级。

（2） 标定
相似性度量，又称标定，就是根据实际情况，按一定准则或某一种方法，给论域X 中的元素两两之间都赋以[0,1]内的一个数，称为相似系数。

它的大小表征两个元素彼此接近或相似的程度。

12{,}n U x x x =),...,2,1;,...,2,1(m k n i s x x x k k ik ik ==-='∑∑==-==n i k ik k n i ik k x x n s x n x 12
1)(1,1其中}
1|min{}1|max{}1|min{n i x n i x n i x x x ik ik ik ik ik ≤≤-≤≤≤≤-='')
,,2,1;,,2,1(lg m k n i x x ik ik ==='
用rij 表示元素xi 与xj 的相似系数，其中： xi = { xi1, xi2, …, xim}, i = 1, 2, …, n xj = { xj1, xj 2, …, xjm}, j = 1, 2, …, n
rij ∈[0,1] ①数量积法
②夹角余弦法
③相关系数法
④最大最小法。

一般海明距离法
一般欧式距离
（3） 聚类
① 传递闭包法
标定得到的模糊矩阵，一般仅具有自反性、对称性，故可用传递闭包运算将模糊相似矩阵R ，改造成模糊等价矩阵t(R)后进行分类。

② 布尔矩阵法 ③ 最大树法
3. 通过在网上查资料学习对MA TLAB 中模糊工具箱的使用。

【实验记录】 1. 数据的标准化 x =
)
(max ,1,111∑
∑
=≠=⋅=⎪⎩⎪⎨⎧≠⋅==m
k jk ik j i m
k jk ik ij
x x M j i x x M j i r 其中∑
∑∑
====
m k jk m k ik m
k jk
ik ij x x x x r 1
2121
∑
∑∑
===----=
m k j jk m k i ik m
k j jk i ik ij x x x x x x x x r 1
2
12
1
)()(|
|||∑∑
==∨∧=m k jk ik
m
k jk ik ij x x
x x r 1
1
)
()
(∑
=-=
m k jk ik j i x x x x d 1
|
|),(R R R ⊆
0.2973 0.1351 0.1712 0.3964
0.2703 0.1532 0.1622 0.4144
0.2703 0.0631 0.2162 0.4505
0.4234 0.2883 0.1081 0.1802
0.2342 0.1081 0.2342 0.4234
0.3514 0.1261 0.1261 0.3964
0.3514 0.1892 0.0991 0.3604
0.2793 0.1892 0.1622 0.3694
0.2072 0.1532 0.2072 0.4324
0.1818 0.1364 0.2727 0.4091
0.3545 0.5000 0.0455 0.1000
0.3273 0.5000 0.0273 0.1455
0.2545 0.5182 0.1000 0.1273
0.3000 0.5000 0.0818 0.1182
0.2909 0.6455 0 0.0636
0.3636 0.4636 0.0818 0.0909
0.3545 0.2636 0.2455 0.1364
0.2909 0.5000 0.1182 0.0909
0.2182 0.5636 0.1455 0.0727
0.2000 0.5636 0.1727 0.0636
0.2743 0.3628 0.1947 0.1681
0.2885 0.2212 0.2404 0.2500
0.1765 0.1863 0.2549 0.3824
0.2087 0.4087 0.1913 0.1913
0.2476 0.2190 0.2286 0.3048
0.2193 0.3860 0.2105 0.1842
0.2308 0.2308 0.2019 0.3365
0.2564 0.4444 0.1453 0.1538
0.1485 0.1881 0.2178 0.4455
0.2897 0.2523 0.2430 0.2150
0.2411 0.3571 0.1786 0.2232
0.1743 0.3303 0.2294 0.2661
0.2703 0.3333 0.1892 0.2072
0.2353 0.1667 0.2353 0.3627
0.2427 0.2039 0.2136 0.3398
0.2286 0.2095 0.3048 0.2571
0.2136 0.2039 0.2524 0.3301
0.2222 0.4359 0.1709 0.1709
0.2736 0.2358 0.2830 0.2075
0.1983 0.4310 0.1983 0.1724
2. 标定
得出40*40的相似矩阵，因位置关系在这里只展示7*7
1.0000 0.9779 0.9398 0.8267 0.9440 0.9590 0.9359
0.9779 1.0000 0.9359 0.8179 0.9465 0.9454 0.9298 0.9398 0.9359 1.0000 0.7699 0.9618 0.9152 0.8785 0.8267 0.8179 0.7699 1.0000 0.7814 0.8382 0.8741 0.9440 0.9465 0.9618 0.7814 1.0000 0.9061 0.8809 0.9590 0.9454 0.9152 0.8382 0.9061 1.0000 0.9553 0.9359 0.9298 0.8785 0.8741 0.8809 0.9553 1.0000 3. 聚类
首先通过 得到等价矩阵，同理如下只展示7*7
1.0000 0.9779 0.9582 0.9067 0.9582 0.9590 0.9553 0.9779 1.0000 0.9582 0.9067 0.9582 0.9590 0.9553 0.9582 0.9582 1.0000 0.9067 0.9618 0.9582 0.9553 0.9067 0.9067 0.9067 1.0000 0.9067 0.9067 0.9067 0.9582 0.9582 0.9618 0.9067 1.0000 0.9582 0.9553 0.9590 0.9590 0.9582 0.9067 0.9582 1.0000 0.9553 0.9553 0.9553 0.9553 0.9067 0.9553 0.9553 1.0000 然后通过取不同的λ来获得等价矩阵的截矩阵来获得聚类。

因为当λ=1时，分为40类，所以因为篇幅限制，只举几个例子： λ=0.9730 C =
1 2 0 0 0 3 0 0 0 0 4 0 0 0 0 5 0 0 0 0 6 0 0 0 0 7 0 0 0 0 8 0 0 0 0 9 0 0 0 0 10 0 0 0 0 11 0 0 0 0 12 0 0 0 0 13 0 0 0 0 14 18 0 0 0 15 0 0 0 0 16 0 0 0 0 17 0 0 0 0 19 20 0 0 0 21 0 0 0 0 22 0 0 0 0 23 0 0 0 0 24 26 28 38 40 25 27 35 0 0 29 0 0 0 0
R R R
30 39 0 0 0
31 33 0 0 0
32 0 0 0 0
34 0 0 0 0
36 0 0 0 0
37 0 0 0 0
λ=0.9600
C =
1 2 8 0 0 0 0 0 0 0
3 5 9 29 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0 0
10 22 23 25 27 30 34 35 37 39
11 12 13 14 16 18 0 0 0 0
15 0 0 0 0 0 0 0 0 0
17 0 0 0 0 0 0 0 0 0
19 20 0 0 0 0 0 0 0 0
21 24 26 28 31 33 38 40 0 0
32 0 0 0 0 0 0 0 0 0
36 0 0 0 0 0 0 0 0 0
λ=0.9582
C =
1 2 3 5 6 8 9 10 22 23 25 27 29 30 34 35 37 39
4
7
11 12 13 14 16 18
15
17
19 20
21 24 26 28 31 33 38 40
32
36
λ=0.9024
C =
Columns 1 through 29
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
Columns 30 through 40
30 31 32 33 34 35 36 37 38 39 40
通过以上过程可以看到整个聚类的过程。

3.关于模糊工具箱的应用。

在这里根据研究需要，仅演示模糊c均值聚类。

上图为随机一个100*2的数组对其做均值聚类。

以上为一个100*4数组对其做均值聚类。

此为与上面聚类相同的四维数据进行的聚类。

【实验总结】
模糊聚类分析中存在大量的矩阵运算，通过MATLAB建立M文件，可以方便计算并快速
得到聚类分析的结果。

模糊工具箱的应用使现实生活中的模糊问题更加可视化。