一、引言

  矩阵是线性代数的核心工具，在工程实践中发挥着不可替代的作用。从机械工程的振动分析到电气工程的电路计算，从计算机图形学的三维变换到机器学习的数据处理，矩阵无处不在。

  本文系统梳理矩阵在三大工程领域（机械、电气、计算机）中的具体应用，每个应用都配有详细的计算步骤，帮助读者建立理论与实践的桥梁。

二、矩阵的基本回顾

2.1 矩阵的定义与表示

  由 $m\times n$ 个数排成的 $m$ 行 $n$ 列的数表称为 $m\times n$ 矩阵：

$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]$$

  记作 $A_{m\times n}=(a_{ij})$，其中 $i=1,2,\cdots ,m$，$j=1,2,\cdots ,n$。

a. 特殊矩阵类型

类型	定义	示例
行矩阵	$1\times n$ 矩阵	$[1,2,3]$
列矩阵	$m\times 1$ 矩阵	$\left[\begin{array}{c}1\\ 2\\ 3\end{array}\right]$
方阵	$m=n$	$A_{n\times n}$
对角矩阵	$a_{ij}=0,i\ne j$	$\mathrm{diag}(a_1,a_2,\cdots ,a_n)$
单位矩阵	$a_{ij}=\{\begin{array}{ll}1& i=j\\ 0& i\ne j\end{array}$	$I_n$
零矩阵	所有元素为0	$O_{m\times n}$

2.2 矩阵的核心运算

a. 矩阵加法

  两个同型矩阵对应元素相加：

$$C_{ij}=A_{ij}+B_{ij}$$

b. 矩阵乘法

  $A_{m\times p}$ 与 $B_{p\times n}$ 的乘积 $C_{m\times n}$：

$$c_{ij}=\sum _{k=1}^p{a}_{ik}\cdot b_{kj}$$

c. 矩阵转置

  行列互换：$(A^T{)}_{ij}=A_{ji}$

d. 矩阵求逆

  若 $A{A}^{-1}=A^{-1}A=I$，则 $A^{-1}$ 为 $A$ 的逆矩阵。

  计算公式（伴随矩阵法）：

$$A^{-1}=\frac{1}{\det (A)}\mathrm{adj}(A)$$

2.3 矩阵在工程中的桥梁作用

  矩阵之所以在工程中广泛应用，源于以下特性：

特性	工程意义
批量运算	一次矩阵运算等价于大量标量运算，适合计算机批处理
线性表示	线性方程组、线性变换都可用矩阵统一描述
坐标变换	刚体运动、坐标变换都归结为矩阵乘法
特征分析	特征值反映系统的固有特性（频率、振型等）

三、机械工程中的应用

3.1 机械振动分析

a. 多自由度系统建模

  对于具有 $n$ 个自由度的振动系统，其运动方程为：

$$M\ddot{x}+C\dot{x}+Kx=F(t)$$

  其中：

• $M$：质量矩阵（$n\times n$ 对称正定矩阵）
• $C$：阻尼矩阵（$n\times n$ 对称或非对称矩阵）
• $K$：刚度矩阵（$n\times n$ 对称正定矩阵）
• $x$：位移向量（$n\times 1$）
• $F(t)$：激振力向量（$n\times 1$）

b. 质量矩阵、刚度矩阵、阻尼矩阵

【计算示例】二自由度弹簧-质量系统

    k₁     k₂
---|---|---|||---|---|||
   m₁     m₂

  质量矩阵（集中质量）：

$$M=\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]$$

  刚度矩阵（根据力平衡推导）：

$$K=\left[\begin{array}{cc}{k}_1+k_2& -k_2\\ -k_2& k_2\end{array}\right]$$

  设 $m_1=1\text{kg}$，$m_2=2\text{kg}$，$k_1=1000\text{N/m}$，$k_2=500\text{N/m}$：

$$M=\left[\begin{array}{cc}1& 0\\ 0& 2\end{array}\right],K=\left[\begin{array}{cc}1500& -500\\ -500& 500\end{array}\right]$$

c. 特征值问题与固有频率求解

  对于无阻尼自由振动：$M\ddot{x}+Kx=0$

  设解的形式为 $x=\varphi e^{i\omega t}$，代入得广义特征值问题：

