MATLAB 中解线性方程组 \(Ax = b\) 的迭代法有多种实现方式。下面我将详细介绍常用的迭代方法、MATLAB 内置函数和自定义实现。
1. 迭代法简介
| 方法 | 特点 | 收敛条件 |
|---|---|---|
| 雅可比迭代 | 简单，并行性好 | 严格对角占优 |
| 高斯-赛德尔 | 收敛通常比雅可比快 | 严格对角占优 |
| 逐次超松弛 | 引入松弛因子加速 | 0 < ω < 2 |
2. MATLAB 内置函数
2.1 直接使用反斜杠（推荐优先尝试）
% 对于一般线性方程组
A = [4, -1, 0; -1, 4, -1; 0, -1, 4];
b = [1; 2; 3];
x = A \ b; % 直接法，通常最有效
2.2 迭代法求解器
% 使用 bicgstab（双共轭梯度稳定法）
x = bicgstab(A, b);% 使用 gmres（广义最小残差法）
x = gmres(A, b);% 使用 pcg（预处理共轭梯度法，对称正定矩阵）
x = pcg(A, b);
3. 自定义迭代法实现
3.1 雅可比迭代法 (Jacobi)
function [x, iter, rel_err] = jacobi(A, b, x0, max_iter, tol)% 雅可比迭代法求解 Ax = b% 输入:% A: 系数矩阵 (n×n)% b: 右端向量 (n×1)% x0: 初始猜测 (n×1)% max_iter: 最大迭代次数% tol: 容许误差% 输出:% x: 解向量% iter: 实际迭代次数% rel_err: 相对误差历史n = length(b);x = x0;x_new = zeros(n, 1);rel_err = zeros(max_iter, 1);D = diag(diag(A)); % 对角线部分R = A - D; % 非对角线部分for iter = 1:max_iter% 雅可比迭代公式: x_new = D^{-1}(b - R*x)x_new = (b - R * x) ./ diag(D);% 计算相对误差rel_err(iter) = norm(x_new - x, 2) / norm(x_new, 2);% 检查收敛if rel_err(iter) < tolx = x_new;rel_err = rel_err(1:iter);fprintf('Jacobi 迭代在 %d 次后收敛\n', iter);return;endx = x_new;endwarning('Jacobi: 达到最大迭代次数 %d 但未收敛\n', max_iter);rel_err = rel_err(1:iter);
end
3.2 高斯-赛德尔迭代法 (Gauss-Seidel)
function [x, iter, rel_err] = gauss_seidel(A, b, x0, max_iter, tol)% 高斯-赛德尔迭代法求解 Ax = b% 输入输出同雅可比法n = length(b);x = x0;rel_err = zeros(max_iter, 1);L = tril(A, -1); % 严格下三角D = diag(diag(A)); % 对角线U = triu(A, 1); % 严格上三角for iter = 1:max_iterx_old = x;% 高斯-赛德尔迭代: (L+D)x_new = b - U*xfor i = 1:nsum1 = 0;for j = 1:i-1sum1 = sum1 + A(i, j) * x(j);endsum2 = 0;for j = i+1:nsum2 = sum2 + A(i, j) * x_old(j);endx(i) = (b(i) - sum1 - sum2) / A(i, i);end% 计算相对误差rel_err(iter) = norm(x - x_old, 2) / norm(x, 2);% 检查收敛if rel_err(iter) < tolrel_err = rel_err(1:iter);fprintf('Gauss-Seidel 迭代在 %d 次后收敛\n', iter);return;endendwarning('Gauss-Seidel: 达到最大迭代次数 %d 但未收敛\n', max_iter);rel_err = rel_err(1:iter);
end
3.3 逐次超松弛迭代法 (SOR)
function [x, iter, rel_err] = sor(A, b, x0, omega, max_iter, tol)% 逐次超松弛迭代法 (SOR) 求解 Ax = b% 输入:% omega: 松弛因子 (0 < omega < 2)% 其他参数同前n = length(b);x = x0;rel_err = zeros(max_iter, 1);if omega <= 0 || omega >= 2error('松弛因子必须在 (0, 2) 范围内');endfor iter = 1:max_iterx_old = x;for i = 1:nsigma = 0;for j = 1:nif j ~= isigma = sigma + A(i, j) * x(j);endend% SOR 迭代公式x(i) = (1 - omega) * x(i) + omega * (b(i) - sigma) / A(i, i);end% 计算相对误差rel_err(iter) = norm(x - x_old, 2) / norm(x, 2);% 检查收敛if rel_err(iter) < tolrel_err = rel_err(1:iter);fprintf('SOR (ω=%.2f) 在 %d 次后收敛\n', omega, iter);return;endendwarning('SOR: 达到最大迭代次数 %d 但未收敛\n', max_iter);rel_err = rel_err(1:iter);
end
4. 测试示例
4.1 定义测试问题
% 创建严格对角占优矩阵（保证迭代法收敛）
n = 100;
A = diag(4 * ones(n, 1)) + diag(-1 * ones(n-1, 1), 1) + diag(-1 * ones(n-1, 1), -1);% 生成精确解
x_exact = ones(n, 1);% 计算右端项
b = A * x_exact;% 设置初始猜测和参数
x0 = zeros(n, 1);
max_iter = 1000;
tol = 1e-6;
omega = 1.2; % SOR 松弛因子
4.2 运行不同迭代法
% 雅可比迭代
tic;
[x_jacobi, iter_j, err_j] = jacobi(A, b, x0, max_iter, tol);
time_j = toc;% 高斯-赛德尔迭代
tic;
[x_gs, iter_gs, err_gs] = gauss_seidel(A, b, x0, max_iter, tol);
time_gs = toc;% SOR 迭代
tic;
[x_sor, iter_sor, err_sor] = sor(A, b, x0, omega, max_iter, tol);
time_sor = toc;% 直接法对比
tic;
x_direct = A \ b;
time_direct = toc;
4.3 比较结果
% 计算误差
err_jacobi = norm(x_jacobi - x_exact, 2);
err_gs = norm(x_gs - x_exact, 2);
err_sor = norm(x_sor - x_exact, 2);
err_direct = norm(x_direct - x_exact, 2);% 显示结果
fprintf('\n=========== 迭代法性能比较 ===========\n');
fprintf('方法\t\t迭代次数\t时间(s)\t\t相对误差\n');
fprintf('--------------------------------------------\n');
fprintf('Jacobi\t\t%d\t\t%.4f\t\t%.2e\n', iter_j, time_j, err_jacobi);
fprintf('Gauss-Seidel\t%d\t\t%.4f\t\t%.2e\n', iter_gs, time_gs, err_gs);
fprintf('SOR (ω=%.1f)\t%d\t\t%.4f\t\t%.2e\n', omega, iter_sor, time_sor, err_sor);
fprintf('直接法\\\t\t-\t\t%.4f\t\t%.2e\n', time_direct, err_direct);% 绘制收敛历史
figure;
semilogy(1:length(err_j), err_j, 'b-', 'LineWidth', 2, 'DisplayName', 'Jacobi');
hold on;
semilogy(1:length(err_gs), err_gs, 'r-', 'LineWidth', 2, 'DisplayName', 'Gauss-Seidel');
semilogy(1:length(err_sor), err_sor, 'g-', 'LineWidth', 2, 'DisplayName', sprintf('SOR (ω=%.1f)', omega));
grid on;
xlabel('迭代次数');
ylabel('相对误差 (log scale)');
title('迭代法收敛历史比较');
legend('Location', 'best');
5. 实用技巧
5.1 预处理技术
% 对角预处理（简单但有效）
D = diag(diag(A));
A_precond = D^(-1/2) * A * D^(-1/2);
b_precond = D^(-1/2) * b;% 使用预处理
x = jacobi(A_precond, b_precond, x0, max_iter, tol);
5.2 收敛性检查函数
function is_convergent = check_convergence(A, method)% 检查迭代法收敛性% method: 'jacobi' 或 'gauss_seidel'n = size(A, 1);if strcmp(method, 'jacobi')D = diag(diag(A));R = A - D;T = -D \ R; % 迭代矩阵elseif strcmp(method, 'gauss_seidel')L = tril(A, -1);D = diag(diag(A));U = triu(A, 1);T = -(D + L) \ U; % 迭代矩阵elseerror('未知方法');endrho = max(abs(eig(T))); % 谱半径is_convergent = (rho < 1);fprintf('谱半径 ρ = %.4f, 收敛: %s\n', rho, string(is_convergent));
end
5.3 自适应 SOR
function [x, omega_opt] = adaptive_sor(A, b, x0, max_iter, tol)% 自适应选择最优松弛因子n = length(b);omegas = 1.0:0.1:1.9; % 测试不同的 ωbest_iter = inf;omega_opt = 1.0;% 测试不同 ωfor omega = omegas[~, iter, ~] = sor(A, b, x0, omega, max_iter, tol);if iter < best_iterbest_iter = iter;omega_opt = omega;endendfprintf('最优松弛因子: ω = %.2f\n', omega_opt);[x, ~, ~] = sor(A, b, x0, omega_opt, max_iter, tol);
end
参考代码 matlab解线性方程组的迭代法 www.youwenfan.com/contentcnt/122244.html
6. 性能优化建议
- 向量化计算：避免内层循环
- 稀疏矩阵：使用
sparse存储 - 预处理：使用
ichol或ilu预处理 - 并行计算：使用
parfor加速
% 向量化的雅可比迭代
function x = jacobi_vectorized(A, b, x0, max_iter, tol)n = length(b);x = x0;d = diag(A);for iter = 1:max_iterx_new = x + (b - A * x) ./ d; % 向量化计算if norm(x_new - x) / norm(x_new) < tolx = x_new;break;endx = x_new;end
end
7. 完整使用示例
% 1. 准备问题
n = 500;
A = gallery('poisson', sqrt(n)); % 生成泊松问题矩阵
b = rand(n, 1);
x0 = zeros(n, 1);
max_iter = 1000;
tol = 1e-8;% 2. 检查收敛性
fprintf('检查收敛性:\n');
check_convergence(A, 'jacobi');
check_convergence(A, 'gauss_seidel');% 3. 求解
[x_j, iter_j] = jacobi(A, b, x0, max_iter, tol);
[x_gs, iter_gs] = gauss_seidel(A, b, x0, max_iter, tol);
[x_sor, omega] = adaptive_sor(A, b, x0, max_iter, tol);% 4. 验证结果
residual_j = norm(b - A * x_j) / norm(b);
residual_gs = norm(b - A * x_gs) / norm(b);
residual_sor = norm(b - A * x_sor) / norm(b);fprintf('\n残差:\n');
fprintf('Jacobi: %.2e\n', residual_j);
fprintf('Gauss-Seidel: %.2e\n', residual_gs);
fprintf('SOR: %.2e\n', residual_sor);
总结
| 情况 | 推荐方法 |
|---|---|
| 小规模稠密矩阵 | 直接法 A\b |
| 大规模稀疏矩阵 | 预处理共轭梯度法 pcg |
| 非对称矩阵 | GMRES 或 BiCGSTAB |
| 教学/简单问题 | 雅可比或高斯-赛德尔 |
| 需要最佳收敛速度 | SOR 或 自适应 SOR |
这些 MATLAB 实现可以直接使用，也可以根据具体问题调整参数。