15、矩阵的更多应用：从种群动态到稀疏矩阵

最新推荐文章于 2026-01-05 12:14:24 发布

原创最新推荐文章于 2026-01-05 12:14:24 发布 · 34 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#矩阵 #种群动态 #Leslie矩阵

MATLAB入门实战精讲专栏收录该内容

17 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

矩阵的更多应用：从种群动态到稀疏矩阵

1. 矩阵在种群动态中的应用

1.1 种群建模思路

在种群动态建模中，以兔子种群为例，可将其按年龄划分为不同的年龄组。设 (X_i) 为第 (i) 个年龄组的种群数量，生存因子 (P_i) 表示第 (i) 个年龄组存活到第 (i + 1) 个年龄组的比例，生育力 (F_i) 表示第 (i) 个年龄组在一个时间间隔内每个成员平均生产的新生个体数量（在生物建模中通常只考虑雌性）。

1.2 兔子种群模型示例

假设兔子种群分为三个年龄组：年轻组、中年组和老年组，时间单位为一个月。已知年轻组不能繁殖（(F_1 = 0)），中年组生育力 (F_2 = 9)，老年组生育力 (F_3 = 12)；年轻组到中年组的生存概率 (P_1 = 1/3)，中年组到老年组的生存概率 (P_2 = 0.5)，且假设老年组兔子在月底全部死亡。

用 (t) 表示当前月份，(t + 1) 表示下一个月，可得到种群更新方案：
- (X_1(t + 1) = F_2X_2(t) + F_3X_3(t))
- (X_2(t + 1) = P_1X_1(t))
- (X_3(t + 1) = P_2X_2(t))

定义种群向量 (X(t) = [X_1(t), X_2(t), X_3(t)]^T)，上述方程可改写为矩阵方程 (X(t + 1) = LX(t))，其中 (L) 为 Leslie 矩阵：
[
L =
\begin{bmatrix}
0 & 9 & 12 \
1/3 & 0 & 0 \
0 & 1/2 & 0
\end{bmatrix}
]

1.3 MATLAB 实现

以下是用 MATLAB 实现的代码：