$$(K-{\omega }^{2}M)\varphi =0$$

  有非零解的条件：$\det (K-{\omega }^{2}M)=0$

【计算步骤】

  第一步：列出特征方程

$$\det \left[\begin{array}{cc}1500-{\omega }^2\cdot 1& -500\\ -500& 500-{\omega }^2\cdot 2\end{array}\right]=0$$

  第二步：计算行列式

$$(1500-{\omega }^2)(500-2{\omega }^2)-(-500)(-500)=0$$

$$750000-3000{\omega }^2-2{\omega }^4+500{\omega }^2-250000=0$$

$$-2{\omega }^4-2500{\omega }^2+500000=0$$

$${\omega }^4+1250{\omega }^2-250000=0$$

  第三步：令 $\lambda ={\omega }^2$，解二次方程

$${\lambda }^2+1250\lambda -250000=0$$

$$\lambda =\frac{-1250\pm \sqrt{1250^2+4\times 250000}}{2}$$

$$=\frac{-1250\pm \sqrt{1562500+1000000}}{2}=\frac{-1250\pm \sqrt{2562500}}{2}$$

$$=\frac{-1250\pm 1600.78}{2}$$

$${\lambda }_1=175.39,{\lambda }_2=-1425.39$$

  第四步：求固有频率

$${\omega }_1=\sqrt{{\lambda }_1}=\sqrt{175.39}=13.24\text{rad/s}$$

$${\omega }_2=\sqrt{1425.39}=37.75\text{rad/s}$$

  固有频率（Hz）：

$$f_1=\frac{{\omega }_1}{2\pi }=2.11\text{Hz},f_2=\frac{{\omega }_2}{2\pi }=6.01\text{Hz}$$

3.2 机器人运动学

a. 齐次坐标变换矩阵

  三维空间中一点 $P$ 的齐次坐标为 $\tilde{P}=[x,y,z,1{]}^T$。

  基本变换矩阵：

平移矩阵（沿 $x,y,z$ 轴平移 $t_x,t_y,t_z$）：

$$T_{\text{translate}}=\left[\begin{array}{cccc}1& 0& 0& t_x\\ 0& 1& 0& t_y\\ 0& 0& 1& t_z\\ 0& 0& 0& 1\end{array}\right]$$

旋转矩阵（绕 $z$ 轴旋转 $\theta$）：

