第13讲：综合应用1
求解线性方程组

《计算机程序设计》

苏醒

计算机学院计算机研究所编译系统研究室

课前准备

  • 在WSL Ubuntu中建立一个专门的目录用于课上练习，如~/course
$ mkdir ~/course
  • 打开VSCode并连接到WSL，打开此目录并新建文件，如~/course/main.cpp
  • 编辑main.cpp，随堂编程练习的代码请直接在此文件中编辑
main.cpp
#include <iostream>
#include <iomanip>
using namespace std;

int main() {
  // ============= begin =============

  // =============  end  =============
  return 0;
}

学习目标

  • 学习内容
    • 多维数组索引的线性化
    • 高斯消元法求解线性方程组
      • 经典消元法
      • LU分解法

学习目标

  • 掌握多维数组索引的线性化寻址方法
  • 应用线性化寻址实现线性方程组求解方法

一、多维数组索引的线性化

1.1 二维数组作为函数参数（回顾）

arrayparam.cpp
#include <iostream>
void printMatrix(double matrix[][7], int rows, int cols) {
  for (int i = 0; i < rows; ++i) {
    for (int j = 0; j < cols; ++j)
      std::cout << matrix[i][j] << " ";
    std::cout << std::endl;
  }
}
  • 多维数组作为函数参数，除最高维外其他维度的长度必须为静态常量

明察秋毫

  • 通过多维索引计算元素地址，依赖于维度长度。
    • 对于MxN矩阵matrix，matrix[i][j]将等价于*(matrix + i * N + j)
  • 若维度长度必须是静态常量，则无法支持动态维度长度！

1.2 多维索引的线性化

解决办法

  • 参照多维索引寻址规则，显式计算多维索引相对于首地址的偏移量，以支持动态维度长度
    • 令指针p指向MxN矩阵matrix，matrix[i][j]将等价于*(p + i * N + j)
linearize.cpp
#include <iostream>
void printMatrix(double *matrix, int rows, int cols) {
  for (int i = 0; i < rows; ++i) {
    for (int j = 0; j < cols; ++j)
      std::cout << matrix[i * cols + j] << " ";
    std::cout << std::endl;
  }
}
int main() {
  int m = 0, n = 0;
  std::cin >> m >> n;
  double *matrix = new double[m * n];
  for (int i = 0; i < m * n; ++i)
    matrix[i] = i + 1;
  printMatrix(matrix, m, n);
  delete[] matrix;
  return 0;
}
$ g++ -o code/linearize code/linearize.cpp
$ echo 3 4 | ./code/linearize
1 2 3 4 
5 6 7 8 
9 10 11 12 
$ echo 4 3 | ./code/linearize
1 2 3 
4 5 6 
7 8 9 
10 11 12

二、经典高斯消元法

2.1 高斯消元法

  • 高斯消元法是一种求解线性方程组的经典方法，其基本步骤如下：
    • 将方程组表示为增广矩阵
    • 通过初等行变换将系数矩阵化为还原矩阵
    • 从最后一行开始，逐行向上回代求解未知数
\[\begin{cases} x_1 + 2x_2 + 3x_3 & = & 14 \\ 5x_1 + 6x_2 + 7x_3 & = & 38 \\ 9x_1 + 10x_2 + 10x_3 & = & 59 \end{cases}\]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & 0 & 1 & 3 \end{array} \right] \]

\[x = 1, y = 2, z = 3\]

2.2 增广矩阵

  • n元一次方程组，可由为其系数\(A\)与右端项\(b\)合并得到的增广矩阵\([A|b]\)表示

\[Ax=b, A \in \mathbb{R}^{n \times n}, b \in \mathbb{R}^n\]

\[\begin{cases} x_1 + 2x_2 + 3x_3 & = & 14 \\ 5x_1 + 6x_2 + 7x_3 & = & 38 \\ 9x_1 + 10x_2 + 10x_3 & = & 59 \end{cases}\]

\[\overset{\text{写作增广矩阵}} \Longrightarrow\]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

实现方式

  • 增广矩阵以行优先方式在内存中连续存储（n行n+1列），前n列存储系数矩阵A，最后一列存储右端项b
  • 可由两个变量表示：（1）double *Ab：增广矩阵首地址，（2）int n：未知数个数（增广矩阵行数）