% Leslie matrix population model
n = 3;
L = zeros(n); % all elements set to zero
L(1,2) = 9;
L(1,3) = 12;
L(2,1) = 1/3;
L(3,2) = 0.5;
x = [0 0 1]'; % remember x must be a column vector!
for t = 1:24
    x = L * x;
    disp( [t x' sum(x)] ) % x’ is a row
end

运行该代码可得到 24 个月内兔子种群的变化情况，部分输出如下：
| Month | Young | Middle | Old | Total |
| — | — | — | — | — |
| 1 | 12 | 0 | 0 | 12 |
| 2 | 0 | 4 | 0 | 4 |
| 3 | 36 | 0 | 2 | 38 |
| 4 | 24 | 12 | 0 | 36 |
| 5 | 108 | 8 | 6 | 122 |
|… |… |… |… |… |
| 22 | 11184720 | 1864164 | 466020 | 13514904 |
| 23 | 22369716 | 3728240 | 932082 | 27030038 |
| 24 | 44739144 | 7456572 | 1864120 | 54059836 |

1.4 种群增长特性

增长因子 ：观察输出可发现，经过几个月后，种群总数每月翻倍，这个因子称为增长因子，它是 Leslie 矩阵的一个特性（即矩阵的主特征值）。在本示例中增长因子为 2，改变 Leslie 矩阵的值，长期增长因子也会改变。
稳定年龄分布 ：三个年龄组的数量会趋于一个极限比例 24:4:1，这个极限比例称为种群的稳定年龄分布，它也是 Leslie 矩阵的一个特性（即属于主特征值的特征向量）。无论初始种群结构如何，给定的 Leslie 矩阵最终都会使种群达到相同的稳定年龄分布，并以相同的增长因子增长。

1.5 绘制种群增长曲线

若要绘制前 15 个月兔子种群总数的增长曲线，可在上述代码的 for 循环中插入 p(t) = sum(x); ，然后运行以下命令：

plot(1:15, p(1:15)), xlabel('months'), ylabel('rabbits')
hold, plot(1:15, p(1:15),'o')

该曲线呈现指数增长的特征。

2. 马尔可夫过程

2.1 马尔可夫链概念

许多过程可以用一系列离散状态来描述，若过程是随机的，在给定时刻处于某个状态有一定概率，且从一个状态转移到另一个状态的概率仅取决于当前状态，而与之前的状态无关，则该过程称为马尔可夫链。

2.2 随机游走示例

假设有一条街道有六个十字路口，一个近视的学生在街道上随机游走，他家在十字路口 1，他喜欢的网吧在十字路口 6。在除家与网吧外的每个十字路口，他向网吧方向移动的概率为 2/3，向家方向移动的概率为 1/3。当他到达家或网吧时，就会“消失”（被吸收）。

学生在随机游走过程中有六种状态，可用六维状态向量 (X(t)) 表示在时刻 (t) 处于各个十字路口的概率，且向量各分量之和为 1。

该马尔可夫过程的转移概率矩阵 (P) 如下：
| | Home | 2 | 3 | 4 | 5 | Cafe |
| — | — | — | — | — | — | — |
| Home | 1 | 1/3 | 0 | 0 | 0 | 0 |
| 2 | 0 | 0 | 1/3 | 0 | 0 | 0 |
| 3 | 0 | 2/3 | 0 | 1/3 | 0 | 0 |
| 4 | 0 | 0 | 2/3 | 0 | 1/3 | 0 |
| 5 | 0 | 0 | 0 | 2/3 | 0 | 0 |
| Cafe | 0 | 0 | 0 | 0 | 2/3 | 1 |

根据转移概率矩阵 (P)，可得到时刻 (t + 1) 处于十字路口 3 的概率为 (X_3(t + 1) = 2/3X_2(t) + 1/3X_4(t))。同样可以用矩阵方程 (X(t + 1) = PX(t)) 来更新状态向量。

2.3 MATLAB 实现

以下是用 MATLAB 实现的代码：

n = 6;
P = zeros(n); % all elements set to zero
for i = 3:6
    P(i,i-1) = 2/3;
    P(i-2,i-1) = 1/3;
end
P(1,1) = 1;
P(6,6) = 1;
x = [0 1 0 0 0 0]'; % remember x must be a column vector!
for t = 1:50
    x = P * x;
    disp( [t x'] )
end

部分输出如下：
| Time | Home | 2 | 3 | 4 | 5 | Cafe |
| — | — | — | — | — | — | — |
| 1 | 0.3333 | 0 | 0.6667 | 0 | 0 | 0 |
| 2 | 0.3333 | 0.2222 | 0 | 0.4444 | 0 | 0 |
| 3 | 0.4074 | 0 | 0.2963 | 0 | 0.2963 | 0 |
| 4 | 0.4074 | 0.0988 | 0 | 0.2963 | 0 | 0.1975 |
| 5 | 0.4403 | 0 | 0.1646 | 0 | 0.1975 | 0.1975 |
| 6 | 0.4403 | 0.0549 | 0 | 0.1756 | 0 | 0.3292 |
| 7 | 0.4586 | 0 | 0.0951 | 0 | 0.1171 | 0.3292 |
| 8 | 0.4586 | 0.0317 | 0 | 0.1024 | 0 | 0.4073 |
|… |… |… |… |… |… |… |
| 20 | 0.4829 | 0.0012 | 0 | 0.0040 | 0 | 0.5119 |
|… |… |… |… |… |… |… |
| 40 | 0.4839 | 0.0000 | 0 | 0.0000 | 0 | 0.5161 |
|… |… |… |… |… |… |… |
| 50 | 0.4839 | 0.0000 | 0 | 0.0000 | 0 | 0.5161 |

运行足够长的时间后，可得到极限概率：学生最终到家的概率约为 48%，到网吧的概率约为 52%。这可能与直觉不太相符，说明在处理统计问题时不能仅凭直觉。

需要注意的是，马尔可夫链方法得到的是理论概率，也可以使用随机数生成器模拟学生的游走过程来验证极限概率。

3. 线性方程组

3.1 线性方程组的矩阵表示

科学应用中常需要求解线性方程组，例如：
[
\begin{cases}
3x + 2y - z = 10 \
-x + 3y + 2z = 5 \
x - y - z = -1
\end{cases}
]
可将其表示为矩阵方程 (Ax = b)，其中
[
A =
\begin{bmatrix}
3 & 2 & -1 \
-1 & 3 & 2 \
1 & -1 & -1
\end{bmatrix},
x =
\begin{bmatrix}
x \
y \
z
\end{bmatrix},
b =
\begin{bmatrix}
10 \
5 \
-1
\end{bmatrix}
]
理论上，解为 (x = A^{-1}b)，其中 (A^{-1}) 是 (A) 的逆矩阵。

3.2 MATLAB 求解线性方程组

3.2.1 左除运算符求解

在 MATLAB 中，可使用左除运算符 \ 求解线性方程组。以下是具体代码：

A = [3 2 -1; -1 3 2; 1 -1 -1];
b = [10 5 -1]';
x = A \ b

运行结果为：

这意味着 (x = -2)，(y = 5)，(z = -6)。可以将矩阵运算 A \ b 理解为“(b) 除以 (A)” 或 “(A) 的逆矩阵乘以 (b)”。

虽然也可以使用 inv(A) * b 来求解，但 A \ b 更准确和高效，因为它使用高斯消元法，所需的浮点运算次数更少。例如，对于一个随机的 30×30 系统， a \ b 需要 24051 次浮点运算，而 inv(a) * b 需要 59850 次。

3.2.2 残差检查

求解线性方程组后，需要检查解的准确性，可通过计算残差 (r = A*x - b) 来进行。理论上，残差 (r) 应为零。在上述示例中，残差为：

r = A*x - b

结果为：

虽然残差的所有元素绝对值都小于 (10^{-15})，但仍可能存在潜在问题。

3.3 超定系统

当方程数量多于未知数数量时，系统称为超定系统，例如：
[
\begin{cases}
x - y = 0 \
y = 2 \
x = 1
\end{cases}
]
MATLAB 仍能给出此类系统的解。以下是求解代码：

A = [1 -1; 0 1; 1 0];
b = [0 2 1]';
x = A \ b

结果为：

x =
    1.3333
    1.6667

残差为：

r = A*x - b

结果为：

在这种情况下，MATLAB 给出的是最小二乘最佳拟合解，即使残差的模 (\sqrt{r(1)^2 + r(2)^2 + r(3)^2}) 最小的值。可以使用 sqrt(r' * r) 或 sqrt(sum(r .* r)) 计算该值（结果为 0.5774）。

3.4 欠定系统

当方程数量少于未知数数量时，系统称为欠定系统，此时有无限多个解，MATLAB 会找到一个使部分未知数为零的解。这类系统中的方程是线性规划问题的约束条件。

3.5 病态系统

有些线性方程组的系数可能是实验结果，存在误差。若系数的微小变化导致解的大幅变化，则称该系统为病态系统。例如：
[
\begin{cases}
10x + 7y + 8z + 7w = 32 \
7x + 5y + 6z + 5w = 23 \
8x + 6y + 10z + 9w = 33 \
7x + 5y + 9z + 10w = 31
\end{cases}
]
使用矩阵左除求解可得 (x = y = z = w = 1)，残差为零。但将方程组右边常数改为 32.1、22.9、32.9 和 31.1 后，解变为 (x = 6)，(y = -7.2)，(z = 2.9)，(w = -0.1)，残差仍然很小。

可使用 MATLAB 函数 rcond 测试系统的病态性，若系数矩阵 (A) 的 rcond(A) 接近零，则系统病态；接近 1 则系统良态。在本示例中，条件估计值约为 (2×10^{-4})，接近零，说明系统病态。也有作者建议，若矩阵的行列式相对于矩阵元素的值较小，则矩阵可能病态，本示例中矩阵 (A) 的行列式为 1，比大多数元素小一个数量级。

3.6 矩阵除法

3.6.1 左除

矩阵左除 A\B 定义为当 (B) 的行数与 (A) 相同时，形式上等同于 inv(A)*B ，但不直接计算逆矩阵。一般来说， x = A \ B 是方程组 (Ax = B) 的解。若 (A) 是方阵，使用高斯消元法进行左除；若 (A) 不是方阵，则以最小二乘法求解超定或欠定方程组。

3.6.2 右除

矩阵右除 B/A 定义为 (A'\B')' 。例如，对于上述线性方程组， x = (b'/a')' 会得到相同的解。

有时，使用 \ 或 / 计算超定或欠定系统的最小二乘解可能会有意外结果，例如：

a = [1 2];
b = [3 4];
a / b

结果为 ans = 0.4400 ，这是因为 a/b 等同于 (b'\a')' ，是 (b^T x^T = a^T) 的最小二乘解。而 a \ b 的结果为：

可通过写出完整的方程来理解。

4. 稀疏矩阵

4.1 稀疏矩阵概念

大型矩阵可能占用大量内存和计算时间，例如求解 (n) 个线性方程组需要 (n^2) 个矩阵元素，计算时间与 (n^3) 成正比。但有些矩阵非零元素相对较少，称为稀疏矩阵。MATLAB 提供了利用稀疏性的功能，可节省大量内存和处理时间。

4.2 创建稀疏矩阵

以第 2 节随机游走问题的转移概率矩阵 (P) 为例，它是稀疏矩阵的一个很好的候选，36 个元素中只有 10 个非零元素。若表示更多十字路口的矩阵会更稀疏，例如 100×100 的矩阵可能只有 1.98% 的非零元素。

在 MATLAB 中，可使用 sparse 函数创建稀疏矩阵。以下是将转移矩阵 (P) 表示为稀疏矩阵的代码：

n = 6;
P = sparse(1, 1, 1, n, n);
P = P + sparse(n, n, 1, n, n);
P = P + sparse(1:n-2, 2:n-1, 1/3, n, n);
P = P + sparse(3:n, 2:n-1, 2/3, n, n)

运行结果（使用 format rat ）为：

P =
   (1,1)       1
   (1,2)     1/3
   (3,2)     2/3
   (2,3)     1/3
   (4,3)     2/3
   (3,4)     1/3
   (5,4)     2/3
   (4,5)     1/3
   (6,5)     2/3
   (6,6)       1

若要显示稀疏矩阵的完整形式，可使用 full(P) ：

full(P)

结果（使用 format rat ）为：

ans =
       1     1/3         0         0         0         0
       0         0     1/3         0         0         0
       0     2/3         0     1/3         0         0
       0         0     2/3         0     1/3         0
       0         0         0     2/3         0         0
       0         0         0         0     2/3         1

sparse 函数的形式为 sparse(rows, cols, entries, m, n) ，用于生成一个 (m×n) 的稀疏矩阵，非零元素的下标为 (rows, cols) 。

4.3 稀疏矩阵的效率测试

创建一个 1000×1000 的单位矩阵 a = eye(1000) ，计算 a^2 大约需要 27.5 秒；而将其表示为稀疏矩阵 s = sparse(1:1000, 1:1000, 1, 1000, 1000) 后，计算 s^2 仅需约 0.05 秒，快了 500 多倍。

还可以使用 full(a) 将稀疏矩阵转换为完整形式，使用 sparse(a) 将完整矩阵转换为稀疏形式。函数 spy 可用于直观显示稀疏矩阵，可对矩阵 (P) 使用该函数，并将 (P) 扩大到约 50×50 后再次使用。

总结

矩阵左除运算符 \ 可直接求解线性方程组，但由于矩阵除法运算符 \ 和 / 在最小二乘解方法中可能产生意外结果，求解方程组时应始终计算残差。
处理大型且非零元素相对较少的矩阵时，应考虑使用 MATLAB 的稀疏矩阵功能。

练习

4.1 马尔可夫过程极限概率计算

计算第 2 节中学生从其余十字路口出发时的极限概率，验证他离网吧越近，最终到达网吧的可能性越大。直接计算 (P^{50})，观察第一行是否能看到极限概率。

4.2 线性方程组求解与检查

使用左除运算符求解方程组：
[
\begin{cases}
2x - y + z = 4 \
x + y + z = 3 \
3x - y - z = 1
\end{cases}
]
计算残差检查解的准确性，同时计算行列式 det 和条件估计值 rcond ，并得出结论。

4.3 病态系统验证

使用左除运算符求解方程组：
[
\begin{cases}
x + 5.000y = 17.0 \
1.5x + 7.501y = 25.503
\end{cases}
]
得到解 (x = 2)，(y = 3)，计算残差。将第二个方程右边的常数改为 25.501，再次求解并计算残差，观察解的变化。还可尝试将该常数改为 25.502、25.504 等，使用 rcond 计算条件估计值，使用 det 计算行列式，验证系统的病态性。也可对系数进行灵敏度分析，依次改变系数相同的小百分比，观察解的变化。

4.4 稀疏矩阵表示 Leslie 矩阵

使用 sparse 表示第 1 节中的 Leslie 矩阵，并通过预测 24 个月的兔子种群数量来测试该表示。

4.5 高斯消元法编程

若熟悉高斯消元法，可编写一个函数 x = mygauss(a, b) 直接对增广系数矩阵的行进行操作来求解一般系统 (Ax = b)。巧妙使用冒号运算符可使代码简化。对随机元素的 (A) 和 (b) 以及第 3 节和练习 4.4 中的系统进行测试，并使用左除检查解的正确性。

通过以上内容，我们了解了矩阵在种群动态、马尔可夫过程、线性方程组求解以及稀疏矩阵处理等方面的应用，掌握了相关的 MATLAB 实现方法和注意事项，同时也通过练习加深了对这些知识的理解和运用能力。

4.3 稀疏矩阵的效率测试

稀疏矩阵在处理大规模数据时，能显著提升计算效率和节省内存。以一个 1000×1000 的单位矩阵为例，我们可以清晰地看到稀疏矩阵和普通矩阵在计算效率上的巨大差异。

以下是具体的操作步骤：
1. 创建普通单位矩阵并计算平方

a = eye(1000);
tic; % 开始计时
a_squared = a^2;
toc; % 结束计时

在这个过程中，计算 a^2 大约需要 27.5 秒。这表明普通矩阵在处理大规模数据时，计算时间会显著增加。

将单位矩阵转换为稀疏矩阵并计算平方

s = sparse(1:1000, 1:1000, 1, 1000, 1000);
tic; % 开始计时
s_squared = s^2;
toc; % 结束计时

将单位矩阵转换为稀疏矩阵后，计算 s^2 仅需约 0.05 秒，比普通矩阵快了 500 多倍。这充分体现了稀疏矩阵在计算效率上的优势。

除了计算效率，稀疏矩阵和普通矩阵之间还可以相互转换：
- 使用 full(a) 可以将稀疏矩阵 a 转换为完整形式，且不会改变其稀疏表示。
- 使用 sparse(a) 可以将完整矩阵 a 转换为稀疏形式。

另外，函数 spy 可以直观地显示稀疏矩阵的结构。以下是使用示例：

spy(P); % 显示矩阵 P 的稀疏结构
P_large = zeros(50); % 创建一个 50×50 的矩阵
% 仿照之前的方式设置 P_large 的元素，使其成为类似的稀疏矩阵
for i = 3:50
    P_large(i,i-1) = 2/3;
    P_large(i-2,i-1) = 1/3;
end
P_large(1,1) = 1;
P_large(50,50) = 1;
spy(P_large); % 显示扩大后的矩阵的稀疏结构

总结

在实际应用中，矩阵的运用广泛且重要，以下是一些关键要点总结：
- 线性方程组求解 ：矩阵左除运算符 \ 是求解线性方程组的有效工具，但由于矩阵除法运算符 \ 和 / 在最小二乘解方法中可能产生意外结果，因此在求解方程组时，始终要计算残差来确保解的准确性。
- 稀疏矩阵应用 ：当处理大型矩阵且非零元素相对较少时，应优先考虑使用 MATLAB 的稀疏矩阵功能。稀疏矩阵能大幅节省内存和处理时间，提高计算效率。

练习

以下是一些相关练习，通过实践可以加深对矩阵知识的理解和运用能力：
1. 马尔可夫过程极限概率计算
- 计算第 2 节中学生从其余十字路口出发时的极限概率，验证他离网吧越近，最终到达网吧的可能性越大。
- 直接计算 (P^{50})，观察第一行是否能看到极限概率。
matlab n = 6; P = zeros(n); for i = 3:6 P(i,i-1) = 2/3; P(i-2,i-1) = 1/3; end P(1,1) = 1; P(6,6) = 1; P_50 = P^50; disp(P_50(1,:)); % 依次改变初始状态向量，计算不同起始位置的极限概率 for start = 2:5 x = zeros(6,1); x(start) = 1; for t = 1:100 % 迭代足够多次以接近极限 x = P * x; end disp(['Starting at intersection ', num2str(start), ':']); disp(x'); end
2. 线性方程组求解与检查
- 使用左除运算符求解方程组：
[
\begin{cases}
2x - y + z = 4 \
x + y + z = 3 \
3x - y - z = 1
\end{cases}
]
- 计算残差检查解的准确性，同时计算行列式 det 和条件估计值 rcond ，并得出结论。

A = [2 -1 1; 1 1 1; 3 -1 -1];
b = [4 3 1]';
x = A \ b;
r = A * x - b; % 计算残差
det_A = det(A); % 计算行列式
cond_est = rcond(A); % 计算条件估计值
disp('Solution:');
disp(x);
disp('Residual:');
disp(r);
disp(['Determinant: ', num2str(det_A)]);
disp(['Condition estimator: ', num2str(cond_est)]);

病态系统验证
- 使用左除运算符求解方程组：
  [
  \begin{cases}
  x + 5.000y = 17.0 \
  1.5x + 7.501y = 25.503
  \end{cases}
  ]
- 得到解 (x = 2)，(y = 3)，计算残差。
- 将第二个方程右边的常数改为 25.501，再次求解并计算残差，观察解的变化。
- 还可尝试将该常数改为 25.502、25.504 等，使用 rcond 计算条件估计值，使用 det 计算行列式，验证系统的病态性。
- 也可对系数进行灵敏度分析，依次改变系数相同的小百分比，观察解的变化。

A = [1 5; 1.5 7.501];
b_original = [17; 25.503];
x_original = A \ b_original;
r_original = A * x_original - b_original;
disp('Original solution:');
disp(x_original);
disp('Original residual:');
disp(r_original);

b_modified = [17; 25.501];
x_modified = A \ b_modified;
r_modified = A * x_modified - b_modified;
disp('Modified solution:');
disp(x_modified);
disp('Modified residual:');
disp(r_modified);

det_A = det(A);
cond_est = rcond(A);
disp(['Determinant: ', num2str(det_A)]);
disp(['Condition estimator: ', num2str(cond_est)]);

% 灵敏度分析
percentage_change = 0.01; % 1% 的变化
for i = 1:size(A, 1)
    for j = 1:size(A, 2)
        A_modified = A;
        A_modified(i,j) = A_modified(i,j) * (1 + percentage_change);
        x_sensitivity = A_modified \ b_original;
        disp(['Changed element (', num2str(i), ',', num2str(j), ') by ', num2str(percentage_change*100), '%:']);
        disp(x_sensitivity);
    end
end

稀疏矩阵表示 Leslie 矩阵
- 使用 sparse 表示第 1 节中的 Leslie 矩阵，并通过预测 24 个月的兔子种群数量来测试该表示。

n = 3;
L = sparse([2 1 3], [1 2 2], [1/3 9 1/2], n, n);
x = sparse([3], [1], [1], n, 1);
for t = 1:24
    x = L * x;
    disp([t full(x') sum(full(x))]);
end

高斯消元法编程
- 若熟悉高斯消元法，可编写一个函数 x = mygauss(a, b) 直接对增广系数矩阵的行进行操作来求解一般系统 (Ax = b)。巧妙使用冒号运算符可使代码简化。
- 对随机元素的 (A) 和 (b) 以及第 3 节和练习 4.4 中的系统进行测试，并使用左除检查解的正确性。

function x = mygauss(a, b)
    n = size(a, 1);
    Ab = [a b]; % 构造增广矩阵
    for i = 1:n-1
        % 选主元
        [~, pivot_row] = max(abs(Ab(i:n, i)));
        pivot_row = pivot_row + i - 1;
        if pivot_row ~= i
            temp = Ab(i, :);
            Ab(i, :) = Ab(pivot_row, :);
            Ab(pivot_row, :) = temp;
        end
        % 消元
        for j = i+1:n
            factor = Ab(j, i) / Ab(i, i);
            Ab(j, :) = Ab(j, :) - factor * Ab(i, :);
        end
    end
    % 回代求解
    x = zeros(n, 1);
    x(n) = Ab(n, n+1) / Ab(n, n);
    for i = n-1:-1:1
        x(i) = (Ab(i, n+1) - Ab(i, i+1:n) * x(i+1:n)) / Ab(i, i);
    end
end

% 测试
A = [3 2 -1; -1 3 2; 1 -1 -1];
b = [10 5 -1]';
x_gauss = mygauss(A, b);
x_left_division = A \ b;
disp('Gauss solution:');
disp(x_gauss);
disp('Left division solution:');
disp(x_left_division);

mermaid 流程图

graph LR
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A([开始]):::startend --> B{选择问题类型}:::decision
    B -->|种群动态| C(使用 Leslie 矩阵建模):::process
    B -->|马尔可夫过程| D(构建转移概率矩阵):::process
    B -->|线性方程组| E(使用矩阵左除求解):::process
    B -->|稀疏矩阵| F(使用 sparse 函数创建):::process
    C --> G(预测种群数量):::process
    D --> H(计算极限概率):::process
    E --> I(计算残差检查解):::process
    F --> J(测试计算效率):::process
    G --> K([结束]):::startend
    H --> K
    I --> K
    J --> K

通过以上的学习和练习，我们可以更深入地理解矩阵在不同领域的应用，掌握相关的 MATLAB 操作技巧，提高解决实际问题的能力。无论是处理种群动态、马尔可夫过程，还是求解线性方程组和处理稀疏矩阵，矩阵都发挥着重要的作用。希望大家通过不断实践，能够熟练运用这些知识，在实际应用中取得更好的效果。