$$R_z(\theta )=\left[\begin{array}{cccc}\cos \theta & -\sin \theta & 0& 0\\ \sin \theta & \cos \theta & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

缩放矩阵（沿各轴缩放 $s_x,s_y,s_z$）：

$$T_{\text{scale}}=\left[\begin{array}{cccc}{s}_x& 0& 0& 0\\ 0& s_y& 0& 0\\ 0& 0& s_z& 0\\ 0& 0& 0& 1\end{array}\right]$$

b. 正运动学求解

  机器人末端执行器的位姿可通过齐次变换矩阵链计算：

$$T_{\text{end}}=T_0^1\cdot T_1^2\cdot T_2^3\cdots T_{n-1}^n$$

【计算示例】二自由度机械臂

      θ₂
       ↗
   L₂ ↗  P (末端)
       │
   L₁ │
       │
   ───┴── θ₁
     基座

  已知：$L_1=4\text{m}$，$L_2=3\text{m}$，${\theta }_1=30°$，${\theta }_2=45°$

  DH参数表：

连杆 $i$	${\theta }_i$	$d_i$	$a_i$	${\alpha }_i$
1	30°	0	4	0°
2	45°	0	3	0°

  第一步：计算各连杆变换矩阵

$$T_0^1=\left[\begin{array}{cccc}\cos 30°& -\sin 30°& 0& 4\cos 30°\\ \sin 30°& \cos 30°& 0& 4\sin 30°\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}0.866& -0.5& 0& 3.464\\ 0.5& 0.866& 0& 2\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

$$T_1^2=\left[\begin{array}{cccc}\cos 45°& -\sin 45°& 0& 3\cos 45°\\ \sin 45°& \cos 45°& 0& 3\sin 45°\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}0.707& -0.707& 0& 2.121\\ 0.707& 0.707& 0& 2.121\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

  第二步：计算末端位姿

$$T_0^2=T_0^1\cdot T_1^2=\left[\begin{array}{cccc}0.866& -0.5& 0& 3.464\\ 0.5& 0.866& 0& 2\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}0.707& -0.707& 0& 2.121\\ 0.707& 0.707& 0& 2.121\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

  第三步：矩阵乘法计算

  $T_0^2(0,3)=0.866\times 2.121+(-0.5)\times 2.121+0\times 0+3.464\times 1=1.837-1.061+3.464=4.24$

  $T_0^2(1,3)=0.5\times 2.121+0.866\times 2.121+0\times 0+2\times 1=1.061+1.837+2=4.898$

$$T_0^2=\left[\begin{array}{cccc}0.259& -0.966& 0& 4.24\\ 0.966& 0.259& 0& 4.90\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

  结果：末端位置为 $(4.24\text{m},4.90\text{m},0\text{m})$

c. 雅可比矩阵与速度分析

  雅可比矩阵描述关节速度与末端速度的关系：

$$\dot{x}=J(q)\dot{q}$$

  对于平面二连杆机械臂：

$$J(q)=\left[\begin{array}{cc}-L_1\sin {\theta }_1-L_2\sin ({\theta }_1+{\theta }_2)& -L_2\sin ({\theta }_1+{\theta }_2)\\ L_1\cos {\theta }_1+L_2\cos ({\theta }_1+{\theta }_2)& L_2\cos ({\theta }_1+{\theta }_2)\end{array}\right]$$

【计算示例】

  已知 ${\theta }_1=30°$，${\theta }_2=45°$，$L_1=4\text{m}$，$L_2=3\text{m}$，${\dot{\theta }}_1=0.5\text{rad/s}$，${\dot{\theta }}_2=0.3\text{rad/s}$

  雅可比矩阵：

$$J=\left[\begin{array}{cc}-4\sin 30°-3\sin 75°& -3\sin 75°\\ 4\cos 30°+3\cos 75°& 3\cos 75°\end{array}\right]$$

$$=\left[\begin{array}{cc}-4\times 0.5-3\times 0.966& -3\times 0.966\\ 4\times 0.866+3\times 0.259& 3\times 0.259\end{array}\right]$$

$$=\left[\begin{array}{cc}-2-2.898& -2.898\\ 3.464+0.777& 0.777\end{array}\right]=\left[\begin{array}{cc}-4.898& -2.898\\ 4.241& 0.777\end{array}\right]$$

  末端速度：

$$\dot{x}=J\cdot \dot{q}=\left[\begin{array}{cc}-4.898& -2.898\\ 4.241& 0.777\end{array}\right]\left[\begin{array}{c}0.5\\ 0.3\end{array}\right]$$

$$=\left[\begin{array}{c}-4.898\times 0.5+(-2.898)\times 0.3\\ 4.241\times 0.5+0.777\times 0.3\end{array}\right]=\left[\begin{array}{c}-3.32\\ 2.36\end{array}\right]\text{m/s}$$

3.3 有限元分析（FEA）

a. 单元刚度矩阵组装

  平面应力问题的三角形单元刚度矩阵：

$$[K{]}^e={\iint }_A[B{]}^T[D][B]tdA$$

  其中 $[B]$ 为应变-位移矩阵，$[D]$ 为弹性矩阵。

【计算示例】三角形单元

    3 (0,2)
       /\
      /  \
     /    \
  1 /______\ 2
   (0,0)   (3,0)

  节点坐标：$(0,0)$、$(3,0)$、$(0,2)$，厚度 $t=0.1\text{m}$，弹性模量 $E=2\times 10^{11}\text{Pa}$，泊松比 $\nu =0.3$

  第一步：计算单元面积

$$A=\frac{1}{2}\left|\begin{array}{ccc}1& 0& 1\\ 1& 3& 1\\ 1& 0& 2\end{array}\right|=\frac{1}{2}(1\times 3\times 2+0\times 1\times 1+1\times 0\times 0-1\times 3\times 1-0\times 1\times 2-1\times 0\times 0)=3{\text{m}}^2$$

  第二步：建立几何矩阵 $[B]$

$$[B]=\frac{1}{2A}\left[\begin{array}{cccccc}{y}_2-y_3& 0& y_3-y_1& 0& y_1-y_2& 0\\ 0& x_3-x_2& 0& x_1-x_3& 0& x_2-x_1\\ x_3-x_2& y_2-y_3& x_1-x_3& y_3-y_1& x_2-x_1& y_1-y_2\end{array}\right]$$

$$=\frac{1}{6}\left[\begin{array}{cccccc}0-2& 0& 2-0& 0& 0-0& 0\\ 0& 3-0& 0& 0-0& 0& 0-3\\ 3-0& 0-2& 0-0& 0-2& 0-3& 0-0\end{array}\right]=\frac{1}{6}\left[\begin{array}{cccccc}-2& 0& 2& 0& 0& 0\\ 0& 3& 0& 0& 0& -3\\ 3& -2& 0& -2& -3& 2\end{array}\right]$$

  第三步：建立弹性矩阵（平面应力）

$$[D]=\frac{E}{1-{\nu }^2}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]=\frac{2\times 10^{11}}{1-0.09}\left[\begin{array}{ccc}1& 0.3& 0\\ 0.3& 1& 0\\ 0& 0& 0.35\end{array}\right]$$

$$=2.198\times 10^{11}\left[\begin{array}{ccc}1& 0.3& 0\\ 0.3& 1& 0\\ 0& 0& 0.35\end{array}\right]$$

  第四步：计算单元刚度矩阵

$$[K{]}^e=[B{]}^T[D][B]\cdot t\cdot A$$

$$=\frac{1}{6}{\left[\begin{array}{cccccc}-2& 0& 2& 0& 0& 0\\ 0& 3& 0& 0& 0& -3\\ 3& -2& 0& -2& -3& 2\end{array}\right]}^T\cdot [D]\cdot \frac{1}{6}\left[\begin{array}{cccccc}-2& 0& 2& 0& 0& 0\\ 0& 3& 0& 0& 0& -3\\ 3& -2& 0& -2& -3& 2\end{array}\right]\times 0.1\times 3$$

b. 整体刚度矩阵

  将所有单元刚度矩阵按节点编号进行组装：

$$[K{]}_{\text{global}}=\sum _e[K{]}^e$$

c. 节点位移与节点力求解

  整体平衡方程：

[K]{u}={F}[K]\{u\} = \{F\}[K]{u}={F}

  引入边界条件后求解：

{u}=[K]−1{F}\{u\} = [K]^{-1}\{F\}{u}=[K]−1{F}

3.4 机械优化设计

a. 灵敏度分析

  目标函数对设计变量的敏感度：

$$S_{ij}=\frac{\partial f}{\partial x_i}\cdot \frac{x_i}{f}$$

b. 海森矩阵与优化算法

  海森矩阵（二阶导数矩阵）：

$$H(x)=\left[\begin{array}{ccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots \\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right]$$

  牛顿法迭代公式：

$$x^{k+1}=x^k-H(x^k{)}^{-1}\nabla f(x^k)$$

四、电气工程中的应用

4.1 电路分析

a. 节点电压法

  对每个独立节点建立KCL方程，写成矩阵形式：

$$[G][V]=[I]$$

  其中 $[G]$ 为节点电导矩阵，$[V]$ 为节点电压向量，$[I]$ 为注入电流向量。

【计算示例】简单电路

    1Ω      2Ω
①──[R1]───[R2]───②
    |       |
   [I₁]    [R₃]
    |       | 1Ω
    |       |
    └───────┘
       ③

  已知：$R_1=1\Omega$，$R_2=2\Omega$，$R_3=1\Omega$，$I_1=3\text{A}$，电流源从节点①流入

  第一步：建立节点电导矩阵

  对节点①、②、③列写KCL方程（节点③接地，$V_3=0$）：

$$G_{11}=\frac{1}{R_1}+\frac{1}{R_2}=1+0.5=1.5\text{S}$$

$$G_{12}=-\frac{1}{R_2}=-0.5\text{S}$$

$$G_{22}=\frac{1}{R_2}+\frac{1}{R_3}=0.5+1=1.5\text{S}$$

$$[G]=\left[\begin{array}{cc}1.5& -0.5\\ -0.5& 1.5\end{array}\right]$$

  第二步：建立电流向量

$$[I]=\left[\begin{array}{c}3\\ 0\end{array}\right]$$

  第三步：求解节点电压

$$\left[\begin{array}{cc}1.5& -0.5\\ -0.5& 1.5\end{array}\right]\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right]=\left[\begin{array}{c}3\\ 0\end{array}\right]$$

  使用克莱姆法则：

$$V_1=\frac{\left|\begin{array}{cc}3& -0.5\\ 0& 1.5\end{array}\right|}{\left|\begin{array}{cc}1.5& -0.5\\ -0.5& 1.5\end{array}\right|}=\frac{4.5}{2}=2.25\text{V}$$

$$V_2=\frac{\left|\begin{array}{cc}1.5& 3\\ -0.5& 0\end{array}\right|}{2}=\frac{-(-1.5)}{2}=0.75\text{V}$$

b. 网孔电流法

  对每个网孔建立KVL方程：

$$[R][I]=[V]$$

【计算示例】

    R₁ = 2Ω
  ┌──────────────┐
  │              │
  │   I₁     I₂  │
  │    ⟲     ⟲   │
  │              │
R₃=1Ω       R₂=3Ω│
  │              │
  │              │
  └──────────────┘
       V = 12V

  第一步：建立网孔电阻矩阵

$$R_{11}=R_1+R_2+R_3=2+3+1=6\Omega$$

$$R_{12}=R_{21}=-(R_2+R_3)=-4\Omega$$

$$[R]=\left[\begin{array}{cc}6& -4\\ -4& 6\end{array}\right]$$

  第二步：建立电压向量

$$[V]=\left[\begin{array}{c}12\\ 0\end{array}\right]$$

  第三步：求解网孔电流

$$[I]=[R{]}^{-1}[V]=\frac{1}{36-16}\left[\begin{array}{cc}6& 4\\ 4& 6\end{array}\right]\left[\begin{array}{c}12\\ 0\end{array}\right]=\frac{1}{20}\left[\begin{array}{c}72\\ 48\end{array}\right]=\left[\begin{array}{c}3.6\\ 2.4\end{array}\right]\text{A}$$

c. 导纳矩阵与阻抗矩阵

  节点导纳矩阵 $[Y]$ 的元素：

$$Y_{ij}=\{\begin{array}{ll}{\sum }_k{y}_{ik}& i=j\\ -y_{ij}& i\ne j\end{array}$$

  阻抗矩阵 $[Z]=[Y{]}^{-1}$

d. 传输线矩阵

  传输线的ABCD参数矩阵：

$$\left[\begin{array}{c}{V}_1\\ I_1\end{array}\right]=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}{V}_2\\ I_2\end{array}\right]$$

