资源描述
第4章 线性代数模型
4.1 行列式与矩阵
本节案例主要涉及线性代数中矩阵与方阵的行列式等概念,通过案例建立数学模型,加深对行列式、矩阵及矩阵运算等相关知识的进一步理解以及了解这些概念的实际应用。
4.1.1 过定点的多项式方程的行列式
1.问题提出
求通过空间中三个点,,的平面方程。
2.模型建立与求解
已知三个点可以确定一个平面,设平面方程为,而三个点在这个平面上,所以它们均满足这个平面方程,因而有
这是一个以为未知量的齐次线性方程组,且不全为0,说明该齐次线性方程组必有非零解,于是系数行列式等于零,即
,
从而得到平面方程为。
计算的MATLAB程序如下:
clc, clear, syms x y z
D=[x,y,z,1;1,2,3,1;3,5,6,1;2,2,4,1];
s=det(D)
3.模型拓展
对于次多项式,其系数为,可由其曲线上个横坐标互不相同的点所唯一确定。
因为个点满足这个多项式,则有
这是一个含有个方程,以为个未知量的线性方程组,其系数行列式为
. (4.1)
这是一个范德蒙行列式。当互不相同时,,由克拉默法则,可解出唯一的,所以次多项式可由其曲线上的个横坐标互不相同的点所唯一确定。
该多项式方程为
.
4.1.2 循环比赛名次模型
例4.1 图4.1给出了6支球队的比赛结果,即1队战胜2,4,5,6队,而输给了3队;5队战胜3,6队,而输给1,2,4队等等。试给出6支球队的一个排名顺序。
图4.1 球队的比赛结果
解 若胜一场记1分,则6支球队的比赛得分情况见表4.1。
表4.1 6支球队的比赛得分情况
球队
1
2
3
4
5
6
1
-
1
0
1
1
1
2
0
-
0
1
1
1
3
1
1
-
1
0
0
4
0
0
0
-
1
1
5
0
0
1
0
-
1
6
0
0
1
0
0
-
通常的方法是计算各队的得分,从高到低排名次,则上面6支球队的得分是
.
按得分计算,球队1排名第一,球队2,3并列第二,……。
这种方法有两个缺点:第一是不能区分球队2与球队3的高低(因为不计算小分);第二是没有考虑所打败对手的强弱,因为不论你战胜的是强队还是弱队,都得1分。那么,如何确定一种方法,使排名更合理呢?
如何考虑对手的强弱,首先要回答一个问题—用什么指标刻画对手的强弱?从直观上讲,所谓对手强,就是对手战胜的队多,对手弱就是对手战胜的队少,那么将对手的得分也记录到自己队的得分中,即所谓的二级得分。这样考虑问题要比直接考虑得分更合理,这里计算一下二级得分
.
按这种计算方法,球队3排名第一,球队1排名第二,……。
这看起来是合理的,但所说的问题并没有根本解决。因为这里只考虑每个队它的对手得分,而没有考虑每个队它对手的对手得分,也就是说,评判对手强弱的计算公式不合理。为了合理的解决这个问题,需要考虑每个队对手的对手得分,即三级得分。那么同样的问题,还有四级得分、五级得分……。看一下各队各级的得分。
,
,
,
.
看来各队名次的排列有点波动,例如球队3和球队1竞争第一名的情况就是这样。如何确定他们的名次呢?再对各队的得分进行分析。将各队比赛的得分情况构造矩阵
同时构造列向量,它可以看成各队的实力,比赛前,认为6个队是相同的。那么第一级得分为
,
也可以被看成各队的实力,比赛以后就不同了。而二级以后的得分如下
,
,
……………………
,
最合理的情况应考虑极限情况,即
.
式中,是矩阵的最大特征值,是对应于的特征向量。可以证明在满足一定的条件下,上述极限是存在的。
因此,名次排列问题就转化为求得分矩阵的最大特征值和对应的特征向量。对于前面的问题,其最大特征值为,相应的特征向量为
,
即第一名到第六名的名次排列为球队1、球队3、球队2、球队5、球队4、球队6。
计算的MATLAB程序如下:
clc, clear
a=load('gtable4_1.txt'); %表格数据中的‘-’替换成0,保存到纯文本文件中
e=ones(6,1); s(:,1)=a*e;
for i=2:6
s(:,i)=a*s(:,i-1); %各级得分保存在矩阵的一列
end
s %显示各级得分
[vec,val]=eigs(a,1) %求最大值及对应的特征向量
vec2=vec/sum(vec) %特征向量归一化
4.1.3 平面图形的几何变换
随着计算机科学技术的发展,计算机图形学的应用领域越来越广,如仿真设计、效果图制作、动画片制作、电子游戏开发等。图形的几何变换包括图形的平移、旋转、放缩等,是计算机图形学中经常遇到的问题。这里只讨论平面图形的几何变换。
1.齐次坐标
平面图形的旋转和放缩都很容易用矩阵乘法实现,但是平面图形的平移并不是线性运算,不能用矩阵乘法表示。可以使用齐次坐标使平移、旋转、放缩能统一用矩阵乘法来实现。
齐次坐标就是一个维向量的维向量表示。例如,二维坐标点的齐次坐标为。
二维坐标与齐次坐标是一对多的关系。通常都采用规格化的齐次坐标,即取。的规格化齐次坐标为。
齐次坐标的几何意义:可理解为在三维空间上第三维为常数的一平面上的二维向量。
齐次坐标的作用是将各种变换用阶数统一的矩阵来表示。提供了用矩阵运算实现图形变换,或者把二维、三维甚至高维空间上的一个点从一个坐标系变换到另一坐标系的有效方法。
例4.2 平移变换可以用齐次坐标写成,这个变换可用矩阵乘法实现,
.
例4.3 旋转变换,这个变换可用矩阵乘法实现,
.
例4.4 放缩变换,这个变换可用矩阵乘法实现,
.
由上述几个例题可以看出,中的任何线性变换都可以用分块矩阵乘以齐次坐标实现,其中是2阶方阵,是列向量。这样,只要把平面图形上点的齐次坐标写成列向量,平面图形的每一次几何变换,都可通过左乘一个3阶变换矩阵来实现。
例4.5 求出矩阵,对应于先乘以3的倍乘变换,然后逆时针旋转,最后对图形的每个点的坐标加上做平移。
解 当时,。我们有
,
所以复合变换的矩阵为
.
几何变换把坐标变换为坐标,记作
,
具体数学表达式为
(4.2)
写成矩阵形式
. (4.3)
2.图像的空间变换
一幅图像经过采样和量化后便可以得到一幅数字图像。通常可以用一个矩阵来表示。
图4.2 数字图像的矩阵表示
一幅数字图像在MATLAB中可以很自然地表示成矩阵
其中,,。
矩阵中的元素称作像素。每一个像素都有和两个坐标,表示其在图像中的位置。另外还有一个值,称灰度值,对应于原始模拟图像在该点处的亮度。量化后的灰度值,代表了相应的色彩浓淡程度,以256色灰度等级的数字图像为例,一般由8位,即一个字节表示灰度值,由0-255对应于由黑到白的颜色变化。对只有黑白二值采用一个比特表示的特定二值图像,就可以用0和1来表示黑白二色。
在MATLAB的图像处理工具箱中提供了一个专门的函数imwarp,用户可以定义参数实现多种类型的空间变换,包括仿射变换(如平移、缩放、旋转、剪切)、投影变换等。函数imwarp具体的调用格式如下。
B=imwarp(A,tform):该函数中A是待变换的图像矩阵;tform表示为执行空间变换的所有参数的结构体;B为按照tform参数变换后的图像矩阵。
在MATLAB中利用函数imwarp实现图像的空间变换时,都需要先定义空间变换的参数。对于空间变换参数的定义,MATLAB也提供了相应的函数affine2d,affined3d,fitgeotrans等函数,它们的作用是创建进行空间变换的参数结构体。affine2d的具体调用方式如下。
tform=affine2d(C):该函数返回一个维的仿射变换参数结构体tform,输入参数C是一个的矩阵。
用户结合使用函数affine2d和函数imwarp,就可以灵活实现图像的线性变换,而变换的结果和变换参数结构体密切相关。以二维仿射变换为例,原图像为,变换后图像为,仿射变换中原图像中某个像素点坐标和变换后该像素点坐标满足关系式,写成矩阵形式即满足。
例4.6 利用函数imwarp,实现图像的旋转和缩放。
计算的MATLAB程序如下:
clc, clear, close all
a=imread('rice.png'); %MATLAB工具箱的图像文件
tf1 = affine2d([cosd(30),-sind(30),0;sind(30),cosd(30),0;0 0 1]) %逆时针旋转30度
ta1=imwarp(a,tf1); %实现图形旋转
tf2=affine2d([5 0 0;0 3 0;0 0 1]) %缩放结构体
ta2=imwarp(a,tf2); %实现图像缩放
subplot(131), imshow(a), subplot(132), imshow(ta1), subplot(133),imshow(ta2)
原图像及变换后的图像见图4.3。
(a)原图像 (b)旋转后的图像 (c)缩放后的图像
图4.3 原图像及变换后的图像
4.1.4 一种矩阵密码问题
密码学在经济和军事方面起着极其重要的作用。现代密码学涉及很多高深的数学知识。密码学中将信息代码称为密码,尚未转换成密码的文字信息称为明文,由密码表示的信息称为密文。从明文到密文的过程称为加密,反之称为解密。1929年,Hill(希尔)通过线性变换对传输信息进行加密处理,提出了密码史上有重要地位的希尔加密算法。下面只介绍希尔加密算法的基本思想。
1.古典密码的基本理论
对于正整数,记集合。
定义4.1 对于一个元素属于集合的阶方阵,若存在一个元素属于集合的方阵,使得
,
称为模可逆,为的模逆矩阵,记为。
的意义是,每一个元素减去的整数倍后,可以化成单位矩阵。例如
.
定义4.2 对的一个整数,若存在的一个整数,使得,称为的模倒数或乘法逆,记作。
可以证明,如果与无公共素数因子,则有唯一的模倒数(素数是指除了1与自身外,不能被其他正整数整除的正整数),反之亦然。例如()。利用这点,可以证明下述定理。
定理4.1 元素属于的方阵模可逆的充要条件是,和没有公共素数因子,即和互素。
显然,所选加密矩阵必须符合该命题的条件。
2.Hill密码的数学模型
一般的加密过程是这样的:
明文加密器密文普通信道解密器明文,
其中的“普通信道解密器”这个环节容易被敌方截获并加以分析。
在这个过程中,运用的数学手段是矩阵运算,加密过程的具体步骤如下:
(1)根据明文字母的表值,将明文信息用数字表示,设明文信息只需要26个拼音大写字母(也可以不止26个,如还有小写字母、数字、标点符号等),通信双方给出这26个字母表值见表4.2。
表4.2 明文字母的表值
A
B
C
D
E
F
G
H
I
J
K
L
M
1
2
3
4
5
6
7
8
9
10
11
12
13
N
O
P
Q
R
S
T
U
V
W
X
Y
Z
14
15
16
17
18
19
20
21
22
23
24
25
0
(2)选择一个三阶可逆整数方阵,称为密码的加密矩阵,它是这个加密体制的“密钥”(是加密的关键,仅通信双方掌握)。
(3)将明文字母逐对分组。密码的加密矩阵为三阶矩阵,则明文字母每3个一组(可以推广到密码,则个明文字母为一组)。若最后一组不足3个字母,则补充没有实际意义的哑字母,这样使每一组都由3个明文字母组成。查出每个明文字母的表值,构成一个三维列向量。
(4)乘以,得一新的三维列向量,由的三个分量反查字母表值得到的三个字母即为密文字母。
以上4步即为密码的加密过程。
解密过程,即为上述过程的逆过程。
例4.7 明文为“ACTIONS”,,求这段明文的密文。
解 将明文相邻字母每3个分为一组:ACT ION SXX,最后两个字母X为哑字母,无实际意义。查表4.2得到每组的表值,并构造三维列向量
,
将上述3个向量左乘矩阵,得到3个三维列向量
,
作模26运算(每个元素都加减26的整数倍,使其化为0~25的一个整数)得
,
反查表4.2得到每对表值对应的字母为:DYW XUG QHW,这就得到了密文“DYWXUGQ”。
计算的MATLAB程序如下
clc, clear
s='ACTIONS';
s=[s,'XX'] %补充两个哑字母‘X’
L=length(s); %计算字符总数
num=double(s)-64; %字母编码
num=mod(num,26) %mod26,变换Z的编码
mm=reshape(num,[3,L/3]); %把行向量变成三行的矩阵
A=[1,1,0;2,1,1;3,2,2]; %输入密钥矩阵
mw=A*mm %求密文的编码值
mw=mod(mw,26) %mod26
mw(mw==0)=26; %变换Z的编码值
mw=reshape(mw,[1,L])+64; %变换到字母的ASCII码值
mwzf=char(mw) %转换成密文的字符
mwzf(end-1:end)=[] %删除最后两个字符
例4.8 甲方收到与之有秘密通信往来的乙方的一个密文信息,密文内容:OWTRBMOQIOC,按照甲方与乙方的约定,他们之间的密文通信采用密码,密钥为三阶矩阵,问这段密文的原文是什么?
解 所选择的明文字母共26个,,26的素数因子为2和13,所以上的方阵可逆的充要条件为不能被2和13整除。若满足上述定理4.1的条件,不难验证,其中,为的伴随矩阵,是的倒数。显然,是中的数。中有模26倒数的整数及其倒数可见表4.3。
表4.3 模26倒数表
1
3
5
7
9
11
15
17
19
21
23
25
1
9
21
15
3
19
7
23
11
5
17
25
模26倒数表4.3可用下列MATLAB程序求得。
clc, clear
m=26;
for a=1:m
for i=1:m
if mod(a*i,m)==1
fprintf('The Inverse (mod %d) of number:%d is:%d\n',m,a,i)
end
end
end
利用表4.3可以反演求出如下
密文“OWTRBMOQIOC”中总共11个字符,为了凑成4组字符,我们在密文末尾再加一个字符“X”。利用可以把密文变换成明文。
,,
,,
即得到明文:MATHEMATICS。
计算的MATLAB程序如下:
clc, clear, m=26;
T=[1 3 5 7 9 11 15 17 19 21 23 25
1 9 21 15 3 19 7 23 11 5 17 25]; %模26的倒数表
a=[1,2,0;0,3,1;0,0,1]; %输入密钥矩阵
ad=det(a) %计算对应的行列式值
ind=find(T(1,:)==ad)
ai=T(2,ind) %求ad模26的倒数
B=mod(ai*det(a)*inv(a),26) %求a的模26逆阵
s='OWTRBMOQIOC';
L=length(s);
if mod(L,3)==1
s=[s,'XX']; %如果字符是奇数个,在最后补一个哑字母
L=L+2; fprintf('添加了2个字符\n')
elseif mod(L,3)==2
s=[s,'X'];
L=L+1; fprintf('添加了1个字符\n')
end
jm=double(s)-64; %求字母对应的编码
jm(jm==26)=0; %如果存在Z,把Z的编码改成0
jm2=reshape(jm,[3,L/3]) %把行向量变成三行的矩阵
mjm=mod(B*jm2,26) %求明文的编码值
mjm(mjm==0)=26; %变换Z的编码值
bm=reshape(mjm,[1,L])+64; %变换到字母的ASCII码值
mzf=char(bm) %转换成明文的字符
mzf(end)=[] %删除最后一个添加的字符
4.2 线性方程组
本节案例主要涉及线性代数中线性方程组和向量组的最大无关组等概念。通过案例建立数学模型,进一步认识线性方程组的实际应用以及初等行变换求解齐次线性方程组和非齐次线性方程组的方法。
4.2.1 化学方程式的平衡问题
1.问题提出
当丙烷()气体燃烧时,丙烷与氧()结合生成二氧化碳()和水()。该反应的化学方程式如下:
.
配平此方程式。
2.问题分析
为了使反应式平衡,我们必须使反应式两端的碳原子、氢原子及氧原子数目对应相等。根据题目中的假设,此问题实际上是线性方程组的求解问题。
3.模型分析与建立
为了使反应式平衡,必须选择合适的才能使反应式两端的碳原子、氢原子及氧原子数目对应相等。由含三个原子,而含一个原子,故而为维持平衡,必须有
,
类似地,为了平衡原子和原子,必须有
,
.
如果将所有未知量移至等号左边,那么将得到一个齐次线性方程组
(4.4)
方程组是含有3个方程、4个未知数的线性方程组。根据线性代数知识,显然此方程组有非零解。但是为了使化学反应式两端平衡,我们必须找到使得均为最小正整数的解。下面利用矩阵的初等行变换求解此线性方程组。
首先对方程组的系数矩阵进行初等行变换,得到如下行最简形矩阵:
.
求矩阵的行最简形的MATLAB程序如下:
clc, clear
a=[3,0,-1,0;0,2,-2,-1;8,0,0,-2];
b=sym(a); %为了精确求解,转化为符号矩阵
c=rref(b) %求矩阵的行最简形
则可得与原方程组同解的方程组:
4.结果分析
根据上面的求解结果,显然此题中的取法不唯一,但配平化学方程式通常取最小可能的正整数,故取,此时,,。从而配平的化学方程式具有如下形式:
.
4.2.2 工资问题
1.问题提出
现有一个木工、一个电工、一个油漆工和一个粉饰工,四人相互同意彼此装修他们自己的房子。在装修之前,他们约定每人工作13天(包括给自己家干活在内),每人的日工资根据一般的市价在210~260元,每人的日工资数应使得每人的总收入与总支出相等。表4.4是他们协商后制定出的工作天数的分配方案,如何计算出他们每人应得的日工资以及每人房子的装修费(只计算工钱,不包括材料费)是多少?
表4.4 工时分配方案
工种
木工
电工
油漆工
粉饰工
在木工家工作天数
4
3
2
3
在电工家工作天数
5
4
2
3
在油漆工家工作天数
2
5
3
3
在粉饰工家工作天数
2
1
6
4
2.问题分析
这是一个收入-支出的闭合模型,为满足“平衡”条件,每人的收支相等,即要求每人在这13天内“总收入=总支出”。
3.模型建立与求解
设木工、电工、油漆工和粉饰工的日工资分别为元。根据每人在这13天内“总收入=总支出”,可建立如下线性方程组:
整理得齐次线性方程组
对方程组的系数矩阵进行初等行变换,得到如下的行最简形矩阵:
,
则可得与原方程组同解的方程组:
即
.
根据题目要求,应使取值在,故可取,此时,,。
从而木工、电工、油漆工和粉饰工的日工资分别为216元,252元,240元和236元。每人房子的装修费用相当于本人13天的工资,因此分别为2808元,3276元,3120元和3068元。
计算的MATLAB程序如下:
clc, clear
a=[-9,3,2,3;5,-9,2,3;2,5,-10,3;2,1,6,-9]; b=sym(a)
c=rref(b)
x=[-c([1:3],4)',1]*59*4
y=13*x
4.2.3 混凝土配料问题
1.问题提出
混凝土由五种主要的原料组成:水泥、水、砂、石和灰,不同的成分影响混凝土的不同特性。例如,水与水泥的比例影响混凝土的最终强度,砂与石的比例影响混凝土的易加工性,灰与水泥的比例影响混凝土的耐久性等。所以不同用途的混凝土需要不同的原料配比。
一个新建成的混凝土生产企业的设备只能生成存储4种类型的混凝土,但是常用的混凝土有5种类型,即超强型A、通用型B、长寿型C、实用型D和普通型E。它们的配方如表4.5所示。
表4.5 常用混凝土配方
原料
超强型A
通用型B
长寿型C
实用型D
普通型E
水泥
20
18
12
16
16
水
10
10
10
10
12
砂
20
25
15
21
19
石
10
5
15
9
9
灰
0
2
8
4
4
请问你能否帮助该企业从上述5种类型中选出4种作为该企业的产品,使得企业能够在只生产4种类型混凝土的情况下,保证5种类型的混凝土都可以从该企业的产品中获得,说明你的理由。
2.问题分析
从数学上,上述5种基本类型的混凝土可抽象为5个向量。于是问题归结为从5个向量中选出4个向量,使得没有被选出的那个向量可以由这4个被选中的向量线性表示,而且表示系数必须全为正数。根据线性代数知识,上面问题的解决需要首先确定5个向量构成的向量组是否线性相关。
3.模型建立与求解
设五种基本类型的混凝土抽象为5个向量,其数值如下:
.
根据判定向量组线性相关性的方法,可以构造矩阵,若矩阵的秩小于5,则向量组线性相关。
下面对矩阵进行初等行变换,得到如下的行最简形矩阵:
.
显然,从而向量组线性相关。
5个向量构成的向量组线性相关,说明向量组中至少有1个向量可以由其余向量线性表示。但是由于此实际问题是混凝土的配比问题,所以要求1个向量由其余4个向量线性表示时的系数必须是正数。
事实上,由矩阵化简的行最简形矩阵可知,企业可以选择超强型A、通用型B、长寿型C和普通型E作为企业的产品,此时没有被选中的实用型D可以用超强型A、通用型B、长寿型C和普通型E如下表示:
.
计算的MATLAB程序如下:
clc, clear
a=[20,18,12,16,16;10,10,10,10,12;20,25,15,21,19
10,5,15,9,9;0,2,8,4,4];
r=rank(a) %求矩阵的秩
b=rref(a) %求行最简形矩阵
4.2.4 调整气象站观测问题
1.问题提出
某地区有12个气象观测站,10年来各观测站的降水量如表4.6。为了节省开支,想要适当减少气象观测站。问题:减少哪些气象观测站可以使所得的降水量的信息量仍然足够大?
表4.6 降水情况统计表 单位:mm
年份
1
2
3
4
5
6
7
8
9
10
11
12
1981
276.2
324.5
158.6
412.5
292.8
258.4
334.1
303.2
282.9
243.2
159.7
331.2
1982
251.6
287.3
349.5
297.6
227.8
453.6
321.5
451
466.2
307.5
421.1
455.1
1983
192.7
436.2
289.9
366.3
466.2
239.1
357.4
219.7
245.7
411.1
357
353.2
1984
246.2
232.4
243.7
372.5
460.4
158.9
298.7
314.5
256.6
327
296.5
423
1985
291.7
311
502.4
254
245.6
324.8
401
266.5
251.3
289.9
255.4
362.1
1986
466.5
158.9
223.5
425.1
251.4
321
315.4
317.4
246.2
277.5
304.2
410.7
1987
258.6
327.4
432.1
403.9
256.6
282.9
389.7
413.2
466.5
199.3
282.1
387.6
1988
453.4
365.5
357.6
258.1
278.8
467.2
355.2
228.5
453.6
315.6
456.3
407.2
1989
158.5
271
410.2
344.2
250
360.7
376.4
179.4
159.2
342.4
331.2
377.7
1990
324.8
406.5
235.7
288.8
192.6
284.9
290.5
343.7
283.4
281.2
243.7
411.1
2.问题分析
上述12个气象观测站观测到的10年来的降水量可抽象为12个列向量。于是问题归结为从12个列向量中选出尽可能少的向量,使得没有被选出的向量可以由被选出的向量线性表示。根据线性代数知识,若能求得这12个向量构成的向量组的最大无关组,则其最大无关组所对应的气象观测站可将其他的气象观测站的气象资料表示出来,因而其他的气象观测站就是可以减少的。
此外,由于这12个向量都是10维向量,故12个向量构成的向量组必定线性相关,从而一定可以减少某些气象观测站,使所得的降水量的信息量仍然足够大。
3.模型建立与求解
设12个气象观测站观测到的10年来的降水量表示为向量。利用线性代数知识,可以求出向量组的一个最大无关组:。且有
计算的MATLAB程序如下:
clc, clear
a=load('gtable4_6.txt');
b=rref(a)
4.结果分析
根据上述求解结果可知减少第11个观测站和第12个观测站可以使得到的降水量的信息仍然足够大。当然也可以减少另外两个观测站,只要代表这两个观测站的列向量可以由代表其他观测站的列向量线性表示即可。
4.2.5 齐王田忌赛马
1.问题提出
战国时期,有一天齐王提出要与田忌赛马,双方约定从各自的上、中、下三个等级的马中各选一匹参赛,每匹马均只能参赛一次,每一次比赛双方各出一匹马,负者要付胜者千金。已经知道,在同等级的马中,田忌的马不如齐王的马,而如果田忌的马比齐王的马高一等级,则田忌的马可取胜。
由于齐王和田忌的出马策略为“上中下”、“上下中”、“中上下”、“中下上”、“下中上”、“下上中”,都共6种。规定在每个等级的比赛中,赢了记1分,输了记-1分,三个等级的总得分为齐王或田忌的赢得值。
记齐王的策略集为,田忌的策略集为,则齐王的赢得矩阵为
,
求齐王和田忌的混合最优策略。
2.模型建立与求解
设齐王和田忌的最优混合策略分别为和。求和归结为求解如下的两个方程组
(4.5)
和
(4.6)
其中,,,。
实际上方程组和都是7个未知数(也是未知数)的非齐次线性方程组,都有无穷多组解,方程组的解为
即对策值,齐王的混合最优策略为
.
类似地,可以求出田忌的混合最优策略为
.
对策值。
因为方程组有无穷多组解,其中的最小范数解为
,,,
即双方都以的概率选取每个纯策略。或者说在六个纯策略中随机地选取一个即为最优策略。总的结局是齐王赢得的期望值是1千金,田忌所失的期望值也是1千金。
求齐王的混合最优策略的MATLAB程序如下:
clc, clear
a=ones(6); a([1:7:end])=3;
a([4,9,17,24,26,31])=-1 %输入齐王的赢得矩阵
A=[a',-ones(6,1); ones(1,6),0]; %输入7个未知数的线性方程组的系数矩阵
b=[zeros(6,1);1];
xv1=rref([A,b]) %把线性方程组的增广阵化成行最简形
xv=pinv(sym(A))*b %求最小范数解
求田忌的混合最优策略的MATLAB程序如下:
clc, clear
a=ones(6); a([1:7:end])=3;
a([4,9,17,24,26,31])=-1 %输入齐王的赢得矩阵
A=[a,-ones(6,1); ones(1,6),0]; %输入7个未知数的线性方程组的系数矩阵
b=[zeros(6,1);1];
xv1=rref([A,b]) %把线性方程组的增广阵化成行最简形
xv=pinv(sym(A))*b %求最小范数解
3.结果分析
从上面的结果可以看出,在公平的比赛情况下,双方同时提交出马顺序策略,齐王可以有多种可能的策略,齐王都能赢得田忌一千两黄金。之前之所以田忌能赢齐王一千两黄金,其原因在于他事先知道了齐王的出马顺序,而后才做出对自己有利的决策。因此在这类对策问题中,在正式比赛之前,对策双方都应该对自己的策略保密,否则不保密的一方将会处于不利的地位。
4.2.6 服务网点的设置问题
1.问题提出
为适应日益扩大的旅游事业的需要,某城市的甲、乙、丙三个照相馆组成一个联营部,联合经营出租相机的业务。游客可由甲、乙、丙三处任何一处租出相机,用完后,还在三处中任意一处即可。估计其转移概率如表4.7所示。今欲选择其中之一附设相机维修点,问该点设在哪一个照相馆为最好?
表4.7 转移概率值
还 相 机 处
甲
乙
丙
租相机处
甲
乙
丙
0.2
0.8
0.1
0.8
0
0.3
0
0.2
0.6
2.模型建立与求解
由于旅客还相机的情况只与该次租机地点有关,而与相机以前所在的店址无关,所以可用表示相机第次被租时所在的店址;“”、“”、“”分别表示相机第次被租用时在甲、乙、丙馆。则是一个马尔可夫链,由表4.7知其状态转移概率矩阵
.
考虑维修点的设置地点问题,实际上要计算这一马尔可夫链的极限概率分布。
可以证明极限概率分布存在,求归结为解线性方程组
即解方程组
得极限概率,,。
由计算看出,经过长期经营后,该联营部的每架照相机还到甲、乙、丙照相馆的概率分别为、、。由于还到甲馆的照相机较多,因此维修点设在甲馆较好。但由于还到乙馆的相机与还到甲馆的相差不多,若是乙的其它因素更为有利的话,比如,交通较甲方便,便于零配件的运输,电力供应稳定等等,亦可考虑设在乙馆。
计算的MATLAB程序如下:
clc, clear
p=[0.2,0.8,0;0.8,0,0.2;0.1,0.3,0.6]
A=[p'-eye(3);ones(1,3)]; b=[zeros(3,1);1];
x=sym(A)\b
4.3 特征值与特征向量
本节案例涉及线性代数中特征值与特征向量概念,通过这些案例的学习,加深对这些概念的理解,并能应用这些知识解决有关实际问题。
4.3.1 污染与工业发展关系问题
1.问题提出
发展与环境问题已成为21世纪各国政府关注的重点,工业发展势必会引起污染,污染程度与工业发展水平的关系如何,请建立工业增长水平与污染水平之间的关系模型。
2.问题分析
首先要考虑用何指标来度量工业增长水平与污染水平。根据常识,污染水平以空气污染或河湖污染为考察对象,以某种污染指数为测量单位,工业发展水平以某种工业指数为测量单位,为使问题简单化并能说明问题,我们可以假设若干年后污染水平与工业发展水平和目前污染水平与工业发展水平间是一种线性关系。
3.模型建立与求解
设是某地区目前的污染水平(以空气或河湖的某种污染指数为测量单位),是目前的工业发展水平(以某种工业指数为测量单位),把这一年作为起点(亦称基年),记作,若干年后(例如5年后)的污染水平和工业发展水平分别为和,它们之间的关系是
,
.
其中常数为专家经验或某种常识而定。写成矩阵形式,就是
,或,
其中,。
一般地,如果以若干年(如5年)作为一个期间,第个期间的污染和工业发展水平记作和,则有
. (4.7)
记,则式的矩阵形式为
. (4.8)
如果已知该地区基年的水平,利用式就可预测第期时该地区的污染程度和工业发展水平,实际上,由式可得
.
如果直接计算的各次幂,计算将十分繁琐,如果利用矩阵特征值和特征向量的有关性质,不但使计算大大简化,而且模型的结构和性质也更为清晰。为此先计算的特征值和特征向量。
的特征多项式为
.
所以,的特征值为,。
对于,解齐次线性方程组,可得的属于的一个特征向量。
对于解齐次线性方程组,可得的属于的一个特征向量,且线性无关。
因而矩阵可以相似对角化,即存在可逆矩阵,使得,,,计算得到
.
易知时,的第一个分量,可以看出环境的污染也将直接威胁人类的生存。
计算的MATLAB程序如下:
clc, clear, syms k
a=[3,1;2,2]; a=sym(a) %转化为符号矩阵
p=charpoly(a) %求特征多项式
t=roots(p) %求特征根
[vec,val]=eig(a) %求特征向量和特征根
Ak=vec*val^k*inv(vec) %求矩阵A的k次幂
alphak=Ak*[1;7] %求差分方程的解
实际上,由差分方程的求解理论知,差分方程的通解
,
由初值条件,得
,
解之,得
,
因而所求差分方程的解为.
计算的MATLAB程序如下:
clc, clear, syms k
a=[3,1;2,2]; a=sym(a) %转化为符号矩阵
p=charpoly(a) %求特征多项式
t=roots(p) %求特征根
[vec,val]=eig(a) %求特征向量和特征根
c=vec\[1;7] %解方程组待定差分方程特解的常数
s=c(1)*t(1)^k*vec(:,1)+c(2)*t(2)^k*vec(:,2) %写出特解
4.3.2 递推数列的计算
差分方程的求解需要利用矩阵的特征值和特征向量,递推数列的计算有时可以转化为差分方程的计算。下面我们以2003年高考天津卷试题为例来说明。
已知数列中,为常数,且,证明对任意,。
证明 因为,所以
.
令,则特征多项式
.
的特征根为,,相应于的特征向量分别记为
.
,
令,则,,,
于是
,
,
从而有
.
计算的MATLAB程序如下:
clc, clear, syms n a0
A=[-2,1;0,3]; b=sym(A); %转换为符号矩阵
c=charpoly(b) %求特征多项式
[p,val]=eig(b) %求特征值和特征向量
pp=inv(p) %求p的逆阵
An=p*val^n*pp %求矩阵A的n次幂
Bn=An*[a0;1]
an=Bn(1) %提出an的表达式
4.3.3 层次分析法案例
1.问题提出
春天来了,张勇、李雨、王刚、赵宇四位大学生相约去寻找那生机勃勃、盎然向上的春天,去呼吸那沁人
展开阅读全文