线性代数在工程中的应用实例

第7章 在科技及工程中的应用实例 ......................................................................... 1

7.1 由拉压杆组成的桁架结构 ........................................................................................................... 1 7.2 格型梯形滤波器系统函数的推导 ............................................................................................... 1 7.3 计算频谱用的DFT矩阵 ............................................................................................................. 2 7.4 显示器色彩制式转换问题 ........................................................................................................... 4 7.5 人员流动问题 ............................................................................................................................... 5 7.6 二氧化碳分子结构的振动频率 ................................................................................................... 5 7.7 二自由度机械振动 ....................................................................................................................... 6 7.8 FIR数字滤波器最优化设计[12] ................................................................................................... 8 7.9 弹性梁的柔度矩阵 ....................................................................................................................... 9 7.10 用二次样条函数插值5个点 ................................................................................................... 11 7.11 飞行器三维空间运动的矩阵描述 ........................................................................................... 12 7.12 金融公司支付基金的流动 ....................................................................................................... 14 7.13 质谱图实验结果分析 ............................................................................................................. 15 7.14 用特征方程解Fibonacci数列问题 ......................................................................................... 16 7.15 简单线性规划问题 ................................................................................................................... 18

·1· 第2章 时域中的离散信号和系统

第7章 在科技及工程中的应用实例

7.1 由拉压杆组成的桁架结构

由13根拉压杆件组成的桁架结构，如图7-1

所示，13个平衡方程已给出，它们来自6个中间节点，每个节点有x,y两个方向的平衡方程，还有一个整体结构的y方向平衡方程。现求其各杆所受的力。

解：按照题给方程组改写成矩阵形式，令

k1cos114/16^214^20.6585k2cos216/16^216^20.7071列k1sin116/16^214^20.7526F2+k2F1=0

图7-1 由拉压杆件组成的桁架结构

方程时假设各杆的受力均为拉力，其相应的方程组及化为矩阵后的形式为：

k2100000000000F10F-F2+F6=0 0-100010000000200010000000000F32000F3=2000F4+k1F5-k2F1=0-k001k0000000012F40k0 10k00000000F-1000k2F1+F3+k3F5=-1000325F7+k1F8-F4=0 000-1001k100000F60-500k3F8+F9= -5000000000k310000F70000-k1-10001000F80F10-k1F5-F6=0 F9+k3F5= 40000000k300010000F94000000000-1000k00F0k2F11-F7=0, 210k2F11+F12=-5000000000000k210F11-5000000000k00010F2000F12+k3F8= 200031200000000000k01FF13+k2F11=0213

(7.1.1)将它看作A*F=B，编成的程序为pla701，核心语句为给A,B赋值，再求F=A\\B，结果为： F=[ -7236; 5117; 2000; -6969; 2812; 5117; -4883; -3167; 1883; 6969; -6906; 4383; 4883 ]

其中负号表示杆受的是压力。

7.2 格型梯形滤波器系统函数的推导

使用计算机解题后，用矩阵模型几乎是最简便的数学方法了。这将给后续课的建模

图7-2 三阶格型梯形滤波器的结构图

1

·2· 第2章 时域中的离散信号和系统

和计算带来性的好处。例如要求出图7-2所示的滤波器的系统函数：

先列出方程，令q=z-1，得到

x1 u k3x4; x2x1; x3k3x2x4; x4qx7; x5 x2 k2x8; x6x5;

x7k2x6x8; x8qx11; x9 x6 k1x12; x10x9; x11k1x10x12; x12qx10; x13y C0x12 C1x11 C2x7 C3x3

这是一组含有13个变量的13个联立方程，用过去的手工方法一个一个消元，理论上是可行的，但它运算极其繁琐，可以预期，95%以上的师生恐怕一个小时也解不出来，而且做对的概率极低。

用矩阵的思路和方法来解就完全不同，它不是通过消元来减少变量，而是想办法补上所有的零元素，把方程扩充为完整的矩阵形式：

x1u-k3x4 x1 0, 0, 0, -k3, 0, 0, 0, 0, 0, 0, 0, 0, 0x11xx2x1 x 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0220x3k3x2x4x3 0, k3, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0x30x4qx7 x 0, 0, 0, 0, 0, 0, q, 0, 0, 0, 0, 0, 04x40x5x2- k2x8 x5 0, 1, 0, 0, 0, 0, 0, -k2, 0, 0, 0, 0, 0x50x6x5x 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 06x60

 0, 0, 0, 0, 0, k, 0, 1, 0, 0, 0, 0, 0x0ux7k2x6x8x727x 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, q, 0, 0x80x8qx1180x9x6-k1x12x9 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -k1, 0x9 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0x10x9x100x10x11k1x10x12x11 0, 0, 0, 0, 0, 0, 0, 0, 0, k1, 0, 1, 0x1100x12qx10x12 0, 0, 0, 0, 0, 0, 0, 0, 0, q, 0, 0, 0x12x 0, 0, C, 0, 0, 0, C, 0, 0, 0, C, C, 0x0x1332101313

X=QX+PU,X/U=inv(I-Q)*P

看似把模型搞复杂了，其实计算却非常容易。程序pla703先对P，Q矩阵赋值，键入

W=inv(I-Q)*P，马上就得出了系统函数。

编程时要注意，本例虽然是数值计算，但计算的内容中带有z变换算子q=z-1，所以P，Q矩阵仍然必须用符号属性，对P，Q赋值时第一个元素必须取含q的算式。熟练后不必列出Q和P的矩阵形式，可以按其下标规律直接进行元素赋值。