4.2 电力系统分析

a. 节点导纳矩阵

  电力系统节点导纳矩阵是稀疏矩阵，定义为：

$$Y_{ij}=\{\begin{array}{ll}{\sum }_{k\ne i}{y}_{ik}& i=j\\ -y_{ij}& i\ne j\end{array}$$

【计算示例】三节点系统

    ①──────②
     \     /
      \   /
       \ /
        ③

  线路参数（标幺值）：

• 线路①②：$y_{12}=-j10$
• 线路②③：$y_{23}=-j8$
• 线路①③：$y_{13}=-j6$

  节点导纳矩阵：

$$Y=\left[\begin{array}{ccc}{Y}_{11}& Y_{12}& Y_{13}\\ Y_{21}& Y_{22}& Y_{23}\\ Y_{31}& Y_{32}& Y_{33}\end{array}\right]=\left[\begin{array}{ccc}-j10-j6& j10& j6\\ j10& -j10-j8& j8\\ j6& j8& -j6-j8\end{array}\right]$$

$$=\left[\begin{array}{ccc}-j16& j10& j6\\ j10& -j18& j8\\ j6& j8& -j14\end{array}\right]$$

b. 潮流计算（牛顿-拉夫逊法）

  牛顿-拉夫逊法的修正方程：

$$\left[\begin{array}{c}\Delta P\\ \Delta Q\end{array}\right]=\left[\begin{array}{cc}H& N\\ J& L\end{array}\right]\left[\begin{array}{c}\Delta \theta \\ \Delta V\end{array}\right]$$

  其中雅可比矩阵元素：

$$H_{ij}=\frac{\partial P_i}{\partial {\theta }_j},N_{ij}=V_j\frac{\partial P_i}{\partial V_j}$$

$$J_{ij}=\frac{\partial Q_i}{\partial {\theta }_j},L_{ij}=V_j\frac{\partial Q_i}{\partial V_j}$$

c. 短路电流计算

  三相短路电流：

$$I_{sc}=[Y_{\text{bus}}{]}^{-1}[V_{\text{prefault}}]$$

4.3 信号处理

a. 离散傅里叶变换矩阵

  $N$ 点DFT的变换矩阵：

$$[W{]}_N=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& w_N& w_N^2& \cdots & w_N^{N-1}\\ 1& w_N^2& w_N^4& \cdots & w_N^{2(N-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& w_N^{N-1}& w_N^{2(N-1)}& \cdots & w_N^{(N-1{)}^2}\end{array}\right]$$

  其中 $w_N=e^{-j2\pi /N}$

b. 协方差矩阵与相关矩阵

  协方差矩阵：

$$[\Sigma ]=E\{(x-\mu )(x-\mu {)}^T\}$$

c. 自适应滤波（最小均方算法）

  LMS算法权重更新：

$$w(n+1)=w(n)+\mu x(n)e(n)$$

  矩阵形式：

$$w(n+1)=w(n)+\mu x(n)[d(n)-w^T(n)x(n)]$$

4.4 控制理论

a. 状态空间表示法

  线性系统的状态空间方程：

$$\dot{x}=Ax+Bu$$

$$y=Cx+Du$$

  其中：

• $A$：系统矩阵（$n\times n$）
• $B$：输入矩阵（$n\times p$）
• $C$：输出矩阵（$q\times n$）
• $D$：直接传递矩阵（$q\times p$）

b. 可控性与可观性矩阵

可控性矩阵：

$$[c]=\left[\begin{array}{ccccc}B& AB& A^{2}B& \cdots & A^{n-1}B\end{array}\right]$$

可观性矩阵：

$$[o]=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^{n-1}\end{array}\right]$$

【计算示例】判断可控性

  系统矩阵和输入矩阵：

$$A=\left[\begin{array}{cc}0& 1\\ -2& -3\end{array}\right],B=\left[\begin{array}{c}0\\ 1\end{array}\right]$$

  第一步：计算 $AB$

$$AB=\left[\begin{array}{cc}0& 1\\ -2& -3\end{array}\right]\left[\begin{array}{c}0\\ 1\end{array}\right]=\left[\begin{array}{c}1\\ -3\end{array}\right]$$

  第二步：构造可控性矩阵

$$[c]=\left[\begin{array}{cc}B& AB\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 1& -3\end{array}\right]$$

  第三步：计算行列式

$$\det ([c])=0\times (-3)-1\times 1=-1\ne 0$$

  结论：系统完全可控

c. 状态反馈与极点配置

  状态反馈律 $u=-Kx$，闭环系统：

$$\dot{x}=(A-BK)x$$

  极点配置到期望位置 $s_i$：

$$\det (sI-A+BK)=\prod _{i=1}^n(s-s_i)$$

4.5 通信系统

a. 多输入多输出（MIMO）信道矩阵

  $N_t$ 发射天线、$N_r$ 接收天线的MIMO信道：

$$y=Hx+n$$

  其中 $H$ 为 $N_r\times N_t$ 信道矩阵，$x$ 为发送信号向量，$y$ 为接收信号向量。

【计算示例】2×2 MIMO

$$H=\left[\begin{array}{cc}0.8& 0.2\\ 0.3& 0.9\end{array}\right],x=\left[\begin{array}{c}1\\ 0\end{array}\right],n=\left[\begin{array}{c}0.1\\ -0.05\end{array}\right]$$

$$y=\left[\begin{array}{cc}0.8& 0.2\\ 0.3& 0.9\end{array}\right]\left[\begin{array}{c}1\\ 0\end{array}\right]+\left[\begin{array}{c}0.1\\ -0.05\end{array}\right]=\left[\begin{array}{c}0.9\\ 0.25\end{array}\right]$$

b. 扩频通信中的正交矩阵

   Walsh-Hadamard 矩阵用于扩频码：

$$H_2=\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right],H_4=\left[\begin{array}{cc}{H}_2& H_2\\ H_2& -H_2\end{array}\right]$$

