26、线性代数方程组的求解方法

线性代数方程组的求解方法

在科学和工程领域中，线性代数方程组的求解是一个常见且重要的问题。本文将详细介绍使用 MATLAB 求解线性代数方程组的多种方法，包括矩阵运算、高斯消元法、高斯 - 约旦消元法以及符号数学方法等。

1. 矩阵形式的线性方程组

线性方程组可以写成矩阵形式 (Ax = b)，其中 (A) 是系数矩阵，(x) 是未知变量的列向量，(b) 是方程右侧值的列向量。例如，对于方程组：
[
\begin{cases}
4x_1 - 2x_2 + x_3 = 7 \
x_1 + x_2 + 5x_3 = 10 \
-2x_1 + 3x_2 - x_3 = 2
\end{cases}
]
可以表示为：
[
\begin{bmatrix}
4 & -2 & 1 \
1 & 1 & 5 \
-2 & 3 & -1
\end{bmatrix}
\begin{bmatrix}
x_1 \
x_2 \
x_3
\end{bmatrix}
=
\begin{bmatrix}
7 \
10 \
2
\end{bmatrix}
]
在 MATLAB 中，可以使用两种简单的方法求解：

A = [4 -2 1; 1 1 5; -2 3 -1];
b = [7; 10; 2];
x = inv(A)*b
x = A\b

2. 2x2 方程组的求解

2.1 可视化求解

对于 2x2 方程组，如：
[
\begin{cases}
x_1 + 2x_2 = 2 \
2x_1 + 2x_2 = 6
\end{cases}
]
可以将其转换为直线方程 (y = mx + b) 的形式：
[
\begin{cases}
y = -0.5x + 1 \
y = -x + 3
\end{cases}
]
在 MATLAB 中可以使用以下脚本绘制这些直线：

% Plot a 2 by 2 system as straight lines
x = -2:5;
y1 = -0.5 * x + 1;
y2 = -x + 3;
plot(x,y1,x,y2)
axis([-2 5 -4 6])

两条直线的交点就是方程组的解，即 (x = 4)，(y = -1)，转换回 (x_1) 和 (x_2) 为 (x_1 = 4)，(x_2 = -1)。

2.2 矩阵求解

该方程组的矩阵形式为：
[
\begin{bmatrix}
1 & 2 \
2 & 2
\end{bmatrix}
\begin{bmatrix}
x_1 \
x_2
\end{bmatrix}
=
\begin{bmatrix}
2 \
6
\end{bmatrix}
]
求解方法是 (x = A^{-1}b)，需要先求矩阵 (A) 的逆。对于 2x2 矩阵 (A = \begin{bmatrix}a_{11} & a_{12} \ a_{21} & a_{22}\end{bmatrix})，其行列式 (D) 定义为 (D = a_{11}a_{22} - a_{12}a_{21})，逆矩阵 (A^{-1}) 为：
[
A^{-1} = \frac{1}{D}
\begin{bmatrix}
a_{22} & -a_{12} \
-a_{21} & a_{11}
\end{bmatrix}
]
对于上述矩阵 (A)，(D = 1\times2 - 2\times2 = -2)，则：
[
A^{-1} = \frac{1}{-2}
\begin{bmatrix}
2 & -2 \
-2 & 1
\end{bmatrix}
=
\begin{bmatrix}
-1 & 1 \
1 & -0.5
\end{bmatrix}
]
因此：
[
\begin{bmatrix}
x_1 \
x_2
\end{bmatrix}
=
\begin{bmatrix}
-1 & 1 \
1 & -0.5
\end{bmatrix}
\begin{bmatrix}
2 \
6
\end{bmatrix}
=
\begin{bmatrix}
4 \
-1
\end{bmatrix}
]
在 MATLAB 中可以使用以下代码验证：

a = [1 2; 2 2];
b = [2; 6];
deta = a(1,1)*a(2,2) - a(1,2)*a(2,1);
inva = (1/deta) * [a(2,2) -a(1,2); -a(2,1) a(1,1)];
x = inva * b;

2.3 行列式和逆矩阵的函数使用

MATLAB 提供了内置函数 det 求行列式， inv 求逆矩阵：

det(a)
inv(a)

3. 高斯消元法和高斯 - 约旦消元法

