《实用大众线性代数》课程参考文献（MATLAB版）线性代数在科技及工程中的应用实例.doc

第 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

第 2 章时域中的离散信号和系统 ·1· 1 第 7 章在科技及工程中的应用实例 7.1 由拉压杆组成的桁架结构由 13 根拉压杆件组成的桁架结构，如图 7-1 所示，13 个平衡方程已给出，它们来自 6 个中间节点，每个节点有 x,y 两个方向的平衡方程，还有一个整体结构的 y 方向平衡方程。现求其各杆所受的力。解：按照题给方程组改写成矩阵形式，令 1 1 2 2 1 1 cos 14 / 16 ^ 2 14 ^ 2 0.6585 cos 16 / 16 ^ 2 16 ^ 2 0.7071 sin 16 / 16 ^ 2 14 ^ 2 0.7526 k k k    = = + = = = + = = = + = 列方程时假设各杆的受力均为拉力，其相应的方程组及化为矩阵后的形式为： 2 2 1 2 2 6 3 4 1 5 2 1 2 1 3 3 5 7 1 8 4 3 8 9 10 1 5 6 9 3 5 2 11 7 2 11 12 12 3 8 13 2 11 F +k F =0 k 1 0 0 0 0 0 0 0 0 0 0 0 -F +F =0 0 - F =2000 F +k F -k F =0 k F +F +k F =-1000 F +k F -F =0 k F +F = -500 F -k F -F =0 F +k F = 4000 k F -F =0, k F +F =-500 F +k F = 2000 F +k F =0                       2 1 2 3 1 3 1 3 2 2 3 2 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -k 0 0 1 k 0 0 0 0 0 0 0 0 k 0 1 0 k 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 1 k 0 0 0 0 0 0 0 0 0 0 0 0 k 1 0 0 0 0 0 0 0 0 -k -1 0 0 0 1 0 0 0 0 0 0 0 k 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 k 0 0 0000000000k 10 0 0 0 0 0 0 0 k 0 0 0 1 0 0000000000k 01                      1 2 3 4 5 6 7 8 9 10 11 12 13 F 0 F 0 F 2000 F 0 F -1000 F 0 F -500 (7.1.1) F 0 F 4000 F 0 F -500 F 2000 F 0                                         =                                           将它看作 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-1 由拉压杆件组成的桁架结构图 7-2 三阶格型梯形滤波器的结构图