五、计算机工程中的应用

5.1 计算机图形学

a. 三维变换矩阵

旋转变换（绕 $x$ 轴旋转 $\alpha$）：

$$R_x(\alpha )=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \alpha & -\sin \alpha & 0\\ 0& \sin \alpha & \cos \alpha & 0\\ 0& 0& 0& 1\end{array}\right]$$

旋转变换（绕 $y$ 轴旋转 $\beta$）：

$$R_y(\beta )=\left[\begin{array}{cccc}\cos \beta & 0& \sin \beta & 0\\ 0& 1& 0& 0\\ -\sin \beta & 0& \cos \beta & 0\\ 0& 0& 0& 1\end{array}\right]$$

旋转变换（绕 $z$ 轴旋转 $\gamma$）：

$$R_z(\gamma )=\left[\begin{array}{cccc}\cos \gamma & -\sin \gamma & 0& 0\\ \sin \gamma & \cos \gamma & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

b. 投影矩阵

透视投影矩阵：

$$P_{\text{perspective}}=\left[\begin{array}{cccc}\frac{1}{\tan (\frac{fov}{2})\cdot aspect}& 0& 0& 0\\ 0& \frac{1}{\tan (\frac{fov}{2})}& 0& 0\\ 0& 0& -\frac{far+near}{far-near}& -\frac{2\cdot far\cdot near}{far-near}\\ 0& 0& -1& 0\end{array}\right]$$

正交投影矩阵：

$$P_{\text{orthographic}}=\left[\begin{array}{cccc}\frac{2}{right-left}& 0& 0& -\frac{right+left}{right-left}\\ 0& \frac{2}{top-bottom}& 0& -\frac{top+bottom}{top-bottom}\\ 0& 0& -\frac{2}{far-near}& -\frac{far+near}{far-near}\\ 0& 0& 0& 1\end{array}\right]$$

七、C语言中的矩阵实现

7.1 矩阵结构体定义

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// 矩阵结构体
typedef struct
 {
    int rows;      // 行数
    int cols;      // 列数
    double** data; // 数据指针
} Matrix;

// 创建矩阵
Matrix* matrix_create(int rows, int cols) 
{
    Matrix* M = (Matrix*)malloc(sizeof(Matrix));
    M->rows = rows;
    M->cols = cols;
    M->data = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        M->data[i] = (double*)calloc(cols, sizeof(double));
    }
    return M;
}