3.1 基本行操作（EROs）

高斯消元法和高斯 - 约旦消元法基于基本行操作（EROs），包括：
1. 替换（Replacement） ：(r_i - sr_j \to r_i)
2. 交换（Interchange） ：(r_i \leftrightarrow r_j)
3. 缩放（Scaling） ：(sr_i \to r_i)

3.2 高斯消元法

高斯消元法的步骤如下：
1. 创建增广矩阵 ([A b])
2. 应用 EROs 将增广矩阵转换为上三角形式（前向消元）
3. 回代求解

例如，对于 2x2 方程组：
[
\begin{cases}
x_1 + 2x_2 = 2 \
2x_1 + 2x_2 = 6
\end{cases}
]
增广矩阵为：
[
\begin{bmatrix}
1 & 2 & 2 \
2 & 2 & 6
\end{bmatrix}
]
应用 EROs 得到上三角形式：
[
\begin{bmatrix}
1 & 2 & 2 \
0 & -2 & 2
\end{bmatrix}
]
回代求解得到 (x_2 = -1)，(x_1 = 4)。

3.3 高斯 - 约旦消元法

高斯 - 约旦消元法与高斯消元法类似，但在得到上三角形式后继续消元，直到得到对角形式。

例如，对于上述 2x2 方程组，继续消元得到：
[
\begin{bmatrix}
1 & 0 & 4 \
0 & -2 & 2
\end{bmatrix}
]
直接得到解 (x_1 = 4)，(x_2 = -1)。

3.4 3x3 方程组示例

对于方程组：
[
\begin{cases}
x_1 + 3x_2 = 1 \
2x_1 + x_2 + 3x_3 = 6 \
4x_1 + 2x_2 + 3x_3 = 3
\end{cases}
]
增广矩阵为：
[
\begin{bmatrix}
1 & 3 & 0 & 1 \
2 & 1 & 3 & 6 \
4 & 2 & 3 & 3
\end{bmatrix}
]
前向消元得到：
[
\begin{bmatrix}
1 & 3 & 0 & 1 \
0 & -5 & 3 & 4 \
0 & 0 & -3 & -9
\end{bmatrix}
]
高斯消元法回代求解得到 (x_3 = 3)，(x_2 = 1)，(x_1 = -2)；高斯 - 约旦消元法继续消元得到对角形式求解。

在 MATLAB 中可以使用以下代码实现：

a = [1 3 0; 2 1 3; 4 2 3];
b = [1 6 3]';
ab = [a b];
ab(2,:) = ab(2,:) - 2*ab(1,:);
ab(3,:) = ab(3,:) - 4 * ab(1,:);
ab(3,:) = ab(3,:) - 2 * ab(2,:);
ab(2,:) = ab(2,:) + ab(3,:);
ab(1,:) = ab(1,:) + 3/5*ab(2,:);

4. 简化行阶梯形（Reduced Row Echelon Form）

简化行阶梯形是高斯 - 约旦消元法的进一步扩展，使得主对角线元素为 1。MATLAB 提供了内置函数 rref 实现：

a = [1 3 0; 2 1 3; 4 2 3];
b = [1 6 3]';
ab = [a b];
rref(ab)

5. 通过增广矩阵求矩阵逆

对于大于 2x2 的方程组，可以通过增广矩阵求矩阵逆，步骤如下：
1. 增广矩阵 ([A I])
2. 化简为 ([I X])，(X) 即为 (A^{-1})

在 MATLAB 中可以使用以下代码实现：

a = [1 3 0; 2 1 3; 4 2 3];
rref([a eye(size(a))])

6. 符号数学

6.1 符号变量和表达式

MATLAB 提供了 sym 类型用于符号变量和表达式，例如：

a = sym('a');
a + a
b = sym('x^2');
c = sym('x^4');
c/b
b^3
c*b
b + sym('4*x^2')

6.2 简化函数

MATLAB 提供了多个简化函数，如 simplify 、 collect 、 expand 、 factor 等：

x = sym('x');
myexpr = cos(x)^2 + sin(x)^2;
simplify(myexpr)
collect(x^2 + 4*x^3 + 3*x^2)
expand((x+2)*(x-1))
factor(ans)

6.3 表达式显示和绘图

