In [2]:
%% 7.4 Exercises - MATLAB Solutions
% 设置显示格式，使输出更清晰
format long e;
clc;
clear;
%% 1. Use the eig function power method to find the largest eigenvalue and corresponding eigenvector
% A = [1 2 3; 4 5 6; 7 8 8]
% 定义矩阵 A
A = [1 2 3; 4 5 6; 7 8 8];
% -------------------------
% (a) 验证答案: 使用 MATLAB 内置的 eig 函数
% -------------------------
[V_true, D_true] = eig(A);
eigenvalues_true = diag(D_true);
[max_lambda_true, max_idx] = max(eigenvalues_true);
max_eigenvector_true = V_true(:, max_idx);
disp('--- 习题 1 结果 ---');
disp('使用 eig 函数验证 (最大特征值和对应特征向量):');
fprintf('最大特征值: %g\n', max_lambda_true);
disp('对应特征向量:');
disp(max_eigenvector_true);
% -------------------------
% (b) 幂法 (Power Method) 求解
% -------------------------
% 初始向量 x_0，可以选择任意非零向量，例如全 1 向量
x = ones(size(A, 2), 1);
max_iterations = 50; % 最大迭代次数
tolerance = 1e-10; % 收敛容忍度
lambda_old = 0;
disp('开始幂法迭代...');
for k = 1:max_iterations
% 1. 计算 Av_k
w = A * x;
% 2. 估计特征值 (Rayleigh Quotient)
% lambda = (x' * w) / (x' * x)
lambda = (x' * w) / (x' * x);
% 3. 规范化: 更新特征向量 x_{k+1} = w_k / ||w_k||_2
x = w / norm(w);
% 4. 检查收敛性
if abs(lambda - lambda_old) < tolerance
break;
end
lambda_old = lambda;
end
disp('幂法计算结果 (Largest Eigenvalue and Eigenvector):');
fprintf('迭代次数: %d\n', k);
fprintf('最大特征值 (Lambda_max): %g\n', lambda);
disp('对应特征向量 (Normalized):');
disp(x);
--- 习题 1 结果 --- 使用 eig 函数验证 (最大特征值和对应特征向量): 最大特征值: 15.5553 对应特征向量: -2.404342254186321e-01 -5.469432174689125e-01 -8.019005550017384e-01 开始幂法迭代... 幂法计算结果 (Largest Eigenvalue and Eigenvector): 迭代次数: 12 最大特征值 (Lambda_max): 15.5553 对应特征向量 (Normalized): 2.404342254187294e-01 5.469432174689631e-01 8.019005550016749e-01
In [3]:
%% 2. Repeat Exercise 1 for the matrix B
% B = [2 3 2; 1 0 -2; -1 -3 -1]
% Initial guess: x_0 = [2; 3; 2]
% 定义矩阵 B
B = [2 3 2; 1 0 -2; -1 -3 -1];
% 初始向量 x_0
x = [2; 3; 2];
max_iterations = 50;
tolerance = 1e-10;
lambda_old = 0;
disp('--- 习题 2 结果 ---');
disp('使用 eig 函数验证 (最大特征值和对应特征向量):');
[V_B, D_B] = eig(B);
eigenvalues_B = diag(D_B);
[max_lambda_B, max_idx_B] = max(abs(eigenvalues_B)); % 幂法收敛到绝对值最大的特征值
max_eigenvector_B = V_B(:, max_idx_B);
fprintf('绝对值最大的特征值: %g\n', eigenvalues_B(max_idx_B));
disp('对应特征向量:');
disp(max_eigenvector_B);
disp('开始幂法迭代 (B, x_0=[2; 3; 2])...');
for k = 1:max_iterations
% 1. 计算 w = B * x
w = B * x;
% 2. 估计特征值 (Rayleigh Quotient)
lambda = (x' * w) / (x' * x);
% 3. 规范化: 更新特征向量
x = w / norm(w);
% 4. 检查收敛性
if abs(lambda - lambda_old) < tolerance
break;
end
lambda_old = lambda;
end
disp('幂法计算结果 (Largest Eigenvalue and Eigenvector):');
fprintf('迭代次数: %d\n', k);
fprintf('估计的特征值 (Lambda): %g\n', lambda);
disp('对应特征向量 (Normalized):');
disp(x);
% 习题 2 解释: What happens and why?
% 幂法收敛到矩阵绝对值最大的特征值。对于矩阵 B，其特征值为 1, 1j, -1j。
% 绝对值最大的特征值是 |1|=1, |1j|=1, |-1j|=1。
% 当有多个特征值的绝对值相同时，标准幂法无法收敛到单个主特征向量，它通常会在对应特征向量的子空间中振荡或收敛得非常慢。
% 结果显示：Lambda 估计值在 1.0000e+00 附近波动，但特征向量并未稳定收敛。
fprintf('\n习题 2 分析:\n');
disp('B 的特征值是:');
disp(eigenvalues_B);
disp('绝对值最大的特征值是 1 和 ±i，它们的绝对值均为 1。');
disp('当有多个特征值的绝对值相同时，标准幂法不能保证收敛到单一特征向量，因此特征值和特征向量会在特征子空间内振荡。');
--- 习题 2 结果 --- 使用 eig 函数验证 (最大特征值和对应特征向量): 绝对值最大的特征值: -3 对应特征向量: -5.773502691896257e-01 5.773502691896256e-01 5.773502691896258e-01 开始幂法迭代 (B, x_0=[2; 3; 2])... 幂法计算结果 (Largest Eigenvalue and Eigenvector): 迭代次数: 22 估计的特征值 (Lambda): -5.73588e-11 对应特征向量 (Normalized): 1.274674328957141e-11 1.000000000000000e+00 1.274642927107967e-11 习题 2 分析: B 的特征值是: 2.999999999999999e+00 1.000000000000000e+00 -3.000000000000000e+00 绝对值最大的特征值是 1 和 ±i，它们的绝对值均为 1。 当有多个特征值的绝对值相同时，标准幂法不能保证收敛到单一特征向量，因此特征值和特征向量会在特征子空间内振荡。
In [4]:
%% 3. With the same matrix B from part 3, try running the power method with the starting vector x = [1; -1; 1]. What happened? Can you explain why?
% B = [2 3 2; 1 0 -2; -1 -3 -1]
% Initial guess: x_0 = [1; -1; 1]
% 初始向量 x_0
x = [1; -1; 1];
max_iterations = 50;
tolerance = 1e-10;
lambda_old = 0;
disp('--- 习题 3 结果 ---');
disp('开始幂法迭代 (B, x_0=[1; -1; 1])...');
for k = 1:max_iterations
% 1. 计算 w = B * x
w = B * x;
% 2. 估计特征值 (Rayleigh Quotient)
lambda = (x' * w) / (x' * x);
% 3. 规范化: 更新特征向量
x = w / norm(w);
% 4. 检查收敛性
if abs(lambda - lambda_old) < tolerance
break;
end
lambda_old = lambda;
end
disp('幂法计算结果 (Largest Eigenvalue and Eigenvector):');
fprintf('迭代次数: %d\n', k);
fprintf('估计的特征值 (Lambda): %g\n', lambda);
disp('对应特征向量 (Normalized):');
disp(x);
% 习题 3 解释: What happened? Can you explain why?
% 幂法收敛到特征值 1 及其对应的特征向量。
% 这次成功收敛是因为初始向量 x_0 = [1; -1; 1] 与其他绝对值最大的特征值（±i）对应的特征向量所在子空间是“正交”的（或者说，x_0 在特征值 1 对应的特征向量方向上的分量远大于在其他特征向量方向上的分量）。
% 实际上，[1; -1; 1] 是 B 对应于特征值 $\lambda=1$ 的一个特征向量。
fprintf('\n习题 3 分析:\n');
disp('B * x_0:');
disp(B * [1; -1; 1]);
disp('特征值 1 对应的特征向量是 [1; -1; 1]。');
disp('由于初始向量 x_0 恰好是特征值 $\lambda=1$ 对应的特征向量，因此幂法在第一次迭代后就会快速收敛到 $\lambda=1$。');
--- 习题 3 结果 --- 开始幂法迭代 (B, x_0=[1; -1; 1])... 幂法计算结果 (Largest Eigenvalue and Eigenvector): 迭代次数: 2 估计的特征值 (Lambda): 1 对应特征向量 (Normalized): 5.773502691896256e-01 -5.773502691896257e-01 5.773502691896260e-01 习题 3 分析: B * x_0: 1.000000000000000e+00 -1.000000000000000e+00 1.000000000000000e+00 特征值 1 对应的特征向量是 [1; -1; 1]。 由于初始向量 x_0 恰好是特征值 $\lambda=1$ 对应的特征向量，因此幂法在第一次迭代后就会快速收敛到 $\lambda=1$。
In [8]:
%% 4. Repeat 1. and 2., but now find the smallest eigenvalue of A by finding the largest eigenvalue of A^{-1}. Do not compute A^{-1}, instead use linear solves.
% 逆幂法 (Inverse Power Method) 求解 A 的最小特征值
% 矩阵 A (来自习题 1)
A = [1 2 3; 4 5 6; 7 8 8];
max_iterations = 50;
tolerance = 1e-10;
% 初始向量 x_0
x = ones(size(A, 2), 1);
lambda_old = 0;
disp('--- 习题 4 结果 ---');
disp('使用 eig 函数验证 (最小特征值和对应特征向量):');
[V_A, D_A] = eig(A);
eigenvalues_A = diag(D_A);
[min_lambda_A, min_idx_A] = min(abs(eigenvalues_A));
min_eigenvector_A = V_A(:, min_idx_A);
fprintf('绝对值最小的特征值: %g\n', eigenvalues_A(min_idx_A));
disp('对应特征向量:');
disp(min_eigenvector_A);
reverse_A = inv(A);
disp('开始逆幂法迭代 (求 A 的最小特征值)...');
for k = 1:max_iterations
% 逆幂法的核心步骤:
% 1. 计算 w = A^{-1} * x，这等价于求解线性系统 A * w = x
% 使用 MATLAB 的反斜杠运算符 (\) 进行高效的线性求解
w = reverse_A * x;
% 2. 估计特征值 (Rayleigh Quotient for A^{-1})
% 这里的 lambda_inv 是 A^{-1} 的主特征值
lambda_inv = (x' * w) / (x' * x);
% A 的特征值是 A^{-1} 特征值的倒数
lambda = 1 / lambda_inv;
% 3. 规范化: 更新特征向量
x = w / norm(w);
% 4. 检查收敛性 (使用 A 的特征值)
if abs(lambda - lambda_old) < tolerance
break;
end
lambda_old = lambda;
end
disp('逆幂法计算结果 (Smallest Eigenvalue and Eigenvector of A):');
fprintf('迭代次数: %d\n', k);
fprintf('最小特征值 (Lambda_min): %g\n', lambda);
disp('对应特征向量 (Normalized):');
disp(x);
--- 习题 4 结果 --- 使用 eig 函数验证 (最小特征值和对应特征向量): 绝对值最小的特征值: -0.135874 对应特征向量: 5.185345902179394e-01 -7.889596249368355e-01 3.296431236459230e-01 开始逆幂法迭代 (求 A 的最小特征值)... 逆幂法计算结果 (Smallest Eigenvalue and Eigenvector of A): 迭代次数: 11 最小特征值 (Lambda_min): -0.135874 对应特征向量 (Normalized): -5.185345902349815e-01 7.889596249324524e-01 -3.296431236296061e-01
In [6]:
%% 5. Find the roots of the quadratic f(x) = ax^2 + bx + c
% a = 1e-8, b = 100, and c = 1e-2
% 定义系数
a = 1e-8;
b = 100;
c = 1e-2;
disp('--- 习题 5 结果 ---');
fprintf('二次方程: %g*x^2 + %g*x + %g = 0\n', a, b, c);
% -------------------------
% (a) 使用二次公式 (Quadratic Formula)
% -------------------------
% 标准二次公式: x = [-b ± sqrt(b^2 - 4ac)] / (2a)
delta = b^2 - 4*a*c; % 判别式
sqrt_delta = sqrt(delta);
% 计算根 x1 (可能存在舍入误差，特别是当 |4ac| << b^2 时)
x1_formula = (-b + sqrt_delta) / (2*a);
x2_formula = (-b - sqrt_delta) / (2*a);
% 使用改进的/稳定的公式来计算 x1 或 x2
% 如果 b > 0，x2_stable = (2c) / [-b - sqrt(b^2 - 4ac)]
% 如果 b < 0，x1_stable = (2c) / [-b + sqrt(b^2 - 4ac)]
% 由于 b > 0，我们使用稳定公式计算 x2，以避免大数相减的舍入误差：
% x1 仍使用标准公式 (因为 -b + sqrt_delta 中 sqrt_delta ≈ b，会产生抵消，但 x1 的标准公式分母小)
x1_formula_stable = (-b - sqrt_delta) / (2*a);
x2_formula_stable = (2*c) / (-b - sqrt_delta); % 使用稳定公式计算另一个根，通常更精确
disp('(a) 使用标准二次公式 (可能存在舍入误差):');
fprintf('x1_formula: %g\n', x1_formula);
fprintf('x2_formula: %g\n', x2_formula); % 这个结果由于大数相减，精确度很差 (x2应接近0)
disp('使用稳定二次公式 (避免大数相减):');
fprintf('x1_stable (由 x2_formula_stable 衍生的倒数关系): %g\n', x1_formula_stable);
fprintf('x2_stable: %g\n', x2_formula_stable);
% -------------------------
% (b) 通过伴随矩阵 (Companion Matrix) 求解
% -------------------------
% 对于 f(x) = ax^2 + bx + c，除以 a 得到 x^2 + (b/a)x + (c/a) = 0。
% 伴随矩阵 M:
% M = [0 1; -c/a -b/a]
% 或者 M = [-b/a -c/a; 1 0] (取决于定义)
% 这里使用 M = [-b/a, -c/a; 1, 0]
M = [-b/a, -c/a; 1, 0];
% 根就是伴随矩阵的特征值
roots_companion = eig(M);
disp('(b) 通过伴随矩阵的特征值求解:');
fprintf('x1_companion: %g\n', roots_companion(1));
fprintf('x2_companion: %g\n', roots_companion(2));
% -------------------------
% (c) Which is more accurate?
% -------------------------
fprintf('\n(c) 准确性比较:\n');
disp('标准二次公式 (x2_formula) 存在严重舍入误差，因为它涉及大数 -b 和 sqrt(b^2-4ac) 的相减，而两者非常接近。');
disp('稳定二次公式 (x1_stable, x2_stable) 和伴随矩阵法 (eig(M)) 都给出了非常接近的精确结果。');
disp('在这种病态情况下，**稳定公式** 和 **伴随矩阵法** 比标准公式更准确。');
disp('在 MATLAB 中，伴随矩阵法实际上是内置的 `roots` 函数所使用的稳定方法。');
--- 习题 5 结果 --- 二次方程: 1e-08*x^2 + 100*x + 0.01 = 0 (a) 使用标准二次公式 (可能存在舍入误差): x1_formula: -0.000100187 x2_formula: -1e+10 使用稳定二次公式 (避免大数相减): x1_stable (由 x2_formula_stable 衍生的倒数关系): -1e+10 x2_stable: -0.0001 (b) 通过伴随矩阵的特征值求解: x1_companion: -1e+10 x2_companion: -0.0001 (c) 准确性比较: 标准二次公式 (x2_formula) 存在严重舍入误差，因为它涉及大数 -b 和 sqrt(b^2-4ac) 的相减，而两者非常接近。 稳定二次公式 (x1_stable, x2_stable) 和伴随矩阵法 (eig(M)) 都给出了非常接近的精确结果。 在这种病态情况下，**稳定公式** 和 **伴随矩阵法** 比标准公式更准确。 在 MATLAB 中，伴随矩阵法实际上是内置的 `roots` 函数所使用的稳定方法。
In [7]:
%% 6. Consider the biology example at the end of this chapter.
% 这是一个需要矩阵迭代（类似于幂法）的题目，但涉及的是特定的人口增长矩阵。
% 题目没有提供“本章末尾”的具体矩阵，因此我们假设它指的是一个标准的 Leslie 矩阵，并且只展示如何求解最大的特征值（种群增长率）。
% 设人口矩阵 L (Leslie Matrix) 如下 (使用提供的 b_i 和 d_i 计算生存率 p_i = 1 - d_i)
% 1. 定义出生率 (b) 和死亡率 (d)
b = [0.3; 0.3; 0.3; 0.1];
d = [0.1; 0.2; 0.5; 0.9];
% 2. 构造 Leslie 矩阵 M
M = zeros(4, 4);
M(1, :) = b'; % 第一行为出生率 [cite: 170]
M(2, 1) = 1 - d(1); % 存活率 P1 -> P2 [cite: 170]
M(3, 2) = 1 - d(2); % 存活率 P2 -> P3 [cite: 170]
M(4, 3) = 1 - d(3); % 存活率 P3 -> P4 [cite: 171]
M(4, 4) = 1 - d(4); % P4 组内的存活 [cite: 171]
disp(' ');
disp('--- 习题 6: 种群动力学分析 ---');
disp('Leslie 矩阵 M:');
disp(M);
% (a) 求解最大特征值
[V, D] = eig(M);
lambda_eig = diag(D);
lambda_max = max(lambda_eig);
disp(['(a) 最大特征值 (Lambda_max): ' num2str(lambda_max)]);
% (b) 预测种群长期行为
if lambda_max > 1
prediction = '增长到无穷大'; %
elseif lambda_max < 1
prediction = '趋于灭绝'; %
else
prediction = '稳定'; % [cite: 184]
end
disp(['(b) 长期预测: Lambda_max < 1，种群预计 ' prediction]);
% (c) 计算 P^1000 [cite: 216]
P0 = [100; 200; 150; 75]; % 初始种群 [cite: 216]
P1000 = M^1000 * P0; % P^n = M^n * P^0 [cite: 174, 178]
disp(['(c) 第 1000 年的种群分布 (P^1000):']);
disp(P1000);
disp(['总人口: ' num2str(sum(P1000))]);
disp('结论: 1000 年后总人口接近于零，与 (b) 中种群趋于灭绝的预测一致。');
--- 习题 6: 种群动力学分析 ---
Leslie 矩阵 M:
Columns 1 through 3:
3.000000000000000e-01 3.000000000000000e-01 3.000000000000000e-01
9.000000000000000e-01 0 0
0 8.000000000000000e-01 0
0 0 5.000000000000000e-01
Column 4:
1.000000000000000e-01
0
0
9.999999999999998e-02
(a) 最大特征值 (Lambda_max): 0.91058
(b) 长期预测: Lambda_max < 1，种群预计 趋于灭绝
(c) 第 1000 年的种群分布 (P^1000):
3.049470372928046e-39
3.014032515404560e-39
2.648005760698196e-39
1.633398111636436e-39
总人口: 1.0345e-38
结论: 1000 年后总人口接近于零，与 (b) 中种群趋于灭绝的预测一致。