2.3 初等行变换

  • 初等行变换
    • 行交换：交换两行
    • 行倍乘：将一行乘以非零常数
    • 行累加：将一行加上另一行的倍数

实现方式

  • 每种初等变换封装为一个函数
函数签名 描述
void rswap(double *Ab, int n, int i, int j) 交换第i、j行
void rscale(double *Ab, int n, int i, double k) 将第i行乘以k
void radd(double *Ab, int n, int i, int j, double k) 将i行加上j行的k倍

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[\overset{\text{交换第2、3行}} \Longrightarrow\]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 9 & 10 & 10 & 59 \\ 5 & 6 & 7 & 38 \end{array} \right] \]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[\overset{\text{第1行倍乘-5}} \Longrightarrow\]

\[ \left[ \begin{array}{ccc:c} -5 & -10 & -15 & -70 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[\overset{\text{第2行加上第1行的-5倍}} \Longrightarrow\]

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & -4 & -8 & -32 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

明察秋毫

初等行变换不改变线性方程组的解（可能打乱解向量中各分量的顺序）

2.3 初等行变换

初等行变换
void rswap(double *Ab, int n, int i, int j) {
  for (int k = 0; k < n + 1; ++k) {
    double temp = Ab[i * (n + 1) + k];
    Ab[i * (n + 1) + k] = Ab[j * (n + 1) + k];
    Ab[j * (n + 1) + k] = temp;
  }
}

void rscale(double *Ab, int n, int i, double k) {
    // ================= begin =================

    // ================= begin =================
}
void radd(double *Ab, int n, int i, int j, double k) {
    // ================= begin =================

    // =================  end  =================
}

2.4 正向消元

  • 不断使用初等行变换进行消元，将增广矩阵转化为还原矩阵
    • 每一行的首个非零元素为1，称为主元，主元所在列称为主元列
    • 主元列的主元必须在其上方的主元列的右侧
    • 主元列下方的所有元素均为0，主元列上方的元素不受限制

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & 0 & 1 & 3 \end{array} \right] \]

明察秋毫

还原矩阵由上三角系数矩阵与右端项组成

2.4 正向消元

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 5 & 6 & 7 & 38 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[ \overset{\text{add(2,1,-5)}}\Longrightarrow \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & -4 & -8 & -32 \\ 9 & 10 & 10 & 59 \end{array} \right] \]

\[ \overset{\text{add(3,1,-9)}}\Longrightarrow \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & -4 & -8 & -32 \\ 0 & -8 & -17 & -67 \end{array} \right] \]

\[ \overset{\text{scale(2,-0.25)}}\Longrightarrow \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & -8 & -17 & -67 \end{array} \right] \]

\[ \overset{\text{add(3,2,8)}}\Longrightarrow \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & 0 & -1 & -3 \end{array} \right] \]

\[ \overset{\text{scale(3,-1)}}\Longrightarrow \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & 0 & 1 & 3 \end{array} \right] \]

2.4 正向消元

  • 假设在消元过程中，主元均不为0，即无需进行交换
通过初等行变换将增广矩阵转化为还原矩阵
void rswap(double *Ab, int n, int i, int j);
void rscale(double *Ab, int n, int i, double k);
void radd(double *Ab, int n, int i, int j, double k);

void triangularize(double *Ab, int n) {
  for (int i = 0; i < n; ++i) {
    // ================= begin =================
    // 倍乘第i行主元使之为1

    // 消除第i行主元下方的元素

    // =================  end  =================
  }
}

实现方式

还原矩阵存储在原始增广矩阵的内存中（前n列，系数部分），不额外分配内存

2.5 回代求解

  • 还原矩阵对应一个三角形线性方程组，该方程组可通过回代方式求解
    • 从最后一行开始，逐行向上代入已求解的未知数

\[ \left[ \begin{array}{ccc:c} 1 & 2 & 3 & 14 \\ 0 & 1 & 2 & 8 \\ 0 & 0 & 1 & 3 \end{array} \right] \]

  • 解第3行，得\(x_3 = 3\)
  • 将\(x_3\)代入第2行，得\(x_2 = 2\)
  • 将\(x_2\)、\(x_3\)代入第1行，得\(x_1 = 1\)