pretty 函数可以以更美观的形式显示表达式， ezplot 函数可以绘制符号表达式的图形：

b = sym('x^2');
pretty(b)
ezplot('x^3 + 3*x^2 - 2')
ezplot('cos(x)',[0 pi])

6.4 方程求解

solve 函数可以求解方程，返回符号解，可使用 double 转换为数值解：

x = sym('x');
solve('2*x^2 + x = 6')
double(ans)
solve('3*x^2 + x')
solve('a*x^2 + b*x')
solve('a*x^2 + b*x','b')

7. 常见陷阱和编程风格指南

7.1 常见陷阱

• 混淆矩阵乘法和数组乘法。数组乘法使用 .* ，矩阵乘法使用 * 。
• 忘记在增广矩阵时，两个矩阵的行数必须相同。

7.2 编程风格指南

在处理符号表达式时，通常将所有变量都设为符号变量更方便。

8. 总结

本文介绍了使用 MATLAB 求解线性代数方程组的多种方法，包括矩阵运算、高斯消元法、高斯 - 约旦消元法、简化行阶梯形和符号数学方法。每种方法都有其适用场景，根据具体问题选择合适的方法可以提高求解效率。

8.1 方法总结

方法	适用场景	步骤
矩阵求逆法	适用于系数矩阵可逆的情况	1. 求系数矩阵的逆 2. 矩阵乘法求解
高斯消元法	适用于各种规模的方程组	1. 创建增广矩阵 2. 前向消元得到上三角形式 3. 回代求解
高斯 - 约旦消元法	适用于各种规模的方程组	1. 创建增广矩阵 2. 前向消元得到上三角形式 3. 继续消元得到对角形式求解
简化行阶梯形法	适用于各种规模的方程组	使用 rref 函数直接求解
符号数学法	适用于需要精确解或处理符号表达式的情况	使用 solve 函数求解

8.2 流程图

graph TD;
    A[开始] --> B{选择求解方法};
    B -->|矩阵求逆法| C[求系数矩阵逆];
    C --> D[矩阵乘法求解];
    B -->|高斯消元法| E[创建增广矩阵];
    E --> F[前向消元];
    F --> G[回代求解];
    B -->|高斯 - 约旦消元法| E;
    E --> F;
    F --> H[继续消元得到对角形式];
    H --> I[求解];
    B -->|简化行阶梯形法| J[使用 rref 函数];
    J --> K[得到解];
    B -->|符号数学法| L[使用 solve 函数];
    L --> M[得到符号解];
    M --> N[可转换为数值解];
    D --> O[结束];
    G --> O;
    I --> O;
    K --> O;
    N --> O;

通过以上方法和工具，我们可以在 MATLAB 中高效地求解各种线性代数方程组。

9. 练习题及解答

9.1 矩阵相关练习

9.1.1 矩阵性质判断

对于矩阵 (A = \begin{bmatrix}1 & 4 & 3 \ 2 & 2 & 1 \ 3 & 1 & 5\end{bmatrix})，(B = \begin{bmatrix}6 & 3 & 6 \ 0 & 3 & 2 \ 5 & 4 & 1\end{bmatrix})，(C = \begin{bmatrix}2 & 5 & 4 \ 1 & 2 & 3 \ 6 & 3 & 6\end{bmatrix})：
- 对称矩阵判断 ：对称矩阵满足 (a_{ij}=a_{ji})。经检查，没有对称矩阵。
- 矩阵迹计算 ：矩阵的迹是主对角线元素之和。(A) 的迹为 (1 + 2 + 5 = 8)；(B) 的迹为 (6 + 3 + 1 = 10)；(C) 的迹为 (2 + 2 + 6 = 10)。
- 数乘矩阵 ：(3 A=\begin{bmatrix}3 & 12 & 9 \ 6 & 6 & 3 \ 9 & 3 & 15\end{bmatrix})
- 矩阵乘法 ：(A C=\begin{bmatrix}1\times2 + 4\times1+3\times6 & 1\times5 + 4\times2+3\times3 & 1\times4 + 4\times3+3\times6 \ 2\times2 + 2\times1+1\times6 & 2\times5 + 2\times2+1\times3 & 2\times4 + 2\times3+1\times6 \ 3\times2 + 1\times1+5\times6 & 3\times5 + 1\times2+5\times3 & 3\times4 + 1\times3+5\times6\end{bmatrix}=\begin{bmatrix}24 & 22 & 34 \ 12 & 17 & 20 \ 37 & 22 & 45\end{bmatrix})
还可以进行的矩阵乘法有 (C*A)。