第 2 章时域中的离散信号和系统 ·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%以上的师生恐怕一个小时也解不出来，而且做对的概率极低。用矩阵的思路和方法来解就完全不同，它不是通过消元来减少变量，而是想办法补上所有的零元素，把方程扩充为完整的矩阵形式： 1 3 4 2 1 3 3 2 4 4 7 5 2 2 8 6 5 7 2 1 2 3 4 6 5 6 7 8 9 8 8 11 9 6 1 12 10 9 11 1 10 12 12 10 13 10 11 12 13 - - - x u k x x x x k x x x qx x x k x x x x k x x x qx x x k x x x x k x x x qx x x x x x x x x x x x x x x = = = + = = = =                                                       + = =   = = + = = 3 3 0, 0, 0, - , 0, 0, 0, 0, 0, 0, 0, 0, 0 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, , 0, 1, 0, 0, 0, k k =         2 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, q, 0, 0, 0, 0, 0, 0 0, 1, 0, 0, 0, 0, 0, - , 0, 0, 0, 0, 0 0, 0, 0, k 2 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, , 0, 1, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, q, 0, k 1 0 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, - , 0 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, k 1 3 2 1 0 0, , 0, 1, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, q, 0, 0, 0 0, 0, C , 0, 0, 0, C , 0, 0, 0, C , C , 0 k                     1 2 3 4 5 6 7 8 9 10 11 12 13 1 0 0 0 0 0 0 0 0 0 0 0 0 x x x x x x x u x x x x x x                                         +                                             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。运行此程序就得到： 3 2 3 2 1 3 2 3 2 1 (13) 24 49 42.9 30.8 24 49 42.9 30.8 (13) 8 15 13 24 8 15 13 24 x q q q z z z W u q q q z z z − − − − − − = + + + + + + = = + + + + + + 用矩阵模型解信号流图的最大优点是一步到位，依靠计算机，既快速，又极易查错。 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章时域中的离散信号和系统 0≤k≤N-1 其中W=xp(2πj/N)是一个相角2π/N为的单位向量，也称为旋转因子。()就称为x()的离散傅里叶变换(DFT),也就是x)的频谱。用矩阵来表示，可写成： 1 1 1 X(0 型形「x0) O) =1 部形w -W.x X(N-1] W(N-1IN-1 x(N-1) W 所以信号频谱的计算，可以简单地用一个矩阵乘法来完成。信号是NX1维数组，矩阵 W称为DFT矩阵，它是NXN维的复数阵。把矩阵乘法看作一个变换，我们就可以把频谱计算看作信号从时域到频域的线性变换。这个矩阵运算可用下列几条MATLAB语句来实现，程序名为pla7O4。 m=input xn=长度为N的数组) .N=1 ongth(m)%输入据 n=[0:1:N-1]:k=[0:1:N-1: %没定n和k的行向量 J*21/0 Nk=N.·nk %求出W 含k值的N乘N维的整数矩阵从=Xn*Nk 多求出离散傅里叶级数系数两条语句无须解释。第三条语句把n转为列向量'，再与行向量k做矩阵乘，它产生的是整数矩阵： 0 2 [0.12.N-11■ 0 2x(N-D) W-1 0N-1.N-1D×N-1D 这个整数方阵恰好是算式X=Wx中W矩阵的指数部分第四条语句就产生了矩阵，而最后一条语句就完成了DFS的全部运算。由于实际上离散傅里叶级数并不太用到。它已被下节将介绍的离散傅里叶变换取代，而且两者的计算程序是完全相同的。所以如果写成 MATLAB子程序，它可以命名为DFT。在这5行DFT子程序前面，加上输入语句： n=input ('xn=')N=length(xn); 在子程序后面，加上绘图语句： 0.5 就得到本题的程序p1a604。 100 200 100 运行D1a604, 并按提示输入x,即可在子图1 50 上得到x序列的波形 ,并在子图2 上得到此信号的幅频特性。例如输入序列为 100 200 Xn=[ones(1,64),zeros(1,192)],则得到的时域图7-3序列x()及其频诺波形及其频谱图如图7-3。计算每一个DFT分量()，需要做W次复数乘法和：1次复数加法。而要完成整个 DFT,求出K,则需要做XW次复数乘法和（心1）次复数加法。另外，它占的数据存储量为仪M。在频谱计算中，N的典型值是1024。所以参与运算的数据达百万，需要的存储量

第 2 章时域中的离散信号和系统 ·3· 3 1 ( ) ( ) 0 1 N kn N n k x n W k N − X =   −  其中 W j N N = exp 2 / (  ) 是一个相角 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); 在子程序后面，加上绘图语句： subplot(2,1,1),stem(xn,'.'); subplot(2,1,2),stem(abs(Xk),'.') 就得到本题的程序 pla604。运行 pla604，并按提示输入 xn，即可在子图 1 上得到 xn 序列的波形，并在子图 2 上得到此信号的幅频特性。例如输入序列为 xn=[ones(1,64),zeros(1,192)]，则得到的时域波形及其频谱图如图 7-3。计算每一个 DFT 分量 X(i)，需要做 N 次复数乘法和 N-1 次复数加法。而要完成整个 DFT，求出 X，则需要做 N×N 次复数乘法和 N(N-1)次复数加法。另外，它占的数据存储量为 N×N。在频谱计算中，N 的典型值是 1024。所以参与运算的数据达百万，需要的存储量 0 100 200 300 0 0.5 1 0 100 200 300 0 50 100 图 7-3 序列 x(n)及其频谱