2.5 回代求解

回代求解
void trsvu(double *Ab, int n) {
  for (int i = n - 1; i >= 0; --i) {
    // ================= begin =================

    // =================  end  =================
  }
}

实现方式

解向量存储在原始增广矩阵的内存中（最后1列，右端项部分），不额外分配内存

三、LU分解法求解线性方程组

3.1 初等行变换的矩阵表示

  • 初等行变换可以表示为矩阵乘法操作
  • 行交换：交换第i行与第j行相当于左乘一个行交换矩阵\(E_{i,j}\)

\[E_{1,3} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}, \]

\[ A = \begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{bmatrix}\Rightarrow \]

\[ E_{1,3}A = \begin{bmatrix} 9 & 10 & 11 & 12 \\ 5 & 6 & 7 & 8 \\ 1 & 2 & 3 & 4 \\ 13 & 14 & 15 & 16 \end{bmatrix} \]

明察秋毫

\(E_{i,j}\)由单位矩阵\(E\)第\(i\)列与第\(j\)列互换得到

3.1 初等行变换的矩阵表示

  • 初等行变换可以表示为矩阵乘法操作
  • 行倍乘：第i行倍乘\(\alpha\)，相当于左乘一个行倍乘矩阵\(E_i^\alpha\)

\[E_2^{-2} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & -2 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

\[ A = \begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{bmatrix} \]

\[ E_{2}^{-2}A = \begin{bmatrix} 1 & 2 & 3 & 4 \\ -10 & -12 & -14 & -16 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{bmatrix} \]

明察秋毫

\(E_i^\alpha\)由单位矩阵\(E\)第\(i\)列倍乘\(\alpha\)得到

3.1 初等行变换的矩阵表示

  • 初等行变换可以表示为矩阵乘法操作
  • 行累加：第i行加上第j行的\(\alpha\)倍，相当于左乘一个行累加矩阵\(E_{i,j}^\alpha\)

\[E_{3,1}^{-2} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ -2 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \]

\[ A = \begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \end{bmatrix} \]

\[ E_{3,1}^{-2}A = \begin{bmatrix} 1 & 2 & 3 & 4 \\ 5 & 6 & 7 & 8 \\ 7 & 6 & 5 & 4 \\ 13 & 14 & 15 & 16 \end{bmatrix} \]

明察秋毫

\(E_{i,j}^\alpha\)由单位矩阵\(E\)第\(i\)行加上第\(j\)行的\(\alpha\)倍得到

3.2 正向消元的矩阵表示

正向消元由一系列初等行变换操作组成

\[L_K \cdots L_2 L_1 [A|b] = [U|b']\]

其中，\(L_{k}\)（\(k=1,2,\cdots,K\)）为初等行变换，\(U\)为上三角矩阵

由\(L_K \cdots L_2 L_1 A = U\)可得

\[A = L_1^{-1} L_2^{-1} \cdots L_K^{-1} U = LU\]

定理

定理 1 矩阵\(L=L_1^{-1} L_2^{-1} \cdots L_K^{-1}\)是一个下三角矩阵

3.3 LU分解法求解线性方程组

原线性方程组\(Ax=b\)可写作

\[LUx=b\]

令\(y=Ux\)，原问题可分解为两个三角线性方程组求解问题

\[Ly=b, Ux=y\]

明察秋毫

LU分解法求解线性方程组的步骤：

  1. 系数矩阵分解\(A=LU\)，其中\(L\)为对角线元素为1的下三角矩阵，\(U\)为上三角矩阵
  2. 通过正代求解\(Ly=b\)，得到中间变量\(y\)
  3. 通过回代求解\(Ux=y\)，得到最终解\(x\)

3.4 LU分解算法

令

\[L= \begin{bmatrix} L_{11} & 0\\ L_{21} & L_{22} \end{bmatrix} \]

\[ U= \begin{bmatrix} U_{11} & U_{12}\\ 0 & U_{22} \end{bmatrix} \]

\[ A= \begin{bmatrix} A_{11} & A_{12}\\ A_{21} & A_{22} \end{bmatrix} \]

由\(A=LU\)可得

\[ \begin{aligned} A_{11} &= L_{11}U_{11} \\ A_{12} &= L_{11}U_{12} \\ A_{21} &= L_{21}U_{11} \\ A_{22} &= L_{21}U_{12} + L_{22}U_{22} \end{aligned} \]

3.4 LU分解算法

设置分块大小使\(L_{11} \in \mathbb{R}^{1 \times 1}\)

\[ \begin{aligned} A_{11} &= L_{11}U_{11} \\ A_{12} &= L_{11}U_{12} \\ A_{21} &= L_{21}U_{11} \\ A_{22} &= L_{21}U_{12} + L_{22}U_{22} \end{aligned} \]

\[\Rightarrow\]

\[ \begin{aligned} l_{11} = 1 &, u_{11} = a_{11} \\ U_{12} &= A_{12} \\ L_{21} &= A_{21} / u_{11} \\ L_{22}U_{22} &= A_{22} - L_{21}U_{12} \end{aligned} \]

明察秋毫

  • 沿对角线从左上角向右下角进行计算
    • 对角线元素\(a_{ii}\)分解为\(l_{ii}=1, u_{ii} =a_{ii}\)
    • \(a_{ii}\)右侧的元素\(a_{ij}\)即为\(u_{ij}\)
    • \(a_{ii}\)下方的元素\(a_{ji}\)除以\(u_{ii}\)即为\(L_{ji}\)

3.4 LU分解算法

\[ \left[ \begin{array}{c:cc} 2 & 3 & 1 \\ \hdashline 4 & 7 & 4 \\ -2 & -2 & 3 \end{array} \right] \]

\[\Rightarrow \left[ \begin{array}{c:cc} 2 & 3 & 1 \\ \hdashline 2 & 7 & 4 \\ -1 & -2 & 3 \end{array} \right] \]

\[\Rightarrow \left[ \begin{array}{c:cc} 2 & 3 & 1 \\ \hdashline 2 & 1 & 2 \\ -1 & 1 & 4 \end{array} \right] \]

\[\Rightarrow \left[ \begin{array}{cc:c} 2 & 3 & 1 \\ 2 & 1 & 2 \\ \hdashline -1 & 1 & 4 \end{array} \right] \]

\[\Rightarrow \left[ \begin{array}{cc:c} 2 & 3 & 1 \\ 2 & 1 & 2 \\ \hdashline -1 & 1 & 2 \end{array} \right] \]

\[\Rightarrow L = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 1 & 1 \end{array} \right], U = \left[ \begin{array}{ccc} 2 & 3 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{array} \right] \]

3.4 LU分解算法

LU分解算法
void lu(double *Ab, int n) {
  for (int i = 0; i < n; ++i) {
    // L11=1，U11=A11，无需计算
    // U12=A12，无需计算
    // 计算L21

    // 计算L22U22

  }
}

实现方式

  • 解向量存储在原始增广矩阵的内存中（最后1列，右端项部分），不额外额外分配内存
  • L矩阵与U矩阵存储在原始增广矩阵的内存中（前n列，系数部分），不额外额外分配内存
    • L矩阵存储在增广矩阵的下三角部分，U矩阵存储在增广矩阵的上三角部分
    • L矩阵的对角线元素均为1，无需存储

总结

本节内容

---
config:
  look: handDrawn
  themeVariables:
    fontSize: 20px
---
mindmap
线性方程组求解
  多维数组索引的线性化
    方式同多维数组索引寻址
  经典高斯消元法
    增广矩阵
    初等行变换
    正向消元
    还原矩阵
    回代求解
  LU分解算法
    初等行变换的矩阵表示
    LU分解法原理
    LU分解算法与实现

学习目标

  • 掌握多维数组索引的线性化寻址方法
  • 应用线性化寻址实现线性方程组求解方法

课后作业

  • 实训（截止时间2025.04.07）
    • 求解线性方程组
  • 复习（准备第二次单元测试）
    • 数组、指针、字符串

计算机程序设计