9.1.2 矩阵运算练习

对于矩阵 (A = \begin{bmatrix}3 & 2 & 1 \ 0 & 5 & 2 \ 1 & 0 & 3\end{bmatrix})，(B = \begin{bmatrix}2 & 1 & 3 \ 1 & 0 & 0 \ 0 & 0 & 1\end{bmatrix})，(C = \begin{bmatrix}0 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1\end{bmatrix})：
- (A * B=\begin{bmatrix}3\times2 + 2\times1+1\times0 & 3\times1 + 2\times0+1\times0 & 3\times3 + 2\times0+1\times1 \ 0\times2 + 5\times1+2\times0 & 0\times1 + 5\times0+2\times0 & 0\times3 + 5\times0+2\times1 \ 1\times2 + 0\times1+3\times0 & 1\times1 + 0\times0+3\times0 & 1\times3 + 0\times0+3\times1\end{bmatrix}=\begin{bmatrix}8 & 3 & 10 \ 5 & 0 & 2 \ 2 & 1 & 6\end{bmatrix})
- (B * A) 也可以进行计算，过程类似。
- (I + A)，这里 (I) 是单位矩阵，(I=\begin{bmatrix}1 & 0 & 0 \ 0 & 1 & 0 \ 0 & 0 & 1\end{bmatrix})，则 (I + A=\begin{bmatrix}4 & 2 & 1 \ 0 & 6 & 2 \ 1 & 0 & 4\end{bmatrix})
- (A.*I=\begin{bmatrix}3\times1 & 2\times0 & 1\times0 \ 0\times0 & 5\times1 & 2\times0 \ 1\times0 & 0\times0 & 3\times1\end{bmatrix}=\begin{bmatrix}3 & 0 & 0 \ 0 & 5 & 0 \ 0 & 0 & 3\end{bmatrix})
- (trace(A)=3 + 5 + 3 = 11)

9.2 自编函数练习

9.2.1 判断方阵函数

function result = issquare(arr)
    [rows, cols] = size(arr);
    if rows == cols
        result = 1;
    else
        result = 0;
    end
end

9.2.2 提取主对角线函数

function diag_vec = mydiag(arr)
    if issquare(arr)
        [n, ~] = size(arr);
        diag_vec = zeros(1, n);
        for i = 1:n
            diag_vec(i) = arr(i, i);
        end
    else
        diag_vec = [];
    end
end

9.2.3 返回对角线和迹的函数

function [diag_vec, tr] = mydiag_trace(arr)
    [n, ~] = size(arr);
    diag_vec = zeros(1, n);
    tr = 0;
    for i = 1:n
        diag_vec(i) = arr(i, i);
        tr = tr + arr(i, i);
    end
end

9.2.4 生成对角矩阵函数

function mat = randdiag(n, low, high)
    mat = zeros(n);
    for i = 1:n
        mat(i, i) = randi([low, high]);
    end
end

9.2.5 生成单位矩阵函数

function mat = myeye(n)
    mat = zeros(n);
    for i = 1:n
        mat(i, i) = 1;
    end
end

9.2.6 生成上三角矩阵函数

function mat = myupp(n)
    mat = randi([1, 10], n);
    for i = 1:n
        for j = 1:i - 1
            mat(i, j) = 0;
        end
    end
end

9.2.7 判断对角矩阵函数

function result = isdiag(mat)
    [n, m] = size(mat);
    if n ~= m
        result = 0;
        return;
    end
    for i = 1:n
        for j = 1:n
            if i ~= j && mat(i, j) ~= 0
                result = 0;
                return;
            end
        end
    end
    result = 1;
end

9.2.8 矩阵相减函数

function result = mymatsub(mat1, mat2)
    [r1, c1] = size(mat1);
    [r2, c2] = size(mat2);
    if r1 == r2 && c1 == c2
        result = zeros(r1, c1);
        for i = 1:r1
            for j = 1:c1
                result(i, j) = mat1(i, j) - mat2(i, j);
            end
        end
    else
        result = [];
    end
