matlab实现CG的代码
时间: 2023-10-13 20:06:41 浏览: 63
以下是使用 MATLAB 实现共轭梯度法(CG)的代码示例:
```matlab
function [x, flag, relres, iter, resvec] = cg(A, b, x0, maxit, tol)
% CG Conjugate Gradient Method for solving Ax = b
% x = CG(A, b) attempts to solve the system of linear equations A*x = b
% for x using the Conjugate Gradient method.
% A can be a matrix or a function handle that returns A*x.
% b is the right-hand side vector.
% x0 is the initial guess (default is zero vector).
% maxit is the maximum number of iterations (default is 100).
% tol is the tolerance for convergence (default is 1e-6).
%
% [x, flag] = CG(A, b) also returns a convergence flag:
% 0 - CG converged to the desired tolerance within maxit iterations.
% 1 - CG iterated maxit times but did not converge.
%
% [x, flag, relres] = CG(A, b) also returns the relative residual norm(b-A*x)/norm(b).
%
% [x, flag, relres, iter] = CG(A, b) also returns the number of iterations.
%
% [x, flag, relres, iter, resvec] = CG(A, b) also returns the residual history.
%
% Example:
% A = gallery('poisson', 100);
% b = randn(size(A, 1), 1);
% x = CG(A, b);
% norm(A*x-b)/norm(b)
%
% See also PCG, GMRES
% Reference: Nocedal and Wright, Numerical Optimization, 2nd ed., 2006, p. 86.
if nargin < 5 || isempty(tol), tol = 1e-6; end
if nargin < 4 || isempty(maxit), maxit = 100; end
if nargin < 3 || isempty(x0), x0 = zeros(size(b)); end
% Initialize
r = b - A(x0);
p = r;
x = x0;
resvec = zeros(maxit+1, 1);
resvec(1) = norm(r);
flag = 0;
for iter = 1:maxit
Ap = A(p);
alpha = (r'*r) / (p'*Ap);
x = x + alpha*p;
r_old = r;
r = r_old - alpha*Ap;
beta = (r'*r) / (r_old'*r_old);
p = r + beta*p;
resvec(iter+1) = norm(r);
relres = resvec(iter+1) / resvec(1);
if relres < tol
flag = 0;
break;
end
end
if flag == 0
% Converged
resvec = resvec(1:iter+1);
else
% Not converged
resvec = resvec(1:maxit+1);
end
end
```
使用示例:
```matlab
% Generate a random matrix A and a random vector b
n = 100;
A = randn(n);
A = A'*A;
b = randn(n, 1);
% Solve Ax = b using CG
[x, flag, relres, iter, resvec] = cg(A, b);
% Plot the residual history
semilogy(resvec);
xlabel('Iteration');
ylabel('Residual norm');
title('CG');
```