// 释放矩阵
void matrix_free(Matrix* M) 
{
    for (int i = 0; i < M->rows; i++)
    {
        free(M->data[i]);
    }
    free(M->data);
    free(M);
}

7.2 基本运算函数

a. 矩阵加法

// 矩阵加法：C = A + B
Matrix* matrix_add(Matrix* A, Matrix* B) 
{
    if (A->rows != B->rows || A->cols != B->cols) 
    {
        printf("Error: Matrix dimensions mismatch!\n");
        return NULL;
    }
    
    Matrix* C = matrix_create(A->rows, A->cols);
    for (int i = 0; i < A->rows; i++) 
    {
        for (int j = 0; j < A->cols; j++) 
        {
            C->data[i][j] = A->data[i][j] + B->data[i][j];
        }
    }
    return C;
}

b. 矩阵乘法

// 矩阵乘法：C = A * B
Matrix* matrix_multiply(Matrix* A, Matrix* B) 
{
    if (A->cols != B->rows) 
    {
        printf("Error: Matrix dimensions incompatible!\n");
        return NULL;
    }
    
    Matrix* C = matrix_create(A->rows, B->cols);
    for (int i = 0; i < A->rows; i++) 
    {
        for (int j = 0; j < B->cols; j++) 
        {
            C->data[i][j] = 0;
            for (int k = 0; k < A->cols; k++) 
            {
                C->data[i][j] += A->data[i][k] * B->data[k][j];
            }
        }
    }
    return C;
}

