共轭梯度matlab代码

下面是一个简单的共轭梯度算法的 MATLAB 代码示例：

function [x, flag, relres, iter] = cg(A, b, x0, tol, maxiter)
% Conjugate Gradient Algorithm
% Solves the system of linear equations Ax = b
% A: matrix of size nxn
% b: column vector of size nx1
% x0: initial guess for x
% tol: tolerance for stopping criterion
% maxiter: maximum number of iterations

n = size(A, 1);
x = x0;
r = b - A * x;
p = r;
normb = norm(b);

for iter = 1:maxiter
    Ap = A * p;
    alpha = (r' * r) / (p' * Ap);
    x = x + alpha * p;
    rprev = r;
    r = r - alpha * Ap;
    beta = (r' * r) / (rprev' * rprev);
    p = r + beta * p;
    relres = norm(r) / normb;
    if relres < tol
        flag = 0;
        return;
    end
end

flag = 1; % If we reach here, we have not converged within maxiter iterations
end

其中，`A` 是矩阵，`b` 是列向量，`x0` 是初始值，`tol` 是停止准则，`maxiter` 是最大迭代次数。函数返回 `x`，`flag`，`relres`，`iter`。`flag` 是标志变量，如果 `flag` 等于 0，则表示算法已经收敛，否则表示算法未能在最大迭代次数内收敛。`relres` 是相对残差，`iter` 是实际迭代次数。

使用示例：

% Define matrix and vector
n = 100;
A = gallery('tridiag', n, -1, 2, -1);
b = ones(n, 1);

% Solve the system using CG
x0 = zeros(n, 1);
tol = 1e-6;
maxiter = 1000;
[x, flag, relres, iter] = cg(A, b, x0, tol, maxiter);

% Print results
if flag == 0
    fprintf('CG converged in %d iterations.\n', iter);
    fprintf('Relative residual: %e\n', relres);
else
    fprintf('CG did not converge in %d iterations.\n', maxiter);
end

其中，`gallery` 函数用于创建三对角矩阵。