第七章 线性代数在后续课程中的应用举例
7.1 电路
7.2 信号与系统课程
7.3 数字信号处理
7.4 静力学
7.5 运动学
7.6 测量学
7.7 文献管理
7.8 经济管理
7.1 电路例 7.1 电阻电路的计算如图 7-1的电路 。 已知,?2?,=4?,
=12?,=4?,=12?,=4?,
=2?,设电压源 =10v,求,,。
图 7-1 应用实例 7.1的电路图
1R 2R
3R 4R 5R 6R
7R S
u 3i 4u 7u
解,用回路电流法进行建模。选回路如图，
设每个网孔。
按图可列出各的回路电流分别为,和
。据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为：
写成矩阵形式为,1 2 3 3
3 3 4 5 5
5 5 6 7
01
0
00
a
bs
c
R R R R i
R R R R R i u
R R R R i
ai bi
ci
1 2 3 3
3 3 4 5 5
5 5 6 7
()
( ) 0
( ) 0
a b s
a b c
bc
R R R i R i u
R i R R R i R i
R i R R R i
也可把参数代入，直接列出数字方程简写成 A*I?b* us
其中 I?[ ； ； ]。今 =10，解此矩阵方程，
即可得问题的解。
因此，可列出 MATLAB程序 ea701
18 12 0 1
12 28 12 0,
0 12 18 0
a
bs
c
i
iu
i
ai bi ci
su
A=[18,-12,0; -12,28,-12; 0,-12,18];
b=[1;0;0]; us=10;
U=rref([A,b*us])
程序运行的结果为：
意味着
1.0000 0 0 0.9259
0 1.0 000 0 0.55 56
0 0 1,000 0 0.37 04
U0
0,9 25 9
0,5 55 6
0,3 70 4
a
b
c
i
i
i
I
7.2 信号与系统课程例 7.3 线性时不变系统的零输入响描述 n阶线性时不变 （ LTI） 连续系统的微分方程为
n≥m
已知 y 及 其 各 阶 导 数 的 初 始 值 为 y(0),
y(1)(0),…,y(n-1)(0),求系统的零输入响应 。
,dddddddddd 111
1
21 ubt
ub
t
ubya
t
ya
t
ya
t
ya
mmm
m
nn
n
n
n
解,?建模当 LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为 0），其形式为（设特征根均为单根）
其中 p1,p2,…,pn是特征方程
a1?n+a2?n-1+…+ an?+ an+1 =0的根，它们可用 roots(a)语句求得。各系数 C1,…,
Cn由 y及其各阶导数的初始值来确定。对此有
1212( ) e e e nptp t p t ny t C C C
C1+ C2+…+ Cn = y0 y0 = y(0)
p1C1+ p2C2+…+ pnCn=Dy0 (Dy0表示 y
的导数的初始值 y(1)(0))
…………………………………
写成矩阵形式为
1 1 1 1
1 1 2 2 0D
n n n n
nnp C p C p C y
0
1
0
0
2
1
11
2
1
1
21
D
D
111
y
y
y
C
C
C
ppp
ppp
n
n
n
n
nn
n
即 V·C = Y0
其解为 C =V \ Y0
式中
V为范德蒙矩阵，在 MATLAB的特殊矩阵库中有 vander函数可直接生成。
11 2 0 0 0[,,,] ; [,D,,D ]nnC C C y y y TT0CY
11
2
1
1
21
111
n
n
nn
n
ppp
ppp
V
MATLAB程序 ea703.m
a=input('输入分母系数向量 a=[a1,a2,...]= ');
n=length(a)-1;
Y0=input(‘输入初始条件向量 Y0=[y0,Dy0，
D2y0,...]= ');
p=roots(a);V=rot90(vander(p));c= V\Y0';
dt=input('dt='); tf=input('tf= ')
t=0:dt:tf; y=zeros(1,length(t));
for k=1:n y= y+c(k)*exp(p(k)*t);end
plot(t,y),grid
程序运行结果用这个通用程序来解一个三阶系统，运行此程序并输入
a=[3,5,7,1]; dt=0.2; tf=8;
而 Y0取 [1,0,0];[0,1,0];[0,0,1]
三种情况，用 hold on语句使三次运行生成的图形画在一幅图上，得到图 7-3。
图 7-3 三阶系统的零输入响应
7.3 数字信号处理例 7.5 数字滤波器系统函数 [12]
数字滤波器的网络结构图实际上也是一种信号流图。它的特点在于所有的相加节点都限定为双输入相加器；另外，数字滤波器器件有一个迟延一个节拍的运算，它也是一个线性算子，它的标注符号为 z?1。根据这样的结构图，也可以用类似于例 7.4的方法，
求它的输入输出之间的传递函数，在数字信号处理中称为系统函数。
图 7-5 某数字滤波器结构图图 7-5表示了某个数字滤波器的结构图，
现在要求出它的系统函数，即输出 y与输入
u之比。先在它的三个中间节点上标注信号的 名称,,，以便对每个节点列写方程。
由于迟延算子 z?1不是数，要用符号代替，
所以取 q? z?1，按照图示情况，可以写出：
1x 2x 3x
写成矩阵形式为经过移项后，系统函数 W可以写成：
12
23
31
2
3 1 1
8 4 4
x qx u
x q x u
xx
11
22
33
00 2
3 1 1
00
8 4 4
01 0 0
q
xx
x q x
xx
x u x = Q x - P u
W = x /u = in v ( I - Q ) * P
现在可以列写计算系统函数的 MATLAB程序 ea705，
syms q
Q(1,2)?q; Q(2,3)=3/8*q?1/4; Q(3,1)=1;
Q(3,3)=0;
P=[2;1/4;0]
W=inv(eye(3)?Q)*P
程序运行的结果为
W = [?16/(?8?3*q^2?2*q)?2*q/(?8?3*q^2?2*q)]
[?2*(3*q?2)/(?8?3*q^2?2*q)?2/(?8?3*q^2?2*q)]
[?16/(?8?3*q^2?2*q)?2*q/(?8?3*q^2?2*q)]
以 y?x3作为输出的系统函数，故再键入
pretty(W(3))
整理后得到
1
2 2 2 1
1 6 2 8 8( 3 )
8 3 2 1,5 4 1,5 4
y q q zW
u q q q q z z
7.4 静力学例 7.6 求双杆系统的支撑反力图 7-6 两杆系统的受力图（左）和分立体受力图（右）
两杆系统受力图见图 7-6，设已知,G1=200;
G2=100; L1= 2; L2=sqrt(2);theta1 =30*pi/180;
theta2 = 45*pi/180;
求所示杆系的支撑反力 Na,Nb,Nc
解,◆画出杆 1和杆 2的分离体图，其中 Na,Nb,Nc
都用其 x,y方向的分量 Nax,Nay,Nbx,Nby,
Ncx,Ncy表示，于是可列出方程如下：
对杆件 1，
ΣX=0 Nax + Ncx = 0
ΣY=0 Nay + Ncy - G1 = 0;
ΣM=0 Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-
G1*L1/2*cos(theta1)=0;
对杆件 2
ΣX=0 Nbx - Ncx = 0;
ΣY=0 Nby - Ncy - G2 = 0;
ΣM=0 Ncy*L2*cos(theta2)+
Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0;
这是一组包含六个未知数 Nax,Nay,Nbx,
Nby,Ncx,Ncy的六个线性代数方程,用手工解时要寻找简化的方法,用了 MATLAB工具，
就可直接列出矩阵方程,（其中 X
为列矩阵,用矩阵除法来解。
◆ MATLAB程序 ea706
G1=200; G2=100; L1= 2; L2 = sqrt (2) ;
% 给原始参数赋值
theta1 = 30*pi/180; theta2 = 45*pi/180;
% 将度化为弧度
AX = B
[ N a x ; N a y ; N b x ; N b y ; N c x ; N c y ]
% 设则按此次序，系数矩阵 A,B可写成下式
A=[1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-
sin(theta1),cos(theta1);..,
0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0,
sin(theta2),cos(theta2)];
B=[0;G1;G1/2*cos(theta1);0;G2;-
G2/2*cos(theta2)];
X = A\B
% 用左除求解线性方程组
◆程序运行的结果为：
即解毕。
1 0 0 0 1,0 0 0 0 0
0 1 0 0 0 1,0 0 0 0
0 0 0 0 - 0,5 0 0 0 0,8 6 6 0
0 0 1 0 - 1,0 0 0 0 0
0 0
A
0 95,09 62
2 0 0,0 0 0 0 1 5 4,9 0 3
8 6,6 0 2 5
,,
0
0 1 0 - 1,0 0 0 0 1 0 0,0000
0 0 0 0 0,7 0 7 1 0,7 0 7 1 - 3 5,3 5 5 3
BX
8
- 9 5,0 9 6 2
1 4 5,0 9 6 2
- 9 5,0 9 6 2
4 5,0 9 6 2
Nax Nay Nbx Nby Ncx Ncy
95.0962 154.9038 -95.0962 145.0962 -95.0962 45.0962
7.5 运动学例 7.7 刚体空间运动学以飞行器为例，它在空中可以围绕三个轴旋转。
假如它在向北飞行，机头正对北方，则它围绕铅垂轴的旋转角称为偏航角（ Yaw），它描述了飞机左右的偏转，用 u表示；围绕翼展轴的旋转角称为倾斜角（ Pitch），它描述了飞机俯仰姿态，用 v表示；
围绕机身轴的旋转角称为滚动角（ Roll），用 w表示；
u,v和 w三个变量统称为欧拉角，它们完全描述了飞机的姿态。
MATLAB中有一个演示程序 quatdemo.m，
专门演示这几个姿态角所造成的飞机状态。
键入：
>>quatdemo
图 7-7 飞行器姿态角演示屏幕上将出现图 7-7的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器。右上方为三个姿态角 u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角。左下方还有 【 静态 】 和 【 动态 】 两个复选钮，我们只介绍 【 静态 】,读者可自行试用 【 动态 】 进行演示。
用键入参数或移动标尺的方法分别给
u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。
用数值来讨论这个程序的实现方法。先把飞行器的三维图像用 N个顶点的三维坐标描述，写成一个 3?N的数据矩阵 G。其顶点次序要适当安排，使得用 plot3命令时按顶点连线能绘制出飞行器的外观。例如以下的程序（ ea707前半部分）即可画出一个最简单的飞行器立体图，如图 7-8所示。
图 7-8 用数据集 G画出的飞行器
Gw=[?4,?3,0;4,?3,0;0,7,0;?4,?3,0]';
Gt=[0,?3,0;0,?3,3;0,2,0;0,?3,0]';
G=[Gw,Gt]
plot3(Gw(1,:),Gw(2,:),Gw(3,:),'r'),hold on
plot3(Gt(1,:),Gt(2,:),Gt(3,:),'g'),
axis equal
运行此程序得出整个飞行器外形的数据集为
（ 7.1）4 4 0 4 0 0 0 0
3 3 7 3 3 3 2 3
0 0 0 0 0 3 0 0
G
飞行器围绕各个轴的旋转的结果，表现为各个顶点坐标发生变化，也就是 G的变化。
只要把三种姿态的变换矩阵 Y,P和 R乘以图形数据矩阵 G即可。其中
（ 7.2）
（ 7.3）
（ 7.4）
c os( ) sin ( ) 0
sin( ) c os( ) 0
0 0 1
uu
uu
Y
1 0 0
0 c os( w ) sin( w )
0 sin ( w ) c os( w )
R
c os( ) 0 sin( )
0 1 0
sin ( ) 0 c os( )
vv
vv
P
单独变化某个姿态角所生成的图形由
G1?Y*G,G2?P*G,G3?R*G算出，如果同时变化三个姿态角，则最后的图像数据成为 Gf?Y*P*R*G?Q*G。这里假定转动的次序为：先滚动 R，再倾斜 P，最后偏航 Y，
由于矩阵乘法不服从交换律，转动次序不同时结果也不同。最后变换矩阵成为
（ 7.5）
c os( ) * c os( ) sin ( ) * c os( w ) c os( ) * sin( ) * sin( w ) sin( ) * sin( w ) c os( ) * sin( ) * c os( w )
sin( ) * c os( ) c os( ) * c os( w ) + sin( ) * sin( ) * sin ( w ) c os( ) * sin( w ) + sin( ) * sin( ) * c os( w )
sin( )
u v u u v u u v
Q u v u u v u u v
v c os( ) * sin( w ),c os( ) * c os( w )
vv
此程序 ea707后半部分如下,
syms u w v
Y=[cos(u),sin(u),0;?sin(u cos(u),0;0,0,1)]
R=[1,0,0;0,cos(w),?sin(w);0,sin(w),cos(w)]
P=[cos(v),0,?sin(v);0,1,0;sin(v),0,cos(v)]
Q=Y*P*R
这里采用了符号运算工具箱以得到普遍的公式，当设定了 u,v,w的具体数值后，比如 u?π/4,v?0,w?π/3，要求出 Q矩阵的数字结果时，要用 subs（代换）命令如下：
A?subs(Q,{u,v,w},{pi/4,0,pi/3})
运行结果为：
0,70 71 0,3 53 6 0,61 24
0,70 71 0,3 53 6 0,61 24
0 0.86 60 0.50 00
A
我们知道,二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列式代表的是三个向量组成的平行六面体的体积。不难算出,det(Y)?1,det(R)?1，
det(P)?1,det(Q)?1。当然也有 det(A)?1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于 1。如果特征值是复数，则它们的绝对值必须为 1。至于对特征向量的要求，读者可从二维情况中的 A5得到启发和推论。因为求特征值的 MATLAB函数 eig还不能用于符号函数，
检验上述结论只能用数值矩阵 A。
键入 [p,lamda]?eig(A) 得到
p0.7701 0.1833? 0.4122i 0.1833?0.4122i
0.3190 0.6702 0.6702
0.5525?0.1315? 0.5745i?0.1315?0.5745i
lamda?1.0000 0 0
0 0.2803? 0.9599i 0
0 0 0.2803? 0.9599i
再键入 abs(lamda)，得到
ans?1.0000,0,0
0 1.0000,0
0,0 1.0000
说明所有的特征值绝对值均为 1。
为了计算各特征向量的范数，可以键入 p'*p，得到
ans?
1.0000 0.0000? 0.0000i 0.0000?0.0000i
0.0000? 0.0000i 1.0000?0.0000?0.0000i
0.0000? 0.0000i?0.0000? 0.0000i 1.0000
可见其特征向量的范数也均为 1。
飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和 3?3的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的 N个顶点。这些顶点用在三维空间中的三个坐标来描述，也就是 3?N
数据集 G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和 4?4变换矩阵。即数据集 G4成为 4?N
矩阵，其最下面一行的元素都取值为 1。
(7.6)
而平移变换矩阵 M和旋转变换矩阵 Y,R,P成为
(7.7)
(7.8)
其中,,为在三个轴,,方向上的平移距离。这种方法在机器人运动学研究中也适用。
3
4
GG=
o ne s( 1,N )
1
2
3
1 0 0
0 1 0
,
0 0 1
0 0 0 1
c
c
c
4M
0
0
0
0 0 0 1
3 3 3
4 4 4
Y,( R,P )
Y,( R,P )
1x 2x 3x1
c 2c 3c
7.6 测量学例 7.8 用坐标测量仪测定圆形工件精密三坐标测量仪可以测量物体表面上任何一点的三维坐标,人们根据这些点的坐标就能推算出物体的其他特征尺寸 。 比如为了测量一个圆形截面的半径,要在 x-y平面内测量其圆周上 n个点的坐标 (,)(i=1,… n)，
然后拟合出其最小二乘圆的半径 。ix iy
设圆周方程为：
其中为圆心的坐标，为半径。对整理上述方程，得到
（ 7.9）
其中,因而,求出就可求出 r。
2 2 212( ) ( )x c y c r
2 2 2 2 21 2 1 2 1 2 32 2 ( ) 2 2x c y c r c c x c y c c x y
2 2 23 1 2c r c c 22
3 1 2r c c c
用 n个测量点坐标 (,)代入，得到
（ 7.10）
这是关于三个未知数的 n个线性方程，所以是一个超定问题。解出就可得知这个最小二乘圆的圆心坐标和半径 r的值：
举一个数字例，设测量了圆周上 7个点，其 x,y
坐标如下：
x = -3.000 -2.000 -1.000 0 1.000 2.000
3.000
y= 3.03 3.90 4.35 4.50 4.40 4.02
3.26
22
1111
1 22
22 22
2
3 22
2 2 1
2 2 1
2 2 1nn
nn
xyxy
c
xy xy
c
c
xy xy
A c = B
iyix
试求此工件的直径、及圆心坐标解,7点应该有 7个行向的方程，其结构相同，
只是数据不同。可以把数据写成列向量，
然后用元素群运算一次列出所有的 7个方程。
其程序为：
x=[-3:3]';
y=[3.03,3.90,4.35,4.50,4.40,4.02,3.26]';
A=[2*x,2*y,ones(size(x))]
B=x.^2+y.^2
c=inv(A'*A)*A'*B,
r=sqrt(c(3)+c(1)^2+c(2)^2)
程序运行的最后结果为：
c = 0.1018 工件圆心的 x坐标
0.4996 工件圆心的 y坐标
15.7533
r = 4.0017 工件的半径 r
为了使更加清楚，把运行此程序前四行得出 A和 B的结果写成方程：
可以看出把 MATLAB的元素群运算和矩阵运算相结合，可以把复杂的计算变得何等简明。
- 6,00 00 6,06 00 1,00 00
- 4,00 00 7,80 00 1,00 00
- 2,00 00 8,70 00 1,00 00
0 9.0000 1.0000
2.0000 8.8000 1.0000
4.0000 8.0400 1.0000
6.0000 6.5200 1
Ac = B
1
2
3
18.1809
19.2100
19.9225
20.2500
20.3600
20.1604
.0 00 0 1 9,62 76
c
c
c
7.7 文献管理例 7.9 情报检索模型假如数据库中包括了 n个文件,而搜索所用的关键词有 m个 。 如果关键词按字母顺序排列,我们就可以把数据库表示为 m?n的矩阵 A。 每个文件用矩阵的列表示 。 A的第 j列的第一个元素是一个数,
它表示第一个关键词出现的相对频率,第二个元素表示第二个关键词出现的相对频率,… 依此类推 。
用于搜索的关键词清单用 空间的向量 x表示 。 如果关键词清单中第 i个关键词在搜索行中出现,则 x的第 i个元素就赋值 1,否则就赋值 0。 为了进行搜索,
只要把 乘以 x。TA
mR
举例来说，假如我们的数据库包含有以下的书名：
B1。应用线性代数
B2，初等线性代数
B3。初等线性代数及其应用
B4，线性代数及其应用
B5。线性代数及应用
B6，矩阵代数及应用
B7。矩阵理论而搜索的 6个关键词组成的集由以下的拼音字母次序排列：
初等，代数，矩阵，理论，线性，应用因为这些关键词在书名中出现最多只有 1次，所以其相对频率数不是 0就是 1。当第 i个关键词出现在第 j本书名上时，元素 A（ i,j）就等于 1，否则就等于 0。这样我们的数据库矩阵就可用下表表示：
关键词书
B1 B2 B3 B4 B5 B6 B7
初等 0 1 1 0 0 0 0
代数 1 1 1 1 1 1 0
矩阵 0 0 0 0 0 1 1
理论 0 0 0 0 0 0 1
线性 1 1 1 1 1 0 0
应用 1 0 1 1 1 1 0
假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为：
搜索结果可以表示为两者的乘积,
于是
0 1 1 0 0 0 0
1 1 1 1 1 1 0
0 0 0 0 0 1 1
0 0 0 0 0 0 1
1 1 1 1 1
A
0
1
0
,
0
0 0 1
1 0 1 1 1 1 0 1
x
0 1 0 0 1 1
1 1 0 0 1 0
1 1 0 0 1 1
0 1 0 0 1 1
0 1 0 0 1 1
0 1
T
y = A x
3
0
2
1
3
0
3
0
3
1
1 0 0 1 2
1
0 0 1 1 0 0 0
Ty A x?
y的各个分量就表示各书与搜索向量匹配的程度。因为3，说明四本书,,,必然包含所有三个关键词。
这四本书就被认为具有最高的匹配度，因而在搜索的结果中把这几本书排在最前面。
本例把线性变换的概念进一步扩展。它不一定是在具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间。在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换。
1y 3y 4y 5y
1B 3B 4B 5B
7.8 经济管理例 7.10 宏观经济模型假定一个国家的经济可以分解为 n个部门,这些部门都有生产产品或服务的独立功能 。 设单列 n元向量 x是这些 n个部门的产出,
组成在 空间的产出向量 。 考虑它有外部的需求,要向外提供产品 。 设 d为外部需求向量,
它也是 空间的单列向量,表示这 n个部门以外的非经济部门对该国经济各部门的需求,
比如政府消费,出口和战略储备等等 。
nR
nR
在各经济部门进行满足外部需求的生产时，
它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂。 Leontiff提出的问题是，
为了满足外部的最终需求向量 d，各生产部门的实际产出 x应该是多少，这对于经济计划的制订当然很有价值。因为
x={内部需求 }?{外部需求 d}
Leontiff的输入输出模型中的一个基本假定是：
对于每个部门，存在着一个在 空间单位消耗列向量,它表示第 i个部门每产出一个单位（比如
100万美金）产品，需要消耗其他部门产出的数量。把这 n个,并列起来，它可以构成一个 n?n
的系数矩阵，可称为内部需求矩阵 V。由于要向外部提供产品，故内部需求矩阵各列向量元素的和必然小于 1（上例中要求它等于 1）。
nR
iv
iv
举一个最简单的例子，假如国民经济由三个部门组成，它们是制造业、农业和服务业。它们的单位消耗列向量如下表。
向下列部门购买每单位输出的输入消耗制造业 农业 服务业制造业 0.5 0.4 0.2
农业 0.2 0.35 0.15
服务业 0.15 0.1 0.3
如果制造业产出了 100个单位的产品，
有 50个单位会被自己消耗,20个单位被农业消耗，而被服务业消耗的是 15个单位，
用算式表示为这就是内部消耗的计算方法，把几个部门都算上，可以写出
内部需求?
0,5 50
10 0 10 0 0,2 20
0,15 15
1v
1
1 1 2 1 3 3 1 2 3 2
3
[,,]
x
x v x v x v v v v x
x
Vx＝
其中于是总的需求方程可以写成为：
x – Vx = d
(I – V )x=d
从而可用 MATLAB语句写出其解的表示式
x = inv( I –V )*d
用一个数字来试算一下，设外部需求为
0.5 0.4 0.2
0.2 0.35 0.15
0.15 0.1 0.3
V
30
20
10
d
用以下程序 ea710解出 x
V=[0.5,0.4,0.2;0.2,0.35,0.15;0.15,0.1,0.3]
d = [30;20;10]
x = inv( eye(3) - V ) * d
程序运行的结果为这个结果是合理的，因为实际产出应该比外部需求大得多，以应付内部的消耗。我们常常说，
某某外部需求可以拉动国民经济增长多少个百分点，就是从这样的模型中得出的。
1 6 0,4 5 6 3
9 4,4 8 6 7
6 2,1 6 7 3
x
内部需求矩阵 V要满足一些基本要求，
一般各列的列向元素总和必须小于 1，否则这个部门就将入不敷出而亏损，但我们仍可能求出上述方程的解。当所有的列向量都出现列向元素总和大于 1的情况时，解 x
中会出现负值，因而是庸解。 请注意分析其实际意义。
7.1 电路
7.2 信号与系统课程
7.3 数字信号处理
7.4 静力学
7.5 运动学
7.6 测量学
7.7 文献管理
7.8 经济管理
7.1 电路例 7.1 电阻电路的计算如图 7-1的电路 。 已知,?2?,=4?,
=12?,=4?,=12?,=4?,
=2?,设电压源 =10v,求,,。
图 7-1 应用实例 7.1的电路图
1R 2R
3R 4R 5R 6R
7R S
u 3i 4u 7u
解,用回路电流法进行建模。选回路如图，
设每个网孔。
按图可列出各的回路电流分别为,和
。据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为：
写成矩阵形式为,1 2 3 3
3 3 4 5 5
5 5 6 7
01
0
00
a
bs
c
R R R R i
R R R R R i u
R R R R i
ai bi
ci
1 2 3 3
3 3 4 5 5
5 5 6 7
()
( ) 0
( ) 0
a b s
a b c
bc
R R R i R i u
R i R R R i R i
R i R R R i
也可把参数代入，直接列出数字方程简写成 A*I?b* us
其中 I?[ ； ； ]。今 =10，解此矩阵方程，
即可得问题的解。
因此，可列出 MATLAB程序 ea701
18 12 0 1
12 28 12 0,
0 12 18 0
a
bs
c
i
iu
i
ai bi ci
su
A=[18,-12,0; -12,28,-12; 0,-12,18];
b=[1;0;0]; us=10;
U=rref([A,b*us])
程序运行的结果为：
意味着
1.0000 0 0 0.9259
0 1.0 000 0 0.55 56
0 0 1,000 0 0.37 04
U0
0,9 25 9
0,5 55 6
0,3 70 4
a
b
c
i
i
i
I
7.2 信号与系统课程例 7.3 线性时不变系统的零输入响描述 n阶线性时不变 （ LTI） 连续系统的微分方程为
n≥m
已知 y 及 其 各 阶 导 数 的 初 始 值 为 y(0),
y(1)(0),…,y(n-1)(0),求系统的零输入响应 。
,dddddddddd 111
1
21 ubt
ub
t
ubya
t
ya
t
ya
t
ya
mmm
m
nn
n
n
n
解,?建模当 LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为 0），其形式为（设特征根均为单根）
其中 p1,p2,…,pn是特征方程
a1?n+a2?n-1+…+ an?+ an+1 =0的根，它们可用 roots(a)语句求得。各系数 C1,…,
Cn由 y及其各阶导数的初始值来确定。对此有
1212( ) e e e nptp t p t ny t C C C
C1+ C2+…+ Cn = y0 y0 = y(0)
p1C1+ p2C2+…+ pnCn=Dy0 (Dy0表示 y
的导数的初始值 y(1)(0))
…………………………………
写成矩阵形式为
1 1 1 1
1 1 2 2 0D
n n n n
nnp C p C p C y
0
1
0
0
2
1
11
2
1
1
21
D
D
111
y
y
y
C
C
C
ppp
ppp
n
n
n
n
nn
n
即 V·C = Y0
其解为 C =V \ Y0
式中
V为范德蒙矩阵，在 MATLAB的特殊矩阵库中有 vander函数可直接生成。
11 2 0 0 0[,,,] ; [,D,,D ]nnC C C y y y TT0CY
11
2
1
1
21
111
n
n
nn
n
ppp
ppp
V
MATLAB程序 ea703.m
a=input('输入分母系数向量 a=[a1,a2,...]= ');
n=length(a)-1;
Y0=input(‘输入初始条件向量 Y0=[y0,Dy0，
D2y0,...]= ');
p=roots(a);V=rot90(vander(p));c= V\Y0';
dt=input('dt='); tf=input('tf= ')
t=0:dt:tf; y=zeros(1,length(t));
for k=1:n y= y+c(k)*exp(p(k)*t);end
plot(t,y),grid
程序运行结果用这个通用程序来解一个三阶系统，运行此程序并输入
a=[3,5,7,1]; dt=0.2; tf=8;
而 Y0取 [1,0,0];[0,1,0];[0,0,1]
三种情况，用 hold on语句使三次运行生成的图形画在一幅图上，得到图 7-3。
图 7-3 三阶系统的零输入响应
7.3 数字信号处理例 7.5 数字滤波器系统函数 [12]
数字滤波器的网络结构图实际上也是一种信号流图。它的特点在于所有的相加节点都限定为双输入相加器；另外，数字滤波器器件有一个迟延一个节拍的运算，它也是一个线性算子，它的标注符号为 z?1。根据这样的结构图，也可以用类似于例 7.4的方法，
求它的输入输出之间的传递函数，在数字信号处理中称为系统函数。
图 7-5 某数字滤波器结构图图 7-5表示了某个数字滤波器的结构图，
现在要求出它的系统函数，即输出 y与输入
u之比。先在它的三个中间节点上标注信号的 名称,,，以便对每个节点列写方程。
由于迟延算子 z?1不是数，要用符号代替，
所以取 q? z?1，按照图示情况，可以写出：
1x 2x 3x
写成矩阵形式为经过移项后，系统函数 W可以写成：
12
23
31
2
3 1 1
8 4 4
x qx u
x q x u
xx
11
22
33
00 2
3 1 1
00
8 4 4
01 0 0
q
xx
x q x
xx
x u x = Q x - P u
W = x /u = in v ( I - Q ) * P
现在可以列写计算系统函数的 MATLAB程序 ea705，
syms q
Q(1,2)?q; Q(2,3)=3/8*q?1/4; Q(3,1)=1;
Q(3,3)=0;
P=[2;1/4;0]
W=inv(eye(3)?Q)*P
程序运行的结果为
W = [?16/(?8?3*q^2?2*q)?2*q/(?8?3*q^2?2*q)]
[?2*(3*q?2)/(?8?3*q^2?2*q)?2/(?8?3*q^2?2*q)]
[?16/(?8?3*q^2?2*q)?2*q/(?8?3*q^2?2*q)]
以 y?x3作为输出的系统函数，故再键入
pretty(W(3))
整理后得到
1
2 2 2 1
1 6 2 8 8( 3 )
8 3 2 1,5 4 1,5 4
y q q zW
u q q q q z z
7.4 静力学例 7.6 求双杆系统的支撑反力图 7-6 两杆系统的受力图（左）和分立体受力图（右）
两杆系统受力图见图 7-6，设已知,G1=200;
G2=100; L1= 2; L2=sqrt(2);theta1 =30*pi/180;
theta2 = 45*pi/180;
求所示杆系的支撑反力 Na,Nb,Nc
解,◆画出杆 1和杆 2的分离体图，其中 Na,Nb,Nc
都用其 x,y方向的分量 Nax,Nay,Nbx,Nby,
Ncx,Ncy表示，于是可列出方程如下：
对杆件 1，
ΣX=0 Nax + Ncx = 0
ΣY=0 Nay + Ncy - G1 = 0;
ΣM=0 Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-
G1*L1/2*cos(theta1)=0;
对杆件 2
ΣX=0 Nbx - Ncx = 0;
ΣY=0 Nby - Ncy - G2 = 0;
ΣM=0 Ncy*L2*cos(theta2)+
Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0;
这是一组包含六个未知数 Nax,Nay,Nbx,
Nby,Ncx,Ncy的六个线性代数方程,用手工解时要寻找简化的方法,用了 MATLAB工具，
就可直接列出矩阵方程,（其中 X
为列矩阵,用矩阵除法来解。
◆ MATLAB程序 ea706
G1=200; G2=100; L1= 2; L2 = sqrt (2) ;
% 给原始参数赋值
theta1 = 30*pi/180; theta2 = 45*pi/180;
% 将度化为弧度
AX = B
[ N a x ; N a y ; N b x ; N b y ; N c x ; N c y ]
% 设则按此次序，系数矩阵 A,B可写成下式
A=[1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-
sin(theta1),cos(theta1);..,
0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0,
sin(theta2),cos(theta2)];
B=[0;G1;G1/2*cos(theta1);0;G2;-
G2/2*cos(theta2)];
X = A\B
% 用左除求解线性方程组
◆程序运行的结果为：
即解毕。
1 0 0 0 1,0 0 0 0 0
0 1 0 0 0 1,0 0 0 0
0 0 0 0 - 0,5 0 0 0 0,8 6 6 0
0 0 1 0 - 1,0 0 0 0 0
0 0
A
0 95,09 62
2 0 0,0 0 0 0 1 5 4,9 0 3
8 6,6 0 2 5
,,
0
0 1 0 - 1,0 0 0 0 1 0 0,0000
0 0 0 0 0,7 0 7 1 0,7 0 7 1 - 3 5,3 5 5 3
BX
8
- 9 5,0 9 6 2
1 4 5,0 9 6 2
- 9 5,0 9 6 2
4 5,0 9 6 2
Nax Nay Nbx Nby Ncx Ncy
95.0962 154.9038 -95.0962 145.0962 -95.0962 45.0962
7.5 运动学例 7.7 刚体空间运动学以飞行器为例，它在空中可以围绕三个轴旋转。
假如它在向北飞行，机头正对北方，则它围绕铅垂轴的旋转角称为偏航角（ Yaw），它描述了飞机左右的偏转，用 u表示；围绕翼展轴的旋转角称为倾斜角（ Pitch），它描述了飞机俯仰姿态，用 v表示；
围绕机身轴的旋转角称为滚动角（ Roll），用 w表示；
u,v和 w三个变量统称为欧拉角，它们完全描述了飞机的姿态。
MATLAB中有一个演示程序 quatdemo.m，
专门演示这几个姿态角所造成的飞机状态。
键入：
>>quatdemo
图 7-7 飞行器姿态角演示屏幕上将出现图 7-7的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器。右上方为三个姿态角 u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角。左下方还有 【 静态 】 和 【 动态 】 两个复选钮，我们只介绍 【 静态 】,读者可自行试用 【 动态 】 进行演示。
用键入参数或移动标尺的方法分别给
u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。
用数值来讨论这个程序的实现方法。先把飞行器的三维图像用 N个顶点的三维坐标描述，写成一个 3?N的数据矩阵 G。其顶点次序要适当安排，使得用 plot3命令时按顶点连线能绘制出飞行器的外观。例如以下的程序（ ea707前半部分）即可画出一个最简单的飞行器立体图，如图 7-8所示。
图 7-8 用数据集 G画出的飞行器
Gw=[?4,?3,0;4,?3,0;0,7,0;?4,?3,0]';
Gt=[0,?3,0;0,?3,3;0,2,0;0,?3,0]';
G=[Gw,Gt]
plot3(Gw(1,:),Gw(2,:),Gw(3,:),'r'),hold on
plot3(Gt(1,:),Gt(2,:),Gt(3,:),'g'),
axis equal
运行此程序得出整个飞行器外形的数据集为
（ 7.1）4 4 0 4 0 0 0 0
3 3 7 3 3 3 2 3
0 0 0 0 0 3 0 0
G
飞行器围绕各个轴的旋转的结果，表现为各个顶点坐标发生变化，也就是 G的变化。
只要把三种姿态的变换矩阵 Y,P和 R乘以图形数据矩阵 G即可。其中
（ 7.2）
（ 7.3）
（ 7.4）
c os( ) sin ( ) 0
sin( ) c os( ) 0
0 0 1
uu
uu
Y
1 0 0
0 c os( w ) sin( w )
0 sin ( w ) c os( w )
R
c os( ) 0 sin( )
0 1 0
sin ( ) 0 c os( )
vv
vv
P
单独变化某个姿态角所生成的图形由
G1?Y*G,G2?P*G,G3?R*G算出，如果同时变化三个姿态角，则最后的图像数据成为 Gf?Y*P*R*G?Q*G。这里假定转动的次序为：先滚动 R，再倾斜 P，最后偏航 Y，
由于矩阵乘法不服从交换律，转动次序不同时结果也不同。最后变换矩阵成为
（ 7.5）
c os( ) * c os( ) sin ( ) * c os( w ) c os( ) * sin( ) * sin( w ) sin( ) * sin( w ) c os( ) * sin( ) * c os( w )
sin( ) * c os( ) c os( ) * c os( w ) + sin( ) * sin( ) * sin ( w ) c os( ) * sin( w ) + sin( ) * sin( ) * c os( w )
sin( )
u v u u v u u v
Q u v u u v u u v
v c os( ) * sin( w ),c os( ) * c os( w )
vv
此程序 ea707后半部分如下,
syms u w v
Y=[cos(u),sin(u),0;?sin(u cos(u),0;0,0,1)]
R=[1,0,0;0,cos(w),?sin(w);0,sin(w),cos(w)]
P=[cos(v),0,?sin(v);0,1,0;sin(v),0,cos(v)]
Q=Y*P*R
这里采用了符号运算工具箱以得到普遍的公式，当设定了 u,v,w的具体数值后，比如 u?π/4,v?0,w?π/3，要求出 Q矩阵的数字结果时，要用 subs（代换）命令如下：
A?subs(Q,{u,v,w},{pi/4,0,pi/3})
运行结果为：
0,70 71 0,3 53 6 0,61 24
0,70 71 0,3 53 6 0,61 24
0 0.86 60 0.50 00
A
我们知道,二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列式代表的是三个向量组成的平行六面体的体积。不难算出,det(Y)?1,det(R)?1，
det(P)?1,det(Q)?1。当然也有 det(A)?1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于 1。如果特征值是复数，则它们的绝对值必须为 1。至于对特征向量的要求，读者可从二维情况中的 A5得到启发和推论。因为求特征值的 MATLAB函数 eig还不能用于符号函数，
检验上述结论只能用数值矩阵 A。
键入 [p,lamda]?eig(A) 得到
p0.7701 0.1833? 0.4122i 0.1833?0.4122i
0.3190 0.6702 0.6702
0.5525?0.1315? 0.5745i?0.1315?0.5745i
lamda?1.0000 0 0
0 0.2803? 0.9599i 0
0 0 0.2803? 0.9599i
再键入 abs(lamda)，得到
ans?1.0000,0,0
0 1.0000,0
0,0 1.0000
说明所有的特征值绝对值均为 1。
为了计算各特征向量的范数，可以键入 p'*p，得到
ans?
1.0000 0.0000? 0.0000i 0.0000?0.0000i
0.0000? 0.0000i 1.0000?0.0000?0.0000i
0.0000? 0.0000i?0.0000? 0.0000i 1.0000
可见其特征向量的范数也均为 1。
飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和 3?3的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的 N个顶点。这些顶点用在三维空间中的三个坐标来描述，也就是 3?N
数据集 G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和 4?4变换矩阵。即数据集 G4成为 4?N
矩阵，其最下面一行的元素都取值为 1。
(7.6)
而平移变换矩阵 M和旋转变换矩阵 Y,R,P成为
(7.7)
(7.8)
其中,,为在三个轴,,方向上的平移距离。这种方法在机器人运动学研究中也适用。
3
4
GG=
o ne s( 1,N )
1
2
3
1 0 0
0 1 0
,
0 0 1
0 0 0 1
c
c
c
4M
0
0
0
0 0 0 1
3 3 3
4 4 4
Y,( R,P )
Y,( R,P )
1x 2x 3x1
c 2c 3c
7.6 测量学例 7.8 用坐标测量仪测定圆形工件精密三坐标测量仪可以测量物体表面上任何一点的三维坐标,人们根据这些点的坐标就能推算出物体的其他特征尺寸 。 比如为了测量一个圆形截面的半径,要在 x-y平面内测量其圆周上 n个点的坐标 (,)(i=1,… n)，
然后拟合出其最小二乘圆的半径 。ix iy
设圆周方程为：
其中为圆心的坐标，为半径。对整理上述方程，得到
（ 7.9）
其中,因而,求出就可求出 r。
2 2 212( ) ( )x c y c r
2 2 2 2 21 2 1 2 1 2 32 2 ( ) 2 2x c y c r c c x c y c c x y
2 2 23 1 2c r c c 22
3 1 2r c c c
用 n个测量点坐标 (,)代入，得到
（ 7.10）
这是关于三个未知数的 n个线性方程，所以是一个超定问题。解出就可得知这个最小二乘圆的圆心坐标和半径 r的值：
举一个数字例，设测量了圆周上 7个点，其 x,y
坐标如下：
x = -3.000 -2.000 -1.000 0 1.000 2.000
3.000
y= 3.03 3.90 4.35 4.50 4.40 4.02
3.26
22
1111
1 22
22 22
2
3 22
2 2 1
2 2 1
2 2 1nn
nn
xyxy
c
xy xy
c
c
xy xy
A c = B
iyix
试求此工件的直径、及圆心坐标解,7点应该有 7个行向的方程，其结构相同，
只是数据不同。可以把数据写成列向量，
然后用元素群运算一次列出所有的 7个方程。
其程序为：
x=[-3:3]';
y=[3.03,3.90,4.35,4.50,4.40,4.02,3.26]';
A=[2*x,2*y,ones(size(x))]
B=x.^2+y.^2
c=inv(A'*A)*A'*B,
r=sqrt(c(3)+c(1)^2+c(2)^2)
程序运行的最后结果为：
c = 0.1018 工件圆心的 x坐标
0.4996 工件圆心的 y坐标
15.7533
r = 4.0017 工件的半径 r
为了使更加清楚，把运行此程序前四行得出 A和 B的结果写成方程：
可以看出把 MATLAB的元素群运算和矩阵运算相结合，可以把复杂的计算变得何等简明。
- 6,00 00 6,06 00 1,00 00
- 4,00 00 7,80 00 1,00 00
- 2,00 00 8,70 00 1,00 00
0 9.0000 1.0000
2.0000 8.8000 1.0000
4.0000 8.0400 1.0000
6.0000 6.5200 1
Ac = B
1
2
3
18.1809
19.2100
19.9225
20.2500
20.3600
20.1604
.0 00 0 1 9,62 76
c
c
c
7.7 文献管理例 7.9 情报检索模型假如数据库中包括了 n个文件,而搜索所用的关键词有 m个 。 如果关键词按字母顺序排列,我们就可以把数据库表示为 m?n的矩阵 A。 每个文件用矩阵的列表示 。 A的第 j列的第一个元素是一个数,
它表示第一个关键词出现的相对频率,第二个元素表示第二个关键词出现的相对频率,… 依此类推 。
用于搜索的关键词清单用 空间的向量 x表示 。 如果关键词清单中第 i个关键词在搜索行中出现,则 x的第 i个元素就赋值 1,否则就赋值 0。 为了进行搜索,
只要把 乘以 x。TA
mR
举例来说，假如我们的数据库包含有以下的书名：
B1。应用线性代数
B2，初等线性代数
B3。初等线性代数及其应用
B4，线性代数及其应用
B5。线性代数及应用
B6，矩阵代数及应用
B7。矩阵理论而搜索的 6个关键词组成的集由以下的拼音字母次序排列：
初等，代数，矩阵，理论，线性，应用因为这些关键词在书名中出现最多只有 1次，所以其相对频率数不是 0就是 1。当第 i个关键词出现在第 j本书名上时，元素 A（ i,j）就等于 1，否则就等于 0。这样我们的数据库矩阵就可用下表表示：
关键词书
B1 B2 B3 B4 B5 B6 B7
初等 0 1 1 0 0 0 0
代数 1 1 1 1 1 1 0
矩阵 0 0 0 0 0 1 1
理论 0 0 0 0 0 0 1
线性 1 1 1 1 1 0 0
应用 1 0 1 1 1 1 0
假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为：
搜索结果可以表示为两者的乘积,
于是
0 1 1 0 0 0 0
1 1 1 1 1 1 0
0 0 0 0 0 1 1
0 0 0 0 0 0 1
1 1 1 1 1
A
0
1
0
,
0
0 0 1
1 0 1 1 1 1 0 1
x
0 1 0 0 1 1
1 1 0 0 1 0
1 1 0 0 1 1
0 1 0 0 1 1
0 1 0 0 1 1
0 1
T
y = A x
3
0
2
1
3
0
3
0
3
1
1 0 0 1 2
1
0 0 1 1 0 0 0
Ty A x?
y的各个分量就表示各书与搜索向量匹配的程度。因为3，说明四本书,,,必然包含所有三个关键词。
这四本书就被认为具有最高的匹配度，因而在搜索的结果中把这几本书排在最前面。
本例把线性变换的概念进一步扩展。它不一定是在具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间。在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换。
1y 3y 4y 5y
1B 3B 4B 5B
7.8 经济管理例 7.10 宏观经济模型假定一个国家的经济可以分解为 n个部门,这些部门都有生产产品或服务的独立功能 。 设单列 n元向量 x是这些 n个部门的产出,
组成在 空间的产出向量 。 考虑它有外部的需求,要向外提供产品 。 设 d为外部需求向量,
它也是 空间的单列向量,表示这 n个部门以外的非经济部门对该国经济各部门的需求,
比如政府消费,出口和战略储备等等 。
nR
nR
在各经济部门进行满足外部需求的生产时，
它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂。 Leontiff提出的问题是，
为了满足外部的最终需求向量 d，各生产部门的实际产出 x应该是多少，这对于经济计划的制订当然很有价值。因为
x={内部需求 }?{外部需求 d}
Leontiff的输入输出模型中的一个基本假定是：
对于每个部门，存在着一个在 空间单位消耗列向量,它表示第 i个部门每产出一个单位（比如
100万美金）产品，需要消耗其他部门产出的数量。把这 n个,并列起来，它可以构成一个 n?n
的系数矩阵，可称为内部需求矩阵 V。由于要向外部提供产品，故内部需求矩阵各列向量元素的和必然小于 1（上例中要求它等于 1）。
nR
iv
iv
举一个最简单的例子，假如国民经济由三个部门组成，它们是制造业、农业和服务业。它们的单位消耗列向量如下表。
向下列部门购买每单位输出的输入消耗制造业 农业 服务业制造业 0.5 0.4 0.2
农业 0.2 0.35 0.15
服务业 0.15 0.1 0.3
如果制造业产出了 100个单位的产品，
有 50个单位会被自己消耗,20个单位被农业消耗，而被服务业消耗的是 15个单位，
用算式表示为这就是内部消耗的计算方法，把几个部门都算上，可以写出
内部需求?
0,5 50
10 0 10 0 0,2 20
0,15 15
1v
1
1 1 2 1 3 3 1 2 3 2
3
[,,]
x
x v x v x v v v v x
x
Vx＝
其中于是总的需求方程可以写成为：
x – Vx = d
(I – V )x=d
从而可用 MATLAB语句写出其解的表示式
x = inv( I –V )*d
用一个数字来试算一下，设外部需求为
0.5 0.4 0.2
0.2 0.35 0.15
0.15 0.1 0.3
V
30
20
10
d
用以下程序 ea710解出 x
V=[0.5,0.4,0.2;0.2,0.35,0.15;0.15,0.1,0.3]
d = [30;20;10]
x = inv( eye(3) - V ) * d
程序运行的结果为这个结果是合理的，因为实际产出应该比外部需求大得多，以应付内部的消耗。我们常常说，
某某外部需求可以拉动国民经济增长多少个百分点，就是从这样的模型中得出的。
1 6 0,4 5 6 3
9 4,4 8 6 7
6 2,1 6 7 3
x
内部需求矩阵 V要满足一些基本要求，
一般各列的列向元素总和必须小于 1，否则这个部门就将入不敷出而亏损，但我们仍可能求出上述方程的解。当所有的列向量都出现列向元素总和大于 1的情况时，解 x
中会出现负值，因而是庸解。 请注意分析其实际意义。
