MATLAB写出共轭梯度法解线性方程组

以下是MATLAB代码实现共轭梯度法解线性方程组：

function [x,flag,relres,iter,resvec] = cg(A,b,tol,maxit,M)
% 共轭梯度法解线性方程组 Ax=b
% 输入参数：
% A: 系数矩阵
% b: 右端向量
% tol: 相对误差容限
% maxit: 最大迭代次数
% M: 预处理矩阵或者函数句柄
% 输出参数：
% x: 迭代得到的解向量
% flag: 表示收敛性的标志
% relres: 相对误差
% iter: 实际迭代次数
% resvec: 残量的迭代历史记录

n = length(b);
x = zeros(n,1);
r = b - A*x;
if nargin < 5
    M = eye(n);
end
z = M \ r;
p = z;
rho = z'*r;
resvec = zeros(maxit+1,1);
resvec(1) = norm(r);
tol = tol * resvec(1);
for iter = 1:maxit
    q = A*p;
    alpha = rho / (p'*q);
    x = x + alpha*p;
    r = r - alpha*q;
    resvec(iter+1) = norm(r);
    if resvec(iter+1) < tol
        break;
    end
    z = M \ r;
    rho1 = rho;
    rho = z'*r;
    beta = rho / rho1;
    p = z + beta*p;
end
if resvec(iter+1) < tol
    flag = 0;
else
    flag = 1;
end
relres = resvec(iter+1) / resvec(1);
resvec(iter+2:end) = [];
end

其中，输入参数 `A` 是系数矩阵，`b` 是右端向量，`tol` 是相对误差容限，`maxit` 是最大迭代次数，`M` 是预处理矩阵或函数句柄。输出参数 `x` 是迭代得到的解向量，`flag` 表示收敛性的标志，`relres` 是相对误差，`iter` 是实际迭代次数，`resvec` 是残量的迭代历史记录。