end

9.2.9 矩阵转置函数

function mat_T = mytranspose(mat)
    [rows, cols] = size(mat);
    mat_T = zeros(cols, rows);
    for i = 1:rows
        for j = 1:cols
            mat_T(j, i) = mat(i, j);
        end
    end
end

9.2.10 生成全 5 矩阵函数

function mat = fives(rows, cols)
    mat = 5 * ones(rows, cols);
end

9.2.11 生成帕斯卡矩阵函数

function mat = mypascal(n)
    mat = zeros(n);
    for i = 1:n
        mat(i, 1) = 1;
        mat(i, i) = 1;
    end
    for i = 3:n
        for j = 2:i - 1
            mat(i, j) = mat(i - 1, j - 1) + mat(i - 1, j);
        end
    end
end

9.3 方程组求解练习

9.3.1 方程组矩阵形式及求解

对于方程组 (\begin{cases}4x_1 - x_2 + 3x_4 = 10 \ -2x_1 + 3x_2 + x_3 - 5x_4 = -3 \ x_1 + x_2 - x_3 + 2x_4 = 2 \ 3x_1 + 2x_2 - 4x_3 = 4\end{cases})
矩阵形式为 (A = \begin{bmatrix}4 & -1 & 0 & 3 \ -2 & 3 & 1 & -5 \ 1 & 1 & -1 & 2 \ 3 & 2 & -4 & 0\end{bmatrix})，(b = \begin{bmatrix}10 \ -3 \ 2 \ 4\end{bmatrix})
在 MATLAB 中可以使用 x = A\b 求解。

9.3.2 2x2 方程组特殊情况

对于方程组 (\begin{cases}-3x_1 + x_2 = -4 \ -6x_1 + 2x_2 = 4\end{cases})
- 直线绘图 ：将方程改写为 (y = 3x - 4) 和 (y = 3x + 2)，两条直线平行，无交点。
- 求解分析 ：由第一个方程得 (x_2 = 3x_1 - 4)，代入第二个方程 (-6x_1 + 2(3x_1 - 4)=4)，化简得 (-8 = 4)，矛盾，无解。
- 行列式计算 ：(A=\begin{bmatrix}-3 & 1 \ -6 & 2\end{bmatrix})，(D=-3\times2 - 1\times(-6)=0)，行列式为 0，方程组无解。

9.3.3 2x2 方程组正常求解

对于方程组 (\begin{cases}-3x_1 + x_2 = 2 \ -6x_1 + 2x_2 = 4\end{cases})
- 直线绘图 ：改写为 (y = 3x + 2) 和 (y = 3x + 2)，两条直线重合，有无数解。
- 求解分析 ：由第一个方程得 (x_2 = 3x_1 + 2)，代入第二个方程恒成立。
- 行列式计算 ：(A=\begin{bmatrix}-3 & 1 \ -6 & 2\end{bmatrix})，(D=-3\times2 - 1\times(-6)=0)，行列式为 0，方程组有无数解。

9.3.4 2x2 方程组多种方法求解

对于方程组 (\begin{cases}3x_1 + x_2 = 2 \ 2x_1 = 4\end{cases})
- 矩阵形式 ：(A=\begin{bmatrix}3 & 1 \ 2 & 0\end{bmatrix})，(b=\begin{bmatrix}2 \ 4\end{bmatrix})
- 行列式计算 ：(D = 3\times0 - 1\times2=-2)
- 逆矩阵求解 ：(A^{-1}=\frac{1}{-2}\begin{bmatrix}0 & -1 \ -2 & 3\end{bmatrix}=\begin{bmatrix}0 & 0.5 \ 1 & -1.5\end{bmatrix})，(x = A^{-1}b=\begin{bmatrix}0 & 0.5 \ 1 & -1.5\end{bmatrix}\begin{bmatrix}2 \ 4\end{bmatrix}=\begin{bmatrix}2 \ -4\end{bmatrix})
- 高斯消元法 ：增广矩阵 (\begin{bmatrix}3 & 1 & 2 \ 2 & 0 & 4\end{bmatrix})，经过行变换得到上三角形式求解。
- 高斯 - 约旦消元法 ：继续消元得到对角形式求解。

9.3.5 3x3 方程组求解