c. 矩阵转置

// 矩阵转置
Matrix* matrix_transpose(Matrix* A) 
{
    Matrix* AT = matrix_create(A->cols, A->rows);
    for (int i = 0; i < A->rows; i++) 
    {
        for (int j = 0; j < A->cols; j++) 
        {
            AT->data[j][i] = A->data[i][j];
        }
    }
    return AT;
}

7.3 工程求解函数

a. 高斯消元法求解线性方程组

// 高斯消元法求解 Ax = b
// 返回解向量 x
double* gaussian_elimination(Matrix* A, double* b, int n) 
{
    double** aug = (double**)malloc(n * sizeof(double*));
    for (int i = 0; i < n; i++) 
    {
        aug[i] = (double*)malloc((n + 1) * sizeof(double));
        for (int j = 0; j < n; j++) 
        {
            aug[i][j] = A->data[i][j];
        }
        aug[i][n] = b[i];
    }
    
    // 前向消元
    for (int k = 0; k < n; k++) 
    {
        // 选主元
        int max_row = k;
        for (int i = k + 1; i < n; i++) 
        {
            if (fabs(aug[i][k]) > fabs(aug[max_row][k])) 
            {
                max_row = i;
            }
        }
        // 交换行
        double* temp = aug[k];
        aug[k] = aug[max_row];
        aug[max_row] = temp;
        
        // 消元
        for (int i = k + 1; i < n; i++) 
        {
            double factor = aug[i][k] / aug[k][k];
            for (int j = k; j <= n; j++) 
            {
                aug[i][j] -= factor * aug[k][j];
            }
        }
    }
    
    // 回代求解
    double* x = (double*)malloc(n * sizeof(double));
    for (int i = n - 1; i >= 0; i--) 
    {
        x[i] = aug[i][n];
        for (int j = i + 1; j < n; j++) 
        {
            x[i] -= aug[i][j] * x[j];
        }
        x[i] /= aug[i][i];
    }
    
    for (int i = 0; i < n; i++) 
    {
        free(aug[i]);
    }
    free(aug);
    
    return x;
}

