3. MATLAB矩阵运算求解线性方程组

问题设定与三种求解路径

本节求解一个来自结构静力分析的 4 节点线性方程组 A·x = b，其中 A 为对称正定刚度矩阵，b 为节点载荷向量。三种方法分别对应 MATLAB 中不同的数值路径：矩阵左除（\）、显式 LU 分解（lu）和 Jacobi 迭代。目标是得到同一解向量 x，并比较残差范数、条件数与耗时。

构造测试矩阵与基准解

为便于验证，先构造一个已知精确解 x_true = [1; 2; 3; 4] 的对称正定矩阵。矩阵 A 由对角占优的刚度形式生成，b = A * x_true，这样三种方法求出的 x 都应与 x_true 对比。

% 参数区
n = 4;
x_true = [1; 2; 3; 4];

% 构造对称正定矩阵（对角占优）
A = [ 10  -1   2   0;
      -1  11  -1   3;
       2  -1  10  -2;
       0   3  -2   9 ];

b = A * x_true;   % 精确右端项

% 检查对称性与正定性
fprintf('对称性误差: %.2e\n', norm(A - A', 'fro'));
fprintf('最小特征值: %.4f\n', min(eig(A)));

运行后应看到对称性误差为 0，最小特征值大于 0，说明矩阵对称正定，三种方法理论上均可收敛。

方法一：矩阵左除

MATLAB 的 \ 运算符对一般方阵自动执行带部分选主元的 LU 分解，是工程上最常用的直接解法。

tic;
x_backslash = A \ b;
t_backslash = toc;

res_backslash = norm(A * x_backslash - b);
err_backslash = norm(x_backslash - x_true);

fprintf('左除解: '); disp(x_backslash');
fprintf('残差: %.2e  误差: %.2e  耗时: %.6f s\n', ...
        res_backslash, err_backslash, t_backslash);

方法二：显式 LU 分解

当需要对同一系数矩阵 A 反复求解不同右端项时，显式 LU 分解可以复用分解结果，避免重复分解开销。

tic;
[L, U, P] = lu(A);          % P*A = L*U
y = L \ (P * b);            % 前代
x_lu = U \ y;               % 回代
t_lu = toc;

res_lu = norm(A * x_lu - b);
err_lu = norm(x_lu - x_true);

fprintf('LU解: '); disp(x_lu');
fprintf('残差: %.2e  误差: %.2e  耗时: %.6f s\n', ...
        res_lu, err_lu, t_lu);

方法三：Jacobi 迭代

Jacobi 迭代将 A 拆分为对角部分 D 和非对角部分 R，迭代格式为 x(k+1) = D⁻¹(b − R·x(k))。收敛条件为迭代矩阵谱半径小于 1，对本例的对角占优矩阵成立。

D = diag(diag(A));
R = A - D;
D_inv = inv(D);

x_jac = zeros(n, 1);
max_iter = 500;
tol = 1e-12;
res_hist = zeros(max_iter, 1);

tic;
for k = 1:max_iter
    x_jac = D_inv * (b - R * x_jac);
    res_hist(k) = norm(A * x_jac - b);
    if res_hist(k) < tol
        break;
    end
end
t_jac = toc;

res_jac = norm(A * x_jac - b);
err_jac = norm(x_jac - x_true);

fprintf('Jacobi解: '); disp(x_jac');
fprintf('迭代次数: %d\n', k);
fprintf('残差: %.2e  误差: %.2e  耗时: %.6f s\n', ...
        res_jac, err_jac, t_jac);

结果对比与选型判断