对于方程组 (\begin{cases}2x_1 + 2x_2 + x_3 = 2 \ x_2 + 2x_3 = 1 \ x_1 + x_2 + 3x_3 = 3\end{cases})
- 增广矩阵 ：(\begin{bmatrix}2 & 2 & 1 & 2 \ 0 & 1 & 2 & 1 \ 1 & 1 & 3 & 3\end{bmatrix})
- 高斯消元法 ：经过前向消元得到上三角形式，回代求解。
- 高斯 - 约旦消元法 ：继续消元得到对角形式求解。
- MATLAB 求解 ：

A = [2 2 1; 0 1 2; 1 1 3];
b = [2; 1; 3];
x = A\b;
det_A = det(A);
inv_A = inv(A);

9.4 实际应用问题求解

9.4.1 胰岛素清除问题

已知 (C = C_0 e^{-\frac{30t}{m}})，(m = 65)，(C_0 = 90)，(C = 10)

syms t;
eq = 10 == 90*exp(-30*t/65);
sol = solve(eq, t);
t_val = double(sol);

9.4.2 电路电压问题

对于方程组 (\begin{cases}2(V_a - V_b)+5(V_a - V_c)-e^{-t}=0 \ 2(V_b - V_a)+2V_b + 3(V_b - V_c)=0 \ V_c = 2\sin(t)\end{cases})

syms Va Vb Vc t;
eq1 = 2*(Va - Vb)+5*(Va - Vc)-exp(-t)==0;
eq2 = 2*(Vb - Va)+2*Vb + 3*(Vb - Vc)==0;
eq3 = Vc == 2*sin(t);
sol = solve([eq1, eq2, eq3], [Va, Vb, Vc]);
Va_sol = sol.Va;
Vb_sol = sol.Vb;
Vc_sol = sol.Vc;

9.4.3 细菌繁殖问题

已知 (\log(N)=\log(N_0)+\frac{t}{T}\log(2))，(N_0 = 10^2)，(N = 10^8)，(t = 8)

syms T;
eq = log10(10^8)==log10(10^2)+(8/T)*log10(2);
sol = solve(eq, T);
T_val = double(sol);

9.5 矩阵性质判断函数

9.5.1 判断反对称矩阵函数

function result = isskewsymmetric(mat)
    [n, m] = size(mat);
    if n ~= m
        result = 0;
        return;
    end
    for i = 1:n
        for j = 1:n
            if mat(i, j) ~= -mat(j, i)
                result = 0;
                return;
            end
        end
    end
    result = 1;
end

10. 总结与回顾

10.1 方法对比

方法	优点	缺点	适用场景
矩阵求逆法	计算简单，适用于系数矩阵可逆的情况	当矩阵不可逆时无法使用	系数矩阵可逆且规模较小的方程组
高斯消元法	通用方法，适用于各种规模的方程组	计算量较大	一般线性方程组求解
高斯 - 约旦消元法	可直接得到解，无需回代	计算量比高斯消元法更大	需要直接得到解的情况
简化行阶梯形法	代码简单，可直接使用 MATLAB 内置函数	对于大规模矩阵计算效率可能不高	快速求解方程组
符号数学法	可得到精确的符号解	计算速度慢，可能无法得到数值解	需要精确解或处理符号表达式的情况

10.2 学习建议

• 对于初学者，先掌握基本的矩阵运算和方程组求解方法，如矩阵求逆法和高斯消元法。
• 随着学习的深入，了解不同方法的优缺点，根据具体问题选择合适的方法。
• 多做练习题，通过实践加深对知识的理解和掌握。

10.3 未来展望

线性代数方程组的求解在科学计算、工程技术等领域有着广泛的应用。随着计算机技术的发展，求解方法也在不断优化和改进。未来，我们可以期待更高效、更精确的求解算法出现，同时也可以将这些方法应用到更多的实际问题中。

graph LR
    A[开始学习] --> B[掌握基本矩阵运算]
    B --> C[学习方程组求解方法]
    C --> D[了解方法优缺点]
    D --> E[根据问题选方法]
    E --> F[多做练习题]
    F --> G[应用到实际问题]
    G --> H[关注新算法发展]

通过对线性代数方程组求解方法的学习和实践，我们可以更好地解决各种实际问题，提高自己的科学计算能力。