b. 雅可比迭代法（稀疏矩阵求解）

// 雅可比迭代法求解线性方程组
// 适用于对角占优的稀疏矩阵
double* jacobi_iteration(Matrix* A, double* b, int n, double tol, int max_iter) 
{
    double* x = (double*)calloc(n, sizeof(double));
    double* x_new = (double*)malloc(n * sizeof(double));
    
    for (int iter = 0; iter < max_iter; iter++) 
    {
        for (int i = 0; i < n; i++) 
        {
            x_new[i] = b[i];
            for (int j = 0; j < n; j++) 
            {
                if (i != j) {
                    x_new[i] -= A->data[i][j] * x[j];
                }
            }
            x_new[i] /= A->data[i][i];
        }
        
        // 检查收敛
        double diff = 0;
        for (int i = 0; i < n; i++) 
        {
            diff += fabs(x_new[i] - x[i]);
        }
        
        for (int i = 0; i < n; i++) 
        {
            x[i] = x_new[i];
        }
        
        if (diff < tol) 
        {
            printf("Converged after %d iterations\n", iter + 1);
            break;
        }
    }
    
    free(x_new);
    return x;
}

7.4 示例程序

int main() 
{
    // 示例：求解二自由度振动系统
    // [K]{x} = {F}
    
    int n = 2;
    Matrix* K = matrix_create(n, n);
    K->data[0][0] = 1500; K->data[0][1] = -500;
    K->data[1][0] = -500; K->data[1][1] = 500;
    
    double F[] = {1000, 500};
    
    printf("刚度矩阵 K:\n");
    for (int i = 0; i < n; i++) 
    {
        for (int j = 0; j < n; j++) 
        {
            printf("%10.2f ", K->data[i][j]);
        }
        printf("\n");
    }
    
    printf("力向量 F: [%.2f, %.2f]\n", F[0], F[1]);
    
    double* x = gaussian_elimination(K, F, n);
    printf("位移解 x: [%.4f, %.4f] m\n", x[0], x[1]);
    
    // 验证
    double* F_calc = matrix_multiply(K, matrix_create_from_array(x, n, 1))->data[0];
    
    matrix_free(K);
    free(x);
    
    return 0;
}