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

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

1. 引言

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

2. 数学原理

2.1 问题描述

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

2.2 目标函数

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

$f(x)=\frac{1}{2}{x}^{T}Ax-b^{T}x$

其梯度为：

$\nabla f(x)=Ax-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\text{for}i\ne j$.

2.5 共轭梯度法的迭代公式

1. 初始化：选择初始点$x_0$，计算 $r_0=b-A{x}_0$，并令 $p_0=r_0$.
2. 迭代步骤：
   • 计算步长${\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$.
3. 检查收敛：如果$r_{k+1}$ 足够小，则停止迭代。

3. 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:maxIter
    Ap = 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) < tol
        res_norms = res_norms(1:iter+1);
        break;
    end

    p = r + (rsnew / rsold) * p;  % 更新方向向量
    rsold = rsnew;
end

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

% 测试案例
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');

画出收敛曲线：

残差设定 $1e-8$量级，几秒就能算出结果，效果不错！

续

如果读者有需求，我们将通过一系列博客展示数值分析课程相关的丰富内容，所有文章均有相应代码实现。请持续关注！

---

作者 ：计算小屋
个人主页 ： 计算小屋的主页