共轭梯度法 MATLAB函数代码

function x = My_CGM(x0,A,b)
% 共轭梯度法(Conjugate Gradient Method)求解求解线性方程 Ax = b
%SYNTAX:
% x = My_CGM(x0,A,b,epsilon)
%DESCRIPTION:
%
% INPUT:
% x0 -- 初始解
% A -- 系数矩阵
% b -- 线性方程右端列向量
% OUTPUT:
% x -- 解列向量
%REFERENCES:
% 《有限单元法》 by 王勖成 P246
%REVISION: 1.0

b = b(:);
x0 = x0(:);

x_old = x0;
gamma_old = b - A*x_old;
p_old = gamma_old;
N = length(b);
for ii = 1:N
alpha = (p_old.'*gamma_old)/(p_old.'*A*p_old);% Eq.6.5.13
x_new = x_old + alpha*p_old;

gamma_new = b - A*x_new;
beta = -(p_old.'*A*p_old)/(p_old.'*A*gamma_new);% Eq.6.5.13
p_new = p_old + beta*gamma_new;

% update
x_old = x_new;
p_old = p_new;

end

x = x_new;
end