第2章时域中的离散信号和系统超过16M,乘法运算的次数超过数百万。因此，在低档的微机上当N较大时，MATLAB 会拒绝执行这个程序，并警告内存不足。实际上，MATLAB内置的是加快计算并节省内存快速傅里叶变换(FFT程序。从这里也可以看出，矩阵给我们提供了一个组织海量数据进行运算的最好方法，对上百万数据进行赋值、运算和绘图，只用了几行程序。线性代数的重要性之所以在近20年来可与微积分相妣美，这是一个重要原因。如果学线性代数而不涉及工程计算实践，那就无法真正体会线性代数的精髓所在。 7.4显示器色彩制式转换问题彩色显示器使用红(R)、绿(G)和蓝(B)光的叠成效应生成颜色。显示器屏幕的内表面由微粒象素组成，每个微粒包括三个荧光点：红、绿、蓝。电子枪位于屏幕的后方，向屏幕上每个点发射电子束。计算机从图形应用程序或扫描仪发出数字信号到电子枪，这些信号控制电子枪设置的电压强度.不同RGB的强度组合将产生不同的颜色.电子枪由偏转电磁场帮助瞄准以确保快速精确地屏幕刷新每个微粒包括三个荧光点：、蓝 8 电子枪吕显示器屏幕的内表数字信号面由微粒象素组成图形应用程序或扫描仪图74彩色显示器的工作原理颜色模型规定一些属性或原色，将颜色分解成不同属性的数字化组合.这就是色彩制式的转换问题。观察者在屏幕上实际看到的色彩要受色彩制式和屏幕上荧光点数量的影响。，因此每家计算机屏幕制造商都必须在(R,G,B)数据和国际通行的C亚色彩标准之间进行转换C正标准使用三原色，分别称为X,y和Z。针对短余辉荧光点的一类典型转换是 (0.610.290.150)R)X) 0.350.59 0.063 G →Aa=B 0.040.120.787JBz 计算机程序把用CE数据(X)表示的色彩信息流发送到屏幕，求屏幕上的电子枪将这些数据转换成(R,G,B)数据的方程。现在要根据CE数据(XYZ)计算对应的(R,G,B)数据，就是等式 Aa=B中A和B已知，求a。如果A是可逆矩阵，则由Aa=B可得a=AB. 解：在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.3040 1.277 可见将C正数整米花6的民用电视信号的发送并不直接采用方程为a=Bp 4

第 2 章时域中的离散信号和系统 ·4· 4 超过 16M，乘法运算的次数超过数百万。因此，在低档的微机上，当 N 较大时，MATLAB 会拒绝执行这个程序，并警告内存不足。实际上，MATLAB 内置的是加快计算并节省内存快速傅里叶变换(FFT)程序。从这里也可以看出，矩阵给我们提供了一个组织海量数据进行运算的最好方法，对上百万数据进行赋值、运算和绘图，只用了几行程序。线性代数的重要性之所以在近 20 年来可与微积分相妣美，这是一个重要原因。如果学线性代数而不涉及工程计算实践，那就无法真正体会线性代数的精髓所在。 7.4 显示器色彩制式转换问题彩色显示器使用红(R)、绿(G)和蓝(B)光的叠成效应生成颜色。显示器屏幕的内表面由微粒象素组成, 每个微粒包括三个荧光点: 红、绿、蓝。电子枪位于屏幕的后方, 向屏幕上每个点发射电子束。计算机从图形应用程序或扫描仪发出数字信号到电子枪, 这些信号控制电子枪设置的电压强度. 不同 RGB 的强度组合将产生不同的颜色. 电子枪由偏转电磁场帮助瞄准以确保快速精确地屏幕刷新。图 7-4 彩色显示器的工作原理颜色模型规定一些属性或原色, 将颜色分解成不同属性的数字化组合. 这就是色彩制式的转换问题。观察者在屏幕上实际看到的色彩要受色彩制式和屏幕上荧光点数量的影响。. 因此每家计算机屏幕制造商都必须在(R, G, B)数据和国际通行的 CIE 色彩标准之间进行转换, CIE 标准使用三原色, 分别称为 X, Y 和 Z。针对短余辉荧光点的一类典型转换是 0.61 0.29 0.150 0.35 0.59 0.063 0.04 0.12 0.787           R G B           = X Y Z            Aα = β . 计算机程序把用 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 在彩色电视技术中，还有许多地方会用到线性变换.比如民用电视信号的发送并不直接采用