如何高效求解大规模Hilbert矩阵？
摘要：博客地址：https:www.cnblogs.comzylyehuo 运行效果 运行代码 clear; clc; close all; % 参数定义 n = 2; % Hilbert矩阵的维数 maxIter = 40; % 最大迭
博客地址：https://www.cnblogs.com/zylyehuo/
运行效果
运行代码
clear; clc; close all;
% ============================================ 参数定义 ============================================
n = 2; % Hilbert矩阵的维数
maxIter = 40; % 最大迭代次数
omega = 0.75; % SOR松弛因子
epsilon = 1e-2; % 迭代绝对误差限
% ============================================ 主函数调用 ============================================
main(n, maxIter, omega, epsilon);
% ============================================ SOR松弛因子优化模块 ============================================
fprintf('\n\n=========================================== SOR松弛因子优化分析 ===========================================\n');
findOptimalOmega(n, maxIter, epsilon, omega);
% ============================================ 主函数模块 ============================================
function main(n, maxIter, omega, epsilon)
fprintf('=========================================== 主程序开始 ===========================================\n');
% 生成 Hilbert 矩阵
H = myHilbert(n);
fprintf('生成的 Hilbert(%d) 矩阵为：\n', n);
disp(H);
fprintf('------------------------------------------------------------------------------------------------------------------------\n');
x_exact = ones(n, 1);
fprintf('精确解 x* = [' );
fprintf('%.16f, ', x_exact(1:end-1));
fprintf('%.16f]^T\n', x_exact(end));
% 计算右端向量 b = H * x_exact
b = H * x_exact;
fprintf('右端向量 b = [' );
fprintf('%.16f, ', b(1:end-1));
fprintf('%.16f]^T\n', b(end));
fprintf('------------------------------------------------------------------------------------------------------------------------\n');
% 高斯消去法
x_gauss = gaussianElimination(H, b);
err_gauss = norm(x_gauss - x_exact);
fprintf('高斯消去法求得解 x = [' );
fprintf('%.16f, ', x_gauss(1:end-1));
fprintf('%.16f]^T\n', x_gauss(end));
fprintf('高斯消去法绝对误差 = %.6e\n', err_gauss);
fprintf('------------------------------------------------------------------------------------------------------------------------\n');
% Jacobi 迭代
[x_jacobi, err_history_jacobi, iter_jacobi, spectral_radius_jacobi] = jacobiIteration(H, b, maxIter, epsilon, x_exact);
% fprintf('Jacobi 迭代法最终解 x = [' );
% fprintf('%.16f, ', x_jacobi(1:end-1));
% fprintf('%.16f]^T\n', x_jacobi(end));
% Gauss-Seidel 迭代
[x_gs, err_history_gs, iter_gs, spectral_radius_gs] = gaussSeidelIteration(H, b, maxIter, epsilon, x_exact);
% fprintf('Gauss-Seidel 迭代法最终解 x = [' );
% fprintf('%.16f, ', x_gs(1:end-1));
% fprintf('%.16f]^T\n', x_gs(end));
% SOR 迭代
[x_sor, err_history_sor, iter_sor, spectral_radius_sor] = SORIteration(H, b, omega, maxIter, epsilon, x_exact);
% fprintf('SOR 迭代法最终解 x = [' );
% fprintf('%.16f, ', x_sor(1:end-1));
% fprintf('%.16f]^T\n', x_sor(end));
% ============================================ 可视化模块 ============================================
visualizeResults(x_gauss, x_jacobi, x_gs, x_sor, x_exact, ...
err_history_jacobi, err_history_gs, err_history_sor, ...
iter_jacobi, iter_gs, iter_sor, ...
spectral_radius_jacobi, spectral_radius_gs, spectral_radius_sor, epsilon);
fprintf('========================================================================================================\n');
fprintf('程序运行结束。\n');
fprintf('========================================================================================================\n');
end
% ====================== 生成 Hilbert 矩阵 ======================
function H = myHilbert(n)
H = zeros(n);
for i = 1:n
for j = 1:n
H(i,j) = 1 / (i + j - 1);
end
end
end
% ====================== 高斯消去法 ======================
function x = gaussianElimination(A, b)
fprintf('=================== 高斯消去法求解过程 ===================\n');
n = length(b);
% 构造增广矩阵
Aug = [A, b];
fprintf('初始增广矩阵 [A|b]:\n');
disp(Aug);
% 前向消元
fprintf('\n--- 前向消元过程 ---\n');
for k = 1:n-1
fprintf('第 %d 步消元:\n', k);
fprintf('主元: A(%d,%d) = %.6f\n', k, k, Aug(k,k));
for i = k+1:n
factor = Aug(i,k)/Aug(k,k);
fprintf('对第 %d 行，消元因子 = %.6f\n', i, factor);
Aug(i,k:n+1) = Aug(i,k:n+1) - factor*Aug(k,k:n+1);
fprintf('消元后的增广矩阵:\n');
disp(Aug);
end
fprintf('\n');
end
fprintf('消元完成后的上三角增广矩阵:\n');
disp(Aug);
% 回代求解
fprintf('\n--- 回代求解过程 ---\n');
x = zeros(n,1);
for i = n:-1:1
if i == n
x(i) = Aug(i,n+1)/Aug(i,i);
fprintf('x(%d) = b(%d)/A(%d,%d) = %.6f/%.6f = %.6f\n', ...
i, i, i, i, Aug(i,n+1), Aug(i,i), x(i));
else
sum_val = Aug(i,i+1:n)*x(i+1:n);
x(i) = (Aug(i,n+1) - sum_val)/Aug(i,i);
fprintf('x(%d) = (b(%d) - A(%d,%d:%d)*x(%d:%d))/A(%d,%d) = (%.6f - %.6f)/%.6f = %.6f\n', ...
i, i, i, i+1, n, i+1, n, i, i, Aug(i,n+1), sum_val, Aug(i,i), x(i));
end
end
fprintf('\n高斯消去法求解完成\n');
fprintf('==========================================================\n');
end
% ====================== 使用增广矩阵求解线性方程组 ======================
function x = solveWithAugmented(A, b)
% 使用增广矩阵和高斯消元法求解 Ax = b
n = length(b);
Aug = [A, b];
% 前向消元
for k = 1:n-1
for i = k+1:n
factor = Aug(i,k)/Aug(k,k);
Aug(i,k:n+1) = Aug(i,k:n+1) - factor*Aug(k,k:n+1);
end
end
% 回代求解
x = zeros(n,1);
for i = n:-1:1
x(i) = (Aug(i,n+1) - Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
end
end
% ====================== 幂方法计算谱半径 ======================
function spectral_radius = powerMethod(M, maxIter, tol)
% 幂方法计算矩阵的谱半径
n = size(M, 1);
x = rand(n, 1);
x = x / norm(x);
lambda_prev = 0;
for k = 1:maxIter
y = M * x;
lambda = norm(y);
x = y / lambda;
if abs(lambda - lambda_prev) < tol
break;
end
lambda_prev = lambda;
end
spectral_radius = lambda;
end
% ====================== 手动计算2x2矩阵特征值 ======================
function eigenvalues = manualEigenvalues2x2(M)
% 手动计算2x2矩阵的特征值
a = M(1,1); b = M(1,2);
c = M(2,1); d = M(2,2);
trace_M = a + d;
det_M = a*d - b*c;
discriminant = trace_M^2 - 4*det_M;
if discriminant >= 0
lambda1 = (trace_M + sqrt(discriminant)) / 2;
lambda2 = (trace_M - sqrt(discriminant)) / 2;
eigenvalues = [lambda1, lambda2];
else
real_part = trace_M / 2;
imag_part = sqrt(-discriminant) / 2;
eigenvalues = [real_part + 1i*imag_part, real_part - 1i*imag_part];
end
end
% ====================== Jacobi 迭代 ======================
function [x, err_history, iter, spectral_radius] = jacobiIteration(A, b, maxIter, epsilon, x_exact)
n = length(b);
x = zeros(n,1);
err_history = zeros(maxIter, 1);
D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,1);
fprintf('---------------- Jacobi 迭代 ----------------\n');
fprintf('系数矩阵 A:\n'); disp(A);
fprintf('对角矩阵 D:\n'); disp(D);
fprintf('严格下三角矩阵 -L:\n'); disp(-L);
fprintf('严格上三角矩阵 -U:\n'); disp(U);
% 计算迭代矩阵 M = D^{-1}(L+U)
M = zeros(n);
for i = 1:n
% 计算 M 的第 i 行：求解 D * m_i = (L+U) 的第 i 行
rhs = L(i,:)' + U(i,:)';
m_i = solveWithAugmented(D, rhs);
M(i,:) = m_i';
end
% 计算常数向量 c = D^{-1}b
c = solveWithAugmented(D, b);
% 计算谱半径
if n == 2
eigenvalues = manualEigenvalues2x2(M);
spectral_radius = max(abs(eigenvalues));
else
spectral_radius = powerMethod(M, 1000, 1e-10);
end
fprintf('Jacobi 迭代矩阵 M = D^{-1}(L+U):\n');
for i = 1:n
fprintf('[');
for j = 1:n
fprintf('%.8f', M(i,j));
if j < n, fprintf(', '); end
end
fprintf(']\n');
end
fprintf('Jacobi 常向量 c = D^{-1}b = [' );
fprintf('%.16f, ', c(1:end-1));
fprintf('%.16f]^T\n', c(end));
fprintf('Jacobi 迭代矩阵谱半径 = %.6f\n', spectral_radius);
for k = 1:maxIter
% Jacobi迭代: x_new = D^{-1}(b + (L+U)x)
rhs = b + (L+U)*x;
x_new = solveWithAugmented(D, rhs);
err = norm(x_new - x_exact);
err_history(k) = err;
fprintf('Jacobi 第 %d 次迭代，绝对误差 = %.6e\n', k, err);
fprintf('Jacobi 迭代法当前解 x = [' );
fprintf('%.16f, ', x_new(1:end-1));
fprintf('%.16f]^T\n', x_new(end));
x = x_new;
if err < epsilon
fprintf('收敛达到绝对误差 epsilon = %.1e，在第 %d 次迭代完成。\n', epsilon, k);
err_history = err_history(1:k);
iter = k;
return;
end
end
fprintf('达到最大迭代次数 %d，最终绝对误差 = %.6e\n', maxIter, err);
iter = maxIter;
end
% ====================== Gauss-Seidel 迭代 ======================
function [x, err_history, iter, spectral_radius] = gaussSeidelIteration(A, b, maxIter, epsilon, x_exact)
n = length(b);
x = zeros(n,1);
err_history = zeros(maxIter, 1);
D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,1);
fprintf('---------------- Gauss-Seidel 迭代 ----------------\n');
fprintf('系数矩阵 A:\n'); disp(A);
fprintf('对角矩阵 D:\n'); disp(D);
fprintf('严格下三角矩阵 -L:\n'); disp(-L);
fprintf('严格上三角矩阵 -U:\n'); disp(U);
% 计算迭代矩阵 M = (D-L)^{-1}U
M = zeros(n);
DL = D - L;
% 按列计算 M
for j = 1:n
% 计算 M 的第 j 列：求解 (D-L) * m_j = U 的第 j 列
rhs = U(:, j);
m_j = solveWithAugmented(DL, rhs);
M(:, j) = m_j;
end
% 计算常数向量 c = (D-L)^{-1}b
c = solveWithAugmented(D-L, b);
% 计算谱半径
if n == 2
eigenvalues = manualEigenvalues2x2(M);
spectral_radius = max(abs(eigenvalues));
% 详细输出2x2矩阵的计算过程
fprintf('\n--- Gauss-Seidel 谱半径详细计算 ---\n');
fprintf('迭代矩阵 M:\n');
for i = 1:n
fprintf('[');
for j = 1:n
fprintf('%.8f', M(i,j));
if j < n, fprintf(', '); end
end
fprintf(']\n');
end
fprintf('特征值: %.6f, %.6f\n', eigenvalues(1), eigenvalues(2));
else
spectral_radius = powerMethod(M, 1000, 1e-10);
end
fprintf('Gauss-Seidel 迭代矩阵谱半径 = %.6f\n', spectral_radius);
fprintf('Gauss-Seidel 常向量 c = (D-L)^{-1}b = [' );
fprintf('%.16f, ', c(1:end-1));
fprintf('%.16f]^T\n', c(end));
for k = 1:maxIter
% Gauss-Seidel迭代: (D-L) * x_new = b + U*x
rhs = b + U*x;
x_new = solveWithAugmented(D-L, rhs);
err = norm(x_new - x_exact);
err_history(k) = err;
fprintf('Gauss-Seidel 第 %d 次迭代，绝对误差 = %.6e\n', k, err);
fprintf('Gauss-Seidel 迭代法当前解 x = [' );
fprintf('%.16f, ', x_new(1:end-1));
fprintf('%.16f]^T\n', x_new(end));
x = x_new;
if err < epsilon
fprintf('收敛达到绝对误差 epsilon = %.1e，在第 %d 次迭代完成。\n', epsilon, k);
err_history = err_history(1:k);
iter = k;
return;
end
end
fprintf('达到最大迭代次数 %d，最终绝对误差 = %.6e\n', maxIter, err);
iter = maxIter;
end
% ====================== SOR 迭代 ======================
function [x, err_history, iter, spectral_radius] = SORIteration(A, b, omega, maxIter, epsilon, x_exact)
n = length(b);
x = zeros(n,1);
err_history = zeros(maxIter, 1);
D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,1);
fprintf('---------------- SOR 迭代 ----------------\n');
fprintf('系数矩阵 A:\n'); disp(A);
fprintf('对角矩阵 D:\n'); disp(D);
fprintf('严格下三角矩阵 -L:\n'); disp(-L);
fprintf('严格上三角矩阵 -U:\n'); disp(U);
fprintf('松弛因子 ω = %.2f\n', omega);
% 计算迭代矩阵 M = (D-ωL)^{-1}[(1-ω)D + ωU]
M = zeros(n);
D_omegaL = D - omega*L;
omegaU_plus = (1-omega)*D + omega*U;
% 按列计算 M
for j = 1:n
% 计算 M 的第 j 列：求解 (D-ωL) * m_j = ((1-ω)D + ωU) 的第 j 列
rhs = omegaU_plus(:, j);
m_j = solveWithAugmented(D_omegaL, rhs);
M(:, j) = m_j;
end
% 计算常数向量 c = ω*(D-ωL)^{-1}b
c = omega * solveWithAugmented(D-omega*L, b);
% 计算谱半径
if n == 2
eigenvalues = manualEigenvalues2x2(M);
spectral_radius = max(abs(eigenvalues));
else
spectral_radius = powerMethod(M, 1000, 1e-10);
end
fprintf('SOR 迭代矩阵谱半径 = %.6f\n', spectral_radius);
fprintf('SOR 常向量 c = ω*(D-ωL)^{-1}b = [' );
fprintf('%.16f, ', c(1:end-1));
fprintf('%.16f]^T\n', c(end));
for k = 1:maxIter
% SOR迭代: (D-ωL) * x_new = ω*b + ((1-ω)D + ωU)*x
rhs = omega*b + ((1-omega)*D + omega*U)*x;
x_new = solveWithAugmented(D-omega*L, rhs);
err = norm(x_new - x_exact);
err_history(k) = err;
fprintf('SOR 第 %d 次迭代，绝对误差 = %.6e\n', k, err);
fprintf('SOR 迭代法当前解 x = [' );
fprintf('%.16f, ', x_new(1:end-1));
fprintf('%.16f]^T\n', x_new(end));
x = x_new;
if err < epsilon
fprintf('收敛达到绝对误差 epsilon = %.1e，在第 %d 次迭代完成。\n', epsilon, k);
err_history = err_history(1:k);
iter = k;
return;
end
end
fprintf('达到最大迭代次数 %d，最终绝对误差 = %.6e\n', maxIter, err);
iter = maxIter;
end
% ====================== 可视化结果 ======================
function visualizeResults(x_gauss, x_jacobi, x_gs, x_sor, x_exact, ...
err_history_jacobi, err_history_gs, err_history_sor, ...
iter_jacobi, iter_gs, iter_sor, ...
spectral_radius_jacobi, spectral_radius_gs, spectral_radius_sor, epsilon)
% 创建一张图，包含两个子图
figure('Position', [100, 100, 1200, 500]);
% 子图1: 绝对误差收敛曲线（对数坐标）
subplot(1, 2, 1);
iter_range_jacobi = 1:length(err_history_jacobi);
iter_range_gs = 1:length(err_history_gs);
iter_range_sor = 1:length(err_history_sor);
% 使用semilogy绘制对数坐标图
semilogy(iter_range_jacobi, err_history_jacobi, 'r^-', 'LineWidth', 1.5, 'MarkerSize', 2, 'DisplayName', 'Jacobi');
hold on;
semilogy(iter_range_gs, err_history_gs, 'gs-', 'LineWidth', 1.5, 'MarkerSize', 2, 'DisplayName', 'Gauss-Seidel');
semilogy(iter_range_sor, err_history_sor, 'bd-', 'LineWidth', 1.5, 'MarkerSize', 2, 'DisplayName', 'SOR');
% 兼容性方法添加收敛阈值线（替代yline）
x_limits = xlim;
plot(x_limits, [epsilon, epsilon], 'k--', 'LineWidth', 1.5, 'DisplayName', sprintf('收敛阈值 (%.0e)', epsilon));
xlabel('迭代次数', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('绝对误差', 'FontSize', 12, 'FontWeight', 'bold');
title('迭代法绝对误差收敛曲线', 'FontSize', 14, 'FontWeight', 'bold');
% 设置y轴范围从1e-5到1e+0
ylim([epsilon, 1e+2]);
% 设置y轴刻度为1e+0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5
yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e+0, 1e+1, 1e+2]);
% 设置y轴刻度标签格式
yticklabels({'10^{-5}', '10^{-4}', '10^{-3}', '10^{-2}', '10^{-1}', '10^{0}', '10^{1}', '10^{2}'});
grid on;
set(gca, 'FontSize', 11);
% 兼容性修改：先创建legend，然后设置字体大小
h_legend = legend('Location', 'best');
try
set(h_legend, 'FontSize', 11);
catch
% 如果设置字体大小失败，忽略错误
end
% 子图2: 迭代次数比较
subplot(1, 2, 2);
iter_counts = [iter_jacobi, iter_gs, iter_sor];
methods_iter = {'Jacobi', 'Gauss-Seidel', 'SOR'};
bar(iter_counts, 'FaceColor', [0.2 0.6 0.8], 'EdgeColor', 'k', 'LineWidth', 1.5);
set(gca, 'XTickLabel', methods_iter, 'FontSize', 12, 'FontWeight', 'bold');
ylabel('迭代次数', 'FontSize', 12, 'FontWeight', 'bold');
title('达到收敛所需的迭代次数', 'FontSize', 14, 'FontWeight', 'bold');
grid on;
% 添加数值标签
for i = 1:length(iter_counts)
text(i, iter_counts(i), sprintf('%d', iter_counts(i)), ...
'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', ...
'FontSize', 14, 'FontWeight', 'bold', 'Color', 'red');
end
% 设置y轴范围，使数值标签显示更清晰
y_max = max(iter_counts) * 1.5;
ylim([0, y_max]);
% 添加总标题
try
sgtitle('Hilbert矩阵方程求解方法性能比较', 'FontSize', 16, 'FontWeight', 'bold');
catch
set(gcf, 'Name', 'Hilbert矩阵方程求解方法性能比较');
end
% 输出对齐的性能统计
fprintf('\n=============== 性能统计 ===============\n');
fprintf('方法\t\t\t最终误差\t\t迭代次数\t谱半径\t\t收敛状态\n');
fprintf('--------\t\t--------\t\t--------\t--------\t--------\n');
fprintf('高斯消去法\t\t%.6e\t\tN/A\t\t\tN/A\t\t\t直接法\n', norm(x_gauss - x_exact));
fprintf('Jacobi迭代\t\t%.6e\t\t%d\t\t\t%.16f\t\t%s\n', norm(x_jacobi - x_exact), iter_jacobi, ...
spectral_radius_jacobi, getConvergenceStatus(err_history_jacobi(end), epsilon));
fprintf('Gauss-Seidel\t%.6e\t\t%d\t\t\t%.16f\t\t%s\n', norm(x_gs - x_exact), iter_gs, ...
spectral_radius_gs, getConvergenceStatus(err_history_gs(end), epsilon));
fprintf('SOR迭代\t\t\t%.6e\t\t%d\t\t\t%.16f\t\t%s\n', norm(x_sor - x_exact), iter_sor, ...
spectral_radius_sor, getConvergenceStatus(err_history_sor(end), epsilon));
end
% ====================== 辅助函数：获取收敛状态 ======================
function status = getConvergenceStatus(final_error, epsilon)
if final_error < epsilon
status = '已收敛';
else
status = '未收敛';
end
end
% ====================== SOR松弛因子优化函数 ======================
function findOptimalOmega(n, maxIter, epsilon, current_omega)
% SOR松弛因子优化分析
% 输入：n - Hilbert矩阵维数，maxIter - 最大迭代次数
fprintf('开始SOR松弛因子优化分析...\n');
% 生成Hilbert矩阵和精确解
H = myHilbert(n);
x_exact = ones(n, 1);
b = H * x_exact;
% 定义测试参数 - 使用与主程序相同的精度要求
omega_range = 0.1:0.05:1.9; % 松弛因子范围
epsilon_target = epsilon; % 修正：使用与主程序相同的精度要求
% 预存储结果
iter_results = zeros(length(omega_range), 1);
spectral_radius_results = zeros(length(omega_range), 1);
final_error_results = zeros(length(omega_range), 1);
fprintf('测试松弛因子范围: %.1f ~ %.1f，步长 %.2f\n', min(omega_range), max(omega_range), omega_range(2)-omega_range(1));
fprintf('测试精度要求: %.0e\n', epsilon_target);
% 对每个松弛因子进行测试
for i = 1:length(omega_range)
omega = omega_range(i);
% 运行SOR迭代
[x_sor, ~, iter_count, spectral_radius] = SORIteration(H, b, omega, maxIter, epsilon_target, x_exact);
iter_results(i) = iter_count;
spectral_radius_results(i) = spectral_radius;
final_error_results(i) = norm(x_sor - x_exact);
if mod(i, 5) == 0
fprintf('已完成 %.0f%% 的测试...\n', i/length(omega_range)*100);
end
end
% 可视化结果
figure('Position', [200, 200, 1400, 600]);
% 子图1: 迭代次数 vs 松弛因子
subplot(1, 1, 1);
plot(omega_range, iter_results, 'bo-', 'LineWidth', 2, 'MarkerSize', 2);
xlabel('松弛因子 ω', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('迭代次数', 'FontSize', 12, 'FontWeight', 'bold');
title(sprintf('SOR方法：迭代次数 vs 松弛因子 (ε=%.0e)', epsilon_target), 'FontSize', 14, 'FontWeight', 'bold');
grid on;
% 标记最优松弛因子（迭代次数最少）
[min_iter, min_iter_idx] = min(iter_results);
optimal_omega_iter = omega_range(min_iter_idx);
hold on;
plot(optimal_omega_iter, min_iter, 'rs', 'MarkerSize', 8, 'LineWidth', 2);
% 标记主程序使用的松弛因子
main_omega = current_omega;
main_omega_idx = find(abs(omega_range - main_omega) < 0.01);
if ~isempty(main_omega_idx)
main_iter = iter_results(main_omega_idx(1));
plot(main_omega, main_iter, 'g^', 'MarkerSize', 8, 'LineWidth', 2);
legend_str = {sprintf('迭代次数'), ...
sprintf('最优 ω=%.2f', optimal_omega_iter), ...
sprintf('主程序 ω=%.2f', main_omega)};
else
legend_str = {sprintf('迭代次数'), ...
sprintf('最优 ω=%.2f', optimal_omega_iter)};
end
legend(legend_str, 'Location', 'best', 'FontSize', 10);
fprintf('精度 ε=%.0e: 最优ω=%.2f, 最少迭代次数=%d\n', epsilon_target, optimal_omega_iter, min_iter);
fprintf('主程序使用的ω=%.2f: 迭代次数=%d\n', main_omega, main_iter);
% 添加总标题
try
sgtitle(sprintf('Hilbert(%d)矩阵SOR方法松弛因子优化分析', n), ...
'FontSize', 16, 'FontWeight', 'bold');
catch
set(gcf, 'Name', sprintf('Hilbert(%d)矩阵SOR方法松弛因子优化分析', n));
end
% 输出详细分析结果
fprintf('\n=============== SOR松弛因子优化结果汇总 ===============\n');
fprintf('矩阵维数: %d\n', n);
fprintf('测试的松弛因子范围: %.1f ~ %.1f\n', min(omega_range), max(omega_range));
fprintf('精度要求: %.0e\n', epsilon_target);
fprintf('最优松弛因子 (最少迭代次数): ω = %.2f\n', optimal_omega_iter);
fprintf('最少迭代次数: %d\n', min_iter);
fprintf('主程序使用的松弛因子: ω = %.2f\n', main_omega);
fprintf('主程序对应的迭代次数: %d\n', main_iter);
% 提供建议
fprintf('\n建议:\n');
if optimal_omega_iter ~= main_omega
fprintf('建议使用 ω = %.2f 代替当前的 ω = %.2f，可以减少迭代次数从 %d 到 %d\n', ...
optimal_omega_iter, main_omega, main_iter, min_iter);
else
fprintf('当前使用的松弛因子 ω = %.2f 是最优选择\n', main_omega);
end
fprintf('========================================================\n');
end