MATLAB 共轭梯度法求解线性方程组（附代码）

共轭梯度法求解线性方程组

在这里插入图片描述

1. 引言

共轭梯度法（Conjugate Gradient Method）是一种用于求解大型稀疏对称正定线性方程组的迭代算法。该方法结合了梯度下降法和共轭方向的概念，以达到更快速的收敛。共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，在有限元求解中经常用到此方法求解椭圆问题。

2. 数学原理

2.1 问题描述

给定一个对称正定矩阵 ( A ) 和一个向量 ( b )，我们需要求解线性方程组 ( Ax = b )。

2.2 目标函数

我们可以通过最小化二次型函数来求解上述线性方程组：

$\frac{1}{2} x^T A x - b^T x$

其梯度为：

$\nabla f(x) = A x - b$

2.3 梯度下降法

梯度下降法的更新公式为：

$x_{k+1} = x_k - \alpha_k \nabla f(x_k)$

其中 $\alpha_k$ 是步长。

2.4 共轭方向

共轭方向是一组特殊的搜索方向，满足以下条件：

$p_i^T A p_j = 0 \quad \text{for} \quad i \neq j$ .

2.5 共轭梯度法的迭代公式

初始化：选择初始点 $x_0$ ，计算 $r_0 = b - A x_0$ ，并令 $p_0 = r_0$ .
迭代步骤：
- 计算步长 $\alpha_k$ ：
  $\alpha_k = \frac{r_k^T r_k}{p_k^T A p_k}$ .
- 更新解向量 $x_{k+1}$ ：
  $x_{k+1} = x_k + \alpha_k p_k$ .
- 更新残差 $r_{k+1}$ ：
  $r_{k+1} = r_k - \alpha_k A p_k$ .
- 计算方向更新系数 $\beta_k$ ：
  $\beta_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}$ .
- 更新搜索方向 $p_{k+1}$ ：
  $p_{k+1} = r_{k+1} + \beta_k p_k$ .
检查收敛：如果 $r_{k+1}$ 足够小，则停止迭代。

3. MATLAB 实现

以下是详细的 MATLAB 程序实现共轭梯度法。

matlab">function [u,iter,res_norms] = conjugate_gradient(A, F, tol, maxIter)
%% 使用共轭梯度法求解 A * u = F
% A - 系数矩阵
% F - 右端向量
% tol - 收敛容差
% maxIter - 最大迭代次数% 初始化解向量
u = zeros(size(F));
r = F - A * u;  % 初始残差
p = r;  % 初始方向向量
rsold = r' * r;  % 初始残差范数res_norms = zeros(maxIter, 1);res_norms(1) = sqrt(rsold);for iter = 1:maxIterAp = A * p;alpha = rsold / (p' * Ap);u = u + alpha * p;  % 更新解r = r - alpha * Ap;  % 更新残差rsnew = r' * r;  % 新的残差范数res_norms(iter+1) = sqrt(rsnew);% 检查收敛性if sqrt(rsnew) < tolres_norms = res_norms(1:iter+1);break;endp = r + (rsnew / rsold) * p;  % 更新方向向量rsold = rsnew;
end

搞个 250000阶的大矩阵测试一下：

matlab">% 测试案例
n = 250000; % 矩阵规模
A = gallery('poisson', sqrt(n)); % 生成对称正定矩阵
b = rand(n, 1); % 生成右侧向量
% x0 = zeros(n, 1); % 初始解
tol = 1e-8;
max_iter = 100000;% 调用共轭梯度法
[x, k, res_norms] = conjugate_gradient(A, b,  tol, max_iter);% 显示结果
fprintf('迭代次数 k: %d\n', k);% 可视化结果
figure;
semilogy(res_norms, 'o-');
xlabel('迭代次数');
ylabel('残差范数');
title('共轭梯度法残差范数收敛曲线');
grid on;
set(gca,'FontName','Monospaced');