用以下参数：k0=1, k1=1/4, k2=1/2, k3=1/3, C0=-0.2, C1=0.8, C2=1.5, C3=1，编成了程序pla703。运行此程序就得到：

x(13)24q349q242.9q30.824z349z242.9z130.8W(13) 32u8q15q13q248z315z213z124用矩阵模型解信号流图的最大优点是一步到位，依靠计算机，既快速，又极易查错。

7.3 计算频谱用的DFT矩阵

有限长序列x(n)(0≤n≤N-1)有N个样本值。它的傅里叶变换X(k)在频率区间（0≤ω＜2π）的N个等间隔分布的点ωk=2πk/N(0≤k≤N-1)上也有N个样本值。这两组有限长的序列之间可以用简单的关系联系起来：

2

·3· 第2章 时域中的离散信号和系统

N1n X(k)x(n)WknN0kN1

其中WNexp2j/N是一个相角2/N为的单位向量，也称为旋转因子。X(k)就称为x(n)的离散傅里叶变换(DFT)，也就是x(n)的频谱。用矩阵来表示,可写成:

所以信号频谱的计算，可以简单地用一个矩阵乘法来完成。信号是N×1维数组，矩阵W称为DFT矩阵，它是N×N维的复数阵。把矩阵乘法看作一个变换，我们就可以把频谱计算看作信号从时域到频域的线性变换。

这个矩阵运算可用下列几条MATLAB语句来实现，程序名为pla704。

xn=input('xn=(长度为N的数组) '), N=length(xn); % 输入数据 n = [0:1:N-1]; k = [0:1:N-1]; % 设定n和k的行向量 WN = exp(-j*2*pi/N); % 设定Wn因子 nk = n'*k; % 产生一个含nk值的N乘N维的整数矩阵 WNnk = WN .^ nk; % 求出W矩阵 Xk = xn * WNnk % 求出离散傅里叶级数系数 plot(k,Xk) % 画出幅频特性

前两条语句无须解释。第三条语句把n转为列向量n’，再与行向量k做矩阵乘，它产生的是一个整数矩阵：

这个整数方阵恰好是算式中W矩阵的指数部分。第四条语句就产生了W矩阵，而最后一条语句就完成了DFS的全部运算。由于实际上离散傅里叶级数并不太用到。它已被下节将介绍的离散傅里叶变换取代，而且两者的计算程序是完全相同的。所以如果写成MATLAB子程序，它可以命名为DFT。

在这5行DFT子程序前面，加上输入语句：

xn=input('xn= '); N=length(xn);

10.5001002003在子程序后面，加上绘图语句：

subplot(2,1,1),stem(xn,'.'); subplot(2,1,2),stem(abs(Xk),'.')

就得到本题的程序pla604。 100运行pla604，并按提示输入xn，即可在子图150上得到xn序列的波形，并在子图2上得到此信号的幅00100200 3频特性。例如输入序列为

图7-3 序列x(n)及其频谱 xn=[ones(1,),zeros(1,192)]，则得到的时域

波形及其频谱图如图7-3。

计算每一个DFT分量X(i)，需要做N次复数乘法和N-1次复数加法。而要完成整个DFT，求出X，则需要做N×N次复数乘法和N(N-1)次复数加法。另外，它占的数据存储量为N×N。在频谱计算中，N的典型值是1024。所以参与运算的数据达百万，需要的存储量

3

·4· 第2章 时域中的离散信号和系统

超过16M， 乘法运算的次数超过数百万。因此，在低档的微机上，当N较大时，MATLAB

会拒绝执行这个程序，并警告内存不足。实际上，MATLAB内置的是加快计算并节省内存快速傅里叶变换(FFT)程序。

从这里也可以看出，矩阵给我们提供了一个组织海量数据进行运算的最好方法，对上百万数据进行赋值、运算和绘图，只用了几行程序。线性代数的重要性之所以在近20年来可与微积分相妣美，这是一个重要原因。如果学线性代数而不涉及工程计算实践，那就无法真正体会线性代数的精髓所在。

7.4 显示器色彩制式转换问题

彩色显示器使用红(R)、绿(G)和蓝(B)光的叠成效应生成颜色。 显示器屏幕的内表面由微粒象素组成, 每个微粒包括三个荧光点: 红、绿、蓝。 电子位于屏幕的后方, 向屏幕上每个点发射电子束。计算机从图形应用程序或扫描仪发出数字信号到电子, 这些信号控制电子设置的电压强度. 不同 RGB 的强度组合将产生不同的颜色. 电子由偏转电磁场帮助瞄准以确保快速精确地屏幕刷新。

图7-4 彩色显示器的工作原理

颜色模型规定一些属性或原色, 将颜色分解成不同属性的数字化组合. 这就是色彩制式的转换问题。

观察者在屏幕上实际看到的色彩要受色彩制式和屏幕上荧光点数量的影响。. 因此每家计算机屏幕制造商都必须在(R, G, B)数据和国际通行的CIE色彩标准之间进行转换, CIE标准使用三原色, 分别称为X, Y和Z。 针对短余辉荧光点的一类典型转换是

0.610.290.150RX0.350.590.063G=YAα=β. 0.040.120.787BZ计算机程序把用CIE数据(X, Y, Z)表示的色彩信息流发送到屏幕，求屏幕上的电子将这些数据

转换成(R, G, B)数据的方程。现在要根据CIE数据(X, Y, Z)计算对应的(R, G, B)数据, 就是等式A =  中A和 已知, 求。如果A是可逆矩阵, 则由A = 可得 = A1. 解：在Matlab命令窗口输入以下命令

A = [0.61,0.29,0.15;0.35,0.59,0.063;0.04,0.12,0.787]; % B = inv(A),

程序执行的结果为：

B = 2.2586 -1.0395 -0.3473 -1.3495 2.3441 0.0696 0.0910 -0.3046 1.2777

可见将CIE数据(X, Y, Z)转换成(R, G, B)数据的方程为 =B

在彩色电视技术中，还有许多地方会用到线性变换.比如民用电视信号的发送并不直接采用

4

·5· 第2章 时域中的离散信号和系统

(R, G, B)数据，而是使用向量(Y, I, Q)来描述每种颜色.其中Y称为亮度信号，I和Q为色差信号，

这样的制式可以达到彩色和黑白的兼容.如果屏幕是黑白的, 则只用到了Y，这比CIE数据能提供更好的单色图象。YIQ与“标准”RGB色彩之间的对应如下

Y0.2990.5870.114I=0.5960.2750.321Q0.2120.5280.311它的逆变换矩阵留给读者自行计算。

RG B7.5 人员流动问题

某试验性生产线每年一月份进行熟练工与非熟练工的人数统计, 然后将1/6熟练工支援其

他生产部门, 其缺额由招收新的非熟练工补齐。新、老非熟练工经过培训及实践至年终考核有2/5成为熟练工，假设第一年一月份统计的熟练工和非熟练工各占一半，求以后每年一月份统计的熟练工和非熟练工所占百分比。

设第n年一月份统计的熟练工和非熟练工所占百分比分别为xn和yn, 记成向量Zn因为第一年统计的熟练工和非熟练工各占一半, 所以Z1xn. ynx11/2。为了求以后每年的熟=y11/2练工和非熟练工所占百分比, 先求Zn+1与Zn的关系式。根据已知条件可得:

xn+1 = (1

12192)xn +(xn + yn) =xn +yn, 6561053211yn+1 = (1)(xn + yn) =xn +yn,

556109/102/5Zn+1Zn=AZn 1/103/5Zn+1= AZn= A2Zn1= … = AnZ1.

即

将A对角化成为APΛP可得到简明易算的结果：Zn+1= PΛP键入命令[P,D] = eig(A)，得P-1-1nZ1PΛnP-1Z1。

0.97010.7071 10，于是便有： ,DΛ0.24250.7071 00.510.97010.70711n00.97010.70710.50.80.3/(2n)xn1 , 0.5nny 0.24250.70710.24250.707100.50.20.3/(2)n1当n不断增加时，xn1，yn1分别趋向于0.8和0.2。

7.6 二氧化碳分子结构的振动频率

二氧化碳分子可看成中间一个碳原子，左右分别以弹簧（化学键）联接两个氧原子，构成

一个三质量振动系统（见图7-5）。其方程为：

5

·6· 第2章 时域中的离散信号和系统

d2x1mo2kx1kx2dtd2x2mc22kx2kx1kx3

dtd2x3mo2kx2kx3dt设三个原子沿轴向的振动具有同样的频率ω，

图7-5 CO2分子振动结构

xjAje，则代入上式后得到：

kk2A1A1A2

momo2kkk2A2A2A1A3

mcmcmckk2A3A2A3

momo写成矩阵形式为：

itkk2AA012mmoo2kkk2A1A2A30mcmcmckk2A2A30momokk20mmooA12kk k2AO2mcmcmcA3kk20mmoo振动的频率的平方就决定于这个系数矩阵的特征值，因为有三个特征值，要取其中的最大

值。由此写出计算程序pla706，核心语句为： k=14.2e2, mo=16*1.6605e-27, mc=12*1.6605e-27, A=[k/mo,-k/mo,0;-k/mc,2*k/mc,-k/mc;0,-k/mo,k/mo], [x,y]=eig(A), omega=sqrt(y)

lamda=2*pi*3*10^8./max(max(omega)) 其运行结果如下。答案为：

omega = 1.0e+014

及 lamda = 4.257951e-006 m ≈ 4.3μm

7.7 二自由度机械振动

图7-9表示了一个由两个质量、两个弹簧和两个阻尼器构成的二自由度振动系统，今要求在给定初始位置和初始速度下的运动。

解：设x1和x2分别表示两个质量关于它们平衡位置的偏移，则此二自由度振动系统的一般方程为：

图7-8 两自由度机械振动模型

6

·7· 第2章 时域中的离散信号和系统

m1x1(c1c2)x1c2x2(k1k2)x1k2x20m2x2c2(x2x1)k2(x2x1)0可写成矩阵形式： 其中

(7.10.1)

(7.10.2) MX+CX+KX=0

m0ccc2k1k2k2x1 M=1,C12,K,Xc2k20m1c2k2x2(7.10.4)

这是一个四阶的常系数齐次微分方程组，给定其初始位置X0和初始速度Xd0：

x10X0,x20x10x1d0 X0Xd0x20x2d0 (7.10.4)

引进符号xd是为了在MATLAB编程中给导数以变量名的需要。用解析法求这个方程的解

非常麻烦，只有当C=0，即无阻尼时，系统可解耦为两种的振动模态，通常书上只给出解耦简化后的解。

有了MATLAB工具，根本无需设C＝0，也无需解耦，就可以用很简洁的程序求出其数值解。其基本思路是把原始非常化成典型的四个一阶方程构成的状态方程组：

(7.10.5)

At这个方程在初始条件Y=Y0下的解为Y=Y0e。用MATLAB表示为Y=Y0*expm(A*t)。其中expm表示把(A*t)看成矩阵来求其指数。已经知道标量x求指数的方法，求矩阵指数虽然要麻烦一些，但概念是相仿的。我们不必都弄清细节。MATLAB提供了expm, expm1, …等多种函数供用户调用。所以我们只要把Y,Y0和A找到就行。

先把方程（2）写成两个一阶矩阵方程：

Y=AY

IXX0=X

-M\\K-M\\CXd=X=-M\\CXd-M\\KXdXdX=Xd于是

(7.10.6)

0IX0XY=,Y0=,A= (7.10.7) -M\\C-M\\KXdXd0xd1x1对于本题的二自由度系统，X和Xd，所以Y和Y0都是4×1的单列变量；x2xd2由于A中的四个元素都是2×2方阵，A是4×4方阵。对于更多自由度的系统，公式(7)都是正

确的。

下面给出二自由度系统的一个数值例，设m1=1; m2=9; k1 = 4; k2=2; c1和c2可由用户输入。求在初始条件x0 = [1;0]和 xd0 = [0;-1]下，系统的输出x1,x2曲线。

根据上面的模型可以写出程序pla710，运行此程序，输入c1=0.2,c2=0.5所得的结果见图7-9。从中可清楚地看到振动的两种模态。特别是x1的运动反映了两种模态的迭合。给出不同的初始条件，各模态的幅度也会变化。输入c1=0,c2=0所得的结果见图7-10，这时两种振动模态是可

7

·8· 第2章 时域中的离散信号和系统

图7-9二自由度振动的输出波形

图7-10 二自由度振动波形c1=c2=0

c1=0.2，c2=0.8

解耦的。可以看出，用计算机解题时，问题和数据的复杂性对解题的过程并无根本影响。通常我们选择比较简单，且有解析解的数据作为检验程序之用。一旦确定程序正确，它就可用在很复杂的情况。例如作者就曾把b本例题的核心语句用在一个五自由度的系统上，解一个十阶的线性方程组，同样得出满意的结果。

7.8 FIR数字滤波器最优化设计[12]

设给定滤波器的预期幅频特性D(ω)在ωωi(i1,2,…,K)的K个频点上的值，若实际滤波器的脉冲响应的长度为N=2L+1，则其幅特性与预期特性相拟合的方程为：

A(i)d(n)cosniD(i)Di

n0L(i1, 2, …, K)

这K个方程联立。由于A(i)是L1个谐波的线性组合，这方程组中的待定参数d(n)有L1个。

如果KL1，即方程数与未知数的数目相等，则这个线性方程组是适定的，恰好可以解了这些系数。

如果K>L1，即方程数比未知数的数目多，形成了所谓超定方程组。那就不可能找到精确满足这些方程组的系数d(n)，只能找到最近似地满足这些方程的最小二乘解。

如果KeDPd

其中：e为单列K元误差向量，D为预期幅特性在样本点列上的K元单列向量，d为L1元待定系数单列向量，而P则是由cos（nωi）组成的K(L1)的系数矩阵。

cos(L1） e1 D1 1cos(1) d0 e D 1cos() d cos(L2222）1e eK ,D  DK ,P  1cos(K) cos(LK） ,d dL 

前面已经得到它的最小二乘解的公式：

8

·9· 第2章 时域中的离散信号和系统

d=PPT-1PTD

用MATLAB计算此式特别方便，可以表示为dinv（P’*P）*P’*D，也可以更简便地

用左除运算dP\\D来完成。

例7.8 举一个数字例，要求设计长度为N9( 有5个待定系数)的FIR滤波器，要使它在0～π之间的八个频点ω上逼近预期的低通幅特性D：

ω[0, 0.33, 0.67, 1, 1.5, 2, 2.5, 3.14]; D[1,1,1,1,0,0,0,0]; 解：列出方程组的完整形式：

 1 1 1 1  1 1 1 1cos1cos21cos31cos41 D1 Dcos2cos22cos32cos422d(0)cos3cos23cos33cos43 D3d(1) cos4cos24cos34cos44 D4 d(2)P*d=Dcos5cos25cos35cos45 D5d(3)cos6cos26cos36cos46 D6d(4)cos7cos27cos37cos47 D7cos8cos28cos38cos48 D8 看来系数矩阵P的赋值比较麻烦，其实，回想求离散傅里叶变换时的做法，可以利用

列矩阵w’[w1; w2; …; w8]乘行矩阵[0:4]形成一个85矩阵，然后求余弦得到P。这可以用MATLAB语句Pcos(w’*[0:4])轻而易举地列出。因此用下列几行MATLAB程序ag1010就可以方便地完成最小二乘最优滤波器的设计。

N=9; L=floor((N-1)/2); % 列出待求系数序数 w=[0,0.33,0.67,1,1.5,2,2.5,3.14]; % 列出频率向量 D=[1,1,1,1,0,0,0,0]; % 列出预期幅特性向量 P=cos(w’*[0:L]); % 列出P矩阵 d=inv(P'*P)*P'*D’ % 求出待求系数

程序运行的结果为

d = 0.3981 0.6039 0.2137 -0.1205

-0.1506

知道d(n)以后，就可以计算滤波器的频率特性进行校核。得到的幅频特性见图7-11。因为L4，最高的谐波次数为4，在0～π之间幅特性摆动最多两个周期，达到峰值的次数应为四次，这从图上也看得很清楚。

图7-11 设计出的最优化滤波器的频率特性

7.9 弹性梁的柔度矩阵

设简支梁如图7-12所示，在梁的三个位置分别施加力f1，f2和f3后，在该处

产生的综合变形为图示的y1，y2和y3，通常称为挠度。根据虎克定律，在材料未失去弹性的范围内，力与它引起的变形呈线性关系，可以写出完全的矩阵形式为：

9

·10· 第2章 时域中的离散信号和系统 y1d11d12d13f1ydddf 22122232y3d31d32d33f3y=D*f不难从这个矩阵乘法关系中得知矩阵D的物理意义。假如只施加一个力f1，其余两个力f2

和f3为零，则引起的挠度为：y1d11*f1，y2d21*f1，y3d31*f1，如果加的力为一个单位，则在1，2，3三个位置引起的挠度分别为d11，d21和d31。可以用同样的方法来理解其他的d，利用这个概念，可以用实测的方法来得到矩阵D中的各个元素。这些挠度元素愈大，表明这个梁愈柔软，所以矩阵D被称作柔度矩阵。

柔度矩阵的逆就是刚度矩阵K，KD 1，在本例中它也是一个33矩阵。

f1k11k12k13y1fkkky 22122232f3k31k32k33y3其物理意义为造成单位挠度所需要的力。这些力愈大，表明这个梁愈刚劲。

在图7-3所示的情况下，柔度矩阵的9个元素中，没有一个为零。这说明它们之间相互耦合。在位置1施加力不但会引起该处的变形，同时也必然引起位置2和3处的变形；在其他位置加力亦然。如果矩阵变为对角矩阵，那也就意味着互相。位置1加的力只引起位置1处的挠度，常识告诉我们，在图7-3所示的简支梁上是不可能出现这种情况的。

.005.002.001，单位为【公分/牛顿】.002.004.003现在举一个数字的例子，设柔度矩阵D。 .001.003.006（1）设在位置1，2，3处施加的力为30，50和20【牛顿】，试求出其挠度。

（2）设要在位置3处产生0.4【公分】的挠度，其他两处挠度为零，试求应施加的力。 解：列出解此题的Matlab程序pla712

D0.001*[5,2,1;2,4,3;1,3,6] f[30;50;20], yD*f y1[0;0;0.4] Kinv(D), f1K*y1

% 输入柔度矩阵

% 解题（1）,给定力求挠度 % 解题（2）给定挠度 % 求刚度矩阵，求力

程序运行的结果如下：

 0.2700 0.3200题（1）的三个挠度为y 0.3000 2.2373 152.24 33.83刚度矩阵为K=inv(D)152.24 491.52 220.3390 33.83 220.3390 271.18 13.5593题（2）的三个力为f188.1356

 108.4746 最后的力中有负数是合理的，因为题目要求有两个位置上的挠度为零，如果三个力都同向，

造成的变形也同向，那就没有一处的挠度会等于零，必须有反向加的力才行。

可以看到，在刚度矩阵中也出现了负值元素，这是什么原因呢？可从其物理意义去想。对照柔度矩阵的物理意义，刚度矩阵中的第一列元素应该是：使挠度y1达到一个单位，而y2和

10

·11· 第2章 时域中的离散信号和系统

y3都等于零时，在三个位置上应施加的力。即同时在位置1施加k11的力，在位置2施加k21的

力，而在位置3施加k31的力。正如上面的分析，要使两个位置上的挠度为零，这三个力必须有正有负。所以表现为刚度矩阵中每列元素中必定有负值元素。只在一个位置加力（意味着其他两个力为零），在三个位置引起挠度是非常直观，符合常识的，所以读者容易接受柔度矩阵的物理概念；只在一个位置加“挠度”（意味着其他两个挠度为零），在三个位置需要加多大力却需要想像，在常识中大概很少有人经历这种状况，所以刚度矩阵的物理概念就要难理解一些。

7.10 用二次样条函数插值5个点

设给出以下五点坐标数据： x 8 11 15 18 (a) 求通过这些点的二次样条函数； y 5 9 10 8 (b) 求在x=12.7处的y的插值值； (c) 画出插值多项式和这些数据点

解：由于有5个点(n=5)，需要4个样条函数(i=1,...,4)。样条函数的方程为：

22 7 fi(x)aix2bixcii1,2,3,4

四个多项式，每个多项式有三个系数，待求的总共为12个系数。按二次样条函数的规定，设a10，其余11个系数由11个线性方程组成的方程组决定。

其中8个方程由各区段的多项式通过数据点决定，即 i1,f1(x1)a1x12b1x1c1b18c15

2f1(x2)a1x2b1x2c2b111c19i2,f2(x2)a2x22b2x2c2a2112b211c29f2(x3)a2x3b2x3c2a215b215c21022

i3,i4,f3(x3)a3x32b3x3c3a3152b315c310f3(x4)a3x4b3x4c3a318b318c38f4(x4)a4x42b4x4c4a4182b418c482222

f4(x5)a4x5b4x5c4a422b422c47

3个方程由相邻多项式在共同节点处的斜率相等的条件决定。即

i22a1x2b12a2x2b2b12a211b2或b12a211b20

i32a2x3b22a3x3b32a215b22a315b3或2a215b22a315b30

i42a3x4b32a4x3b42a318b32a418b4或2a318b32a418b40

把11个联立方程写成矩阵形式，有：

11

·12· 第2章 时域中的离散信号和系统

81000000000b5b1=1.33331110000000001c1=-5.6667c1900112111000000aa2=-0.27089220015151000000b2b7.29171022c-38.4375 100000015151000c2220000018181000a8用X=A\\B解得a0.05563320000000018181b38b3-2.5000c72c335.000000000000222213010-22-10000000a4a0.062b-2.75000003010-30-10000b44000003610-36-100cc32.250044可以用MATLAB解出这11个系数如右，由此写出这几段多项式方程为： f1(x)1.3333x5.66678x11

f2(x)0.27x27.29x38.44f3(x)0.0556x22.5x35f4(x)0.0625x22.75x37.25(b) 在x=12.7处的y值为：

11x15 15x18 18x22

图7-13 二次样条函数插值结果

f2(x)0.2712.727.2912.738.4410.49

(c) 曲线也要分段画， 画出的曲线如右图，可以看出其第一段是一根直线，对应于a10。

7.11 飞行器三维空间运动的矩阵描述

这里主要举例说明线性代数在刚体空间 运动学中的应用。

机械手、飞行器和机器人都是在空间运动的刚体或多刚体组。对它的运动的描述非常复杂，因为一个刚体要由三个转动和三个平移，即用6个自由度才能定量地决定它的位置。

飞行器在空中可以围绕三个轴旋转。假如它在向北飞行时，机头正对北方，则它围绕铅垂轴的旋转角称为偏航角（Yaw），它描述了飞机左右的偏转，用u表示；围绕翼展轴的旋转 角称为倾斜角（Pitch），它描述了飞机俯仰姿态，图7-14 飞行器姿态角演示 用v表示；围绕机身轴的旋转角称为滚动角（Roll），用w表示；u,v和w三个变量统称为欧拉角，它们完全描述了飞机的姿态。

MATLAB中有一个演示程序quatdemo.m，专门演示这几个姿态角所造成的飞机状态。键入:

quatdemo

屏幕上将出现如图6-7所示的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器；右上方为三个姿态角u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿

12

·13· 第2章 时域中的离散信号和系统

态角: 方位角、俯仰角和倾侧角。左下方还有【静态】和【动态】两个复选钮，我们只介绍【静

态】，读者可自行试用【动态】进行演示。

用键入参数或移动标尺的方法分别给u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。

本例用数值来讨论这个程序的实现方法。先把飞行器的三维图像用N个顶点的三维坐标描述，写成一个3N的数据矩阵G。其顶点次序要适当安排，使得用plot3命令时按顶点连线能绘制出飞行器的外观。例如，以下的程序（pla607前半部分）即可画出一个最简单的飞行器立体图，如图6-8所示。

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]';

图7-15 用G画出的飞行器 尾翼的顶点坐标

G=[Gw,Gt] % 整个飞行器外形的数据集

plot3(Gw(1,:),Gw(2,:),Gw(3,:),'r'),hold on plot3(Gt(1,:),Gt(2,:),Gt(3,:),'g'), axis equal

运行此程序得出整个飞行器外形的数据集为

4 4 0 4 0 0 0 0 3 3 7 3 3 3 2 3 G (7.14.1） 0 0 0 0 0 3 0 0飞行器围绕各个轴旋转的结果，表现为各个顶点坐标发生变化，也就是G的变化。只要把

三种姿态的变换矩阵Y, P和R乘以图形数据矩阵G即可。其中

00cosusinu01cosv0sinv,R0coswsinw,P0 (7.14.2） sinucosu010 Y0100sinwcoswsinv0cosv如果三个姿态角都变化，转动的次序为: 先滚动R，再倾斜P，最后偏航Y，则最后的图像

数据成为Gf  Y G2  Y P G1  YP RG  QG。最后变换矩阵成为

cosucosvsinucoswcosusinvsinwsinusinwcosusinvcoswQsinucosvcosucoswsinusinvsinwcosusinwsinusinvcoswsinvcosvsinwcosvcosw

(7.14.3）

由于矩阵乘法不服从交换律，转动次序不同时结果也不同。这个过程可以用手工计算，也可以用MATLAB程序来实现。程序pla714半部分如下:

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.7071运行结果为 0.35360.6124A0.70710.35360.6124

0.86600.5000我们知道，二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换0 13

·14· 第2章 时域中的离散信号和系统

矩阵的行列式代表的是三个向量组成的平行六面体的体积。不难算出，det(Y)  1，det(R)  1，det(P)  1，det(Q)  1。当然也有det(A)  1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于1。

飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和3×3的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的N个顶点。这些顶点可用三维空间中的三个坐标来描述，也就是3×N数据集G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和4×4变换矩阵。即数据集G4成为4×N矩阵，其最下面一行元素的取值都为1。

G3G4=ones(1,N) (7.14.4）

而平移变换矩阵M和旋转变换矩阵Y, R, P成为

10M40001000010c1c2c31

(7.14.5）

0Y(R,P)0333Y4(R4,P4)00001 (7.14.6）

其中c1,c2,c3为在三个轴x1,x2,x3方向上的平移距离。这种方法也适用于机器人运动学的研究。

7.12 金融公司支付基金的流动

金融机构为保证现金充分支付, 设立一笔总额00万的基金, 分开放置在位于A城和B城的两家公司, 基金在平时可以使用, 但每周末结算时必须确保总额仍然为00万. 经过相当长的一段时期的现金流动, 发现每过一周, A城公司有10%支付基金流动到B城公司, B城公司则有12%支付基金流动到A城公司. 起初A城公司基金为2600万, B城公司基金为2800万. 按此规律, 两公司支付基金数额变化趋势如何? 如果金融专家认为每个公司的支付基金不能少于2200万, 那么是否需要在必要时调动基金?

设第k+1周末结算时, A城公司B城公司的支付基金数分别为ak+1, bk+1 (单位:万元), 则有a0=2600, b0=2800,

ak10.9ak0.12bk. bk10.1ak0.88bk原问题转化为:

(1) 把ak+1, bk+1表示成k的函数, 并确定limak和limbk .

kk(2) 看limak和limbk 是否小于2200.

kk 14

·15· 第2章 时域中的离散信号和系统

ak0.90.122600 ,A,X00.10.882800bk2k+1由此可得：Xk+1=AXk=AXk-1==AX0

解：令Xk为了计算Ak+12600,固然可以用矩阵乘法硬算，但不是高明的办法。较好的方法是将A2800对角化。输入MATLAB命令:[P,D] = eig([0.9,0.12;0.1,0.88])，得到

00.76820.70711.00P,D00.78

0.020.7071因此 Xk+1ak1260001P*00.78k1*inv(P)*2800= bk1键入以下语句即可完成这一计算：Xkp1=P*[D(1,1)^k,0;0,D(2,2)^k]*inv(P)*X0

k1ak132400/110.783800/11, 得到： Xkp1=k1bk127000/110.783800/11可见ak单调递增, bk 单调递减, 且limak =

k32400270002945,limbk =24.

k1111两者都大于2200, 所以不需要调动基金。

7.13 质谱图实验结果分析

某一样本的质谱图给出了一系列的尖峰，

它们表示了样本中所含的各种离子成分的质量。各尖峰的高度Ij取决于不同成分的质量。

峰点 1 j物质的Cij CH4 0.2 28 C2H4 0.2 1 18 C2H6 0.3 0 12 C3H6 0.2 0 2.4 C3H8 0.2 0.1 16 1 2 18 Cniji1N

2 3 其中，Cij是物质i的离子对尖峰j所作贡献，

4 10 0 nj是物质j的离子总数，每个尖峰的系数Cij如右

5 10 表。

6 设某样本产生的质谱图的尖峰为：

I1=3.4, I2=20.5, I3=170, I4=49, I5=39.8, I6=96.3, 试求出样本中各种物质的含量。 解：设五种物质的含量分别为x1,x2,x3,x4,x5，则可列出下列方程组：

0.2*x1+0.2*x2+0.3*x3+0.2*x4+0.2*x5=3.4 28*x1+x2+0.1*x5=20.5

18*x2+12*x3+2.4*x4+16*x5=170 10*x3+x5=49 10*x4+2*x5=39.8 18*x5=96.3

这里有六个方程和五个未知数，是一个超定问题，可用左除求解。写成矩阵形式为：

15

·16· 第2章 时域中的离散信号和系统

0.22800000.20.30.20.23.4x11000.120.5x217018122.416x3

0100149x00102439.8x596.300018解题的程序pla702如下：

A=[0.2,0.2,0.3,0.2,0.2;28,1,0,0,0.1;0,18,12,2.4,16;0,0,10,0,1;0,0,0,10,2;0,0,0,0,18]

B=[3.4;20.5;170;49;39.8;96.3], X=A\\B

解得X的五个分量为：x1=0.6634，x2=1.3909，x3=4.365，x4=2.91，x5=5.35。

7.14 用特征方程解Fibonacci数列问题

Fibonacci数列于公元1200年左右被发现，在现代物理、准晶体结构、化学等领域都有直接的应用，它把0和1作为初始项F0和F1，以后的每一项为其前两项的和，Fk2Fk1Fk，

于是得到数列1,1,2,3,5,8,13,21,…。数列递推的输入是两项，输出是一项，所以需要用两维矩阵建模。为了求出k无限增大时的Fk，相当于研究稳态的趋势，就需要应用特征值和特征向量。这是一个研究矩阵和特征值应用的极好例子。解的过程如下。

一、建立矩阵模型：

Fk2Fk1FkFk1先把Fibonacci规则写成联立方程，设变量为两维向量uk，

FFFk1k1k可写成矩阵方程：

Fk211Fk1uk+1=Auk Fk110FkF1111u其中A，初始条件0F0。 100如果要求出第k项的Fibonacci数，就可以用矩阵连乘的方法，而写成

uk=AAkAu0=Aku0

矩阵连乘运算也很麻烦，注意A是对称实矩阵，把它对角化，可以大大简化计算。令

A=PΛP1，则上式可写成uk=PΛkP-1u0=Aku0，因为对角矩阵乘方就等于主对角线各元

2a0a0a0素自行乘方：，以此类推，所以其计算就方便得多。 20b0b0b二、将系数矩阵对角化

① 建立特征方程并求根。为了将A对角化，首先求A的特征方程和特征值： 令 det(AI)111210

16

·17· 第2章 时域中的离散信号和系统

用二项式定理，得到两个特征根为1151-51.618,2-0.618, 22② 求特征向量，将1代入特征矩阵的系数，进行消元，得出行最简型 （注意121,及12= -1,因而1/(1-1)1）

111/111111111 101/11001100这是秩r=1，n=2的欠定方程，它的一个基础解即特征向量为p1的特征向量为p21。同理得对应于21212112-1因而，其逆，P没有,Pp,pP1212111111规范化并不要紧，因为P-1中的系数可以补偿。于是得出矩阵形式公式：

121ukPΛkP1u012Fk11kuFk12k12101210 11k1102 (a)

也可把矩阵展开，取其第二行，得出标量公式

1k12k11k2k k，第二行为Fk=k1212(b)

三、程序编写及验证

用下列程序pla614可以很快求出任意k阶的Fibonacci数。 A=[1,1;1,0], [p,lamda]=eig(A) for k=[1,3,5,10,20,50,100]

F=p*lamda^k*inv(p)*[1;0]; % 用矩阵公式(a)计算 % 下一语句是用标量公式(b)计算

F1=(lamda(1,1)^k-lamda(2,2)^k)/(lamda(1,1)-lamda(2,2)); [k,F(2,:),F1] % 显示计算结果 end

算出的结果见下表，F(1,:)和F1当然是一样的。 k F 3 5 10 20 50 100 2 5 55 6765 12586269025 3.2248481792631e+020 k=73以上时，Ｆ的有效位数已超过了MATLAB的精度（15位），靠MATLAB已经无法算得精确到整数位的Fibonacci数了。

有趣的是：这样一个完全是自然数的数列，通项公式却是用无理数来表达的。 而且当k趋向于无穷大时，后一项与前一项的比值越来越逼近黄金分割1.618（其倒数是0.618）.。因为λ2小于1，当k无限增大时，公式(b)中的后一项趋于零，Fk的变化趋势仅取决于λ1。Fk+1//Fk的值随k=3到9变化如下：1.5，1.667，1.6，1.625，1.61，1.6190，1.6176，…，愈来愈趋近1.618。

17

·18· 第2章 时域中的离散信号和系统

7.15 简单线性规划问题

线性规划是一门研究如何使用最少的人力、物力和财力去最优地完成科学研究、工业设计、

经济管理中实际问题的专门学科.主要在以下两类问题中得到 应用：一是在人力、物力、财务等资源一定的条件下，如何使用它们来完成最多的任务；二是给一项任务，如何合理安排和规划，能以最少的人力、物力、资金等资源来完成该项任务。它比线性代数多了两个概念：第一、在欠定方程组的基础上，又引进了不等式（包括解必须大于零），从而建立了解的可行域的概念；第二、在可行区内求极大值的概念。最简单的线性规划问题在中学数学中就介绍了。

某工厂用A、B两种配件生产甲、乙两种产品，每生产一件甲产品使用4个A配件并耗100工时，每生产一件乙产品使用4个B配件并耗时200工时，该厂每天最多可从配件厂获得16个A配件和12个B配件，按每天工作800工时计算，该厂所有可能的日生产安排是什么？若生产一件甲产品获利2万元，生产一件

乙产品获利3万元，采用哪种生产安排获得的利润最大？

解：设甲、乙两种产品的日生产分别为x,y件时，工厂获得的利润为z万元，则要求

最大化目标函数z2x3y

而x,y应满足约束条件为

100x200y8004x16, 4y12x,y0

图6-11 例6.15的可行域绘制和目标函数的最大化

在线性规划的问题中，约束条件通常由一组欠定线性方程组和一些不等式构成。在本例中，第一个方程x2y8如果取等式，它就是r=1，n=2的二元欠定方程（组），其解集就是一根斜线x2y80。由于是不等式，其解集（即可行域）是该斜线左下方的半平面。其他三个不等式的约束条件比较简单，作出约束条件所表示的可行域，如图6-11所示。

设目标函数为任意常数c，其等值线2x3yc是一系列斜率为-2/3的平行线。当此直线平移经过可行域时，在点M（4,2）处达到z的最大值.zmax2x3y14。即每天安排生产4件甲产品，2件乙产品时，工厂获利最大为14万元。

规模较大的线性规划问题就需要用到线性代数的理论。它的一般数学模型如下

max(min)zc1x1c2x2cjxjcnxn

a11x1a12x2a1nxn(,)b1AX(,)B s..taxaxax(,)bmnnmm11m22x,x,,x0n12max(min) z表示使目标函数极大（极小）化，s.t.(subject to受约束于)…。可以看出，约

束条件的主体部分就是一个线性方程组，去掉其中的不等式，余下的线性方程组必定是欠定的，这样它才会有无穷多个解并组成一个子空间。加上其中的不等式以后，可行域就会向半平面扩

18

·19· 第2章 时域中的离散信号和系统

大，所以在学线性代数时需要弄清欠定线性方程组的解的特性。约束条件中最后的

x1,x2,,xn0，则是对自变量取值的，通常是默认的。

MATLAB中有解线性规划问题的专用子程序LINPROG(f,A,B)，它用来解决如下的问题： 最小化 z=f*x

约束条件: Ax <= b (x的所有分量为正的条件是默认的，不必列入)

对于例6.14，如果要套用这个函数，就要做一些修改。首先是把求极大改为求极小，那只要把z的系数f号的正负号反转即可。

xzf*X2,32x3y

y约束条件应写成一个完全的矩阵，要把大于改成小于，也可用改变系数符号的方法：

 100 200800x16B

AX 4 0y 0 412于是相应的MATLAB程序pla613n核心语句可写成：

f=[-2,-3], A=[100,200;4,0;0,4], B=[800;16;12],

X=linprog(f,A,B), zmax= -f*X 程序运行的结果为：X,zmax14

42 19