matlab中实现共轭梯度法

实现共轭梯度法的Matlab代码如下：

function[x, flag, relres, iter, resvec] = cg(A, b, x0, tol, maxit, M)

% A: n×n 矩阵

% b: n维向量

% x0: 初始值

% tol: 相对误差容限

% maxit: 迭代次数上限

% M: 预处理矩阵

% x: 解向量

% flag: 表示迭代收敛性的标志

% relres: 相对残差

% iter: 实际迭代次数

% resvec: 历史残差

if nargin < 3 || isempty(x0)

x0 = zeros(size(b)); % 默认初始值为0向量

end

if nargin < 4 || isempty(tol)

tol = 1e-6; % 默认相对误差容限为1e-6

end

if nargin < 5 || isempty(maxit)

maxit = min(length(b), 20); % 默认迭代次数上限为n或20中的较小值

end

if nargin < 6 || isempty(M)

M = []; % 默认不使用预处理

end

if ~all(size(A) == size(b))

error('矩阵A与向量b的维度不匹配！');

end

if max(abs(diag(A))) == 0

error('矩阵A的对角线元素不能全为0！');

end

flag = 1; % 初始化收敛性标志为1

x = x0; % 初始化解向量

r = b - A*x; % 初始化残差向量

if ~isempty(M)

z = M \ r; % 使用预处理求解

else

z = r;

end

p = z; % 初始化搜索方向向量

rho = z'*r; % 初始化rho

resvec = zeros(maxit+1, 1); % 存储历史残差

resvec(1) = norm(r);

for iter = 1:maxit % 开始迭代

q = A*p;

alpha = rho/(p'*q);

x = x + alpha*p;

r = r - alpha*q;

if ~isempty(M)

z = M \ r; % 使用预处理求解

else

z = r;

end

rho_new = z'*r;

beta = rho_new/rho;

rho = rho_new;

p = z + beta*p;

resvec(iter+1) = norm(r);

if norm(r) <= tol*norm(b)

flag = 0; % 迭代收敛

break;

end

end

resvec = resvec(1:iter+1);

relres = norm(r)/(norm(A)*norm(x) + norm(b));

if flag

warning('共轭梯度法未收敛！');

end

end

希望这个答案可以帮到您，如有疑问请随时追问！