Методы сопряженных градиентов придали методу квадратную форму
x = cgs (A, b)
cgs (A, b, tol)
cgs (A, b, tol, maxit)
cgs (A, b, tol, maxit, M)
cgs (A, b, tol, maxit, M1, M2)
cgs (A, b, tol, maxit, M1, M2, x0)
[x, флаг] = cgs (A, b...)
[x, флаг, relres] = cgs (A, b...)
[x, флаг, relres, проход] = cgs (A, b...)
[x, флаг, relres, проход, resvec] = cgs (A, b...)
x = cgs(A,b) пытается решить систему линейных уравнений A*x = b для x. n-by-n матрица коэффициентов A должен быть квадратным и должен быть большим и разреженным. Вектор - столбец b должен иметь длину n. Можно задать A как указатель на функцию, afun, такой, что afun(x) возвращает A*x.
Параметризация Функций объясняет, как предоставить дополнительные параметры функциональному afun, а также функцию перед формирователем mfun, описанный ниже, при необходимости.
Если cgs сходится, сообщение к тому эффекту отображено. Если cgs не удается сходиться после максимального количества итераций или остановов по любой причине, предупреждающее сообщение распечатано, отобразив относительный остаточный norm(b-A*x)/norm(b) и номер итерации в который метод, остановленный или не пройдено.
cgs(A,b,tol) задает допуск метода, tol. Если tol является [], то cgs использует значение по умолчанию, 1e-6.
cgs(A,b,tol,maxit) задает максимальное количество итераций, maxit. Если maxit является [] затем, cgs использует значение по умолчанию, min(n,20).
cgs(A,b,tol,maxit,M) и cgs(A,b,tol,maxit,M1,M2) используют предварительный формирователь M или M = M1*M2 и эффективно решают систему inv(M)*A*x = inv(M)*b для x. Если M является [] затем, cgs не применяет предварительного формирователя. M может быть указателем на функцию mfun, таким образом, что mfun(x) возвращает M\x.
cgs(A,b,tol,maxit,M1,M2,x0) задает исходное предположение x0. Если x0 является [], то cgs использует значение по умолчанию, все-нулевой вектор.
[x,flag] = cgs(A,b,...) возвращает решение x и флаг, который описывает сходимость cgs.
Флаг | Сходимость |
|---|---|
|
|
|
|
| Предварительный формирователь |
|
|
| Один из скаляров, вычисленных во время |
Каждый раз, когда flag не является 0, решение, возвращенный x то, что с минимальной невязкой нормы, вычисленной по всем итерациям. Никакие сообщения не отображены, если flag вывод задан.
[x,flag,relres] = cgs(A,b,...) также возвращает относительный остаточный norm(b-A*x)/norm(b). Если flag является 0, то relres <= tol.
[x,flag,relres,iter] = cgs(A,b,...) также возвращает номер итерации, в котором x был вычислен, где 0 <= iter <= maxit.
[x,flag,relres,iter,resvec] = cgs(A,b,...) также возвращает вектор остаточных норм в каждой итерации, включая norm(b-A*x0).
A = gallery('wilk',21);
b = sum(A,2);
tol = 1e-12; maxit = 15;
M1 = diag([10:-1:1 1 1:10]);
x = cgs(A,b,tol,maxit,M1);отображает сообщение
cgs converged at iteration 13 to a solution with relative residual 2.4e-016.
Этот пример заменяет матричный A в предыдущем примере с указателем на функцию матричного векторного произведения afun и предварительный формирователь M1 с указателем на функцию backsolve mfun. Пример содержится в файле run_cgs это
Вызовы cgs с указателем на функцию @afun в качестве его первого аргумента.
Содержит afun как вложенную функцию, так, чтобы все переменные в run_cgs были доступны afun и myfun.
Следующее показывает код для run_cgs:
function x1 = run_cgs
n = 21;
b = afun(ones(n,1));
tol = 1e-12; maxit = 15;
x1 = cgs(@afun,b,tol,maxit,@mfun);
function y = afun(x)
y = [0; x(1:n-1)] + ...
[((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...
[x(2:n); 0];
end
function y = mfun(r)
y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
end
endКогда вы входите
x1 = run_cgs
MATLAB возвращается
cgs converged at iteration 13 to a solution with relative residual 2.4e-016.
Этот пример демонстрирует использование предварительного формирователя.
Загрузите west0479, действительное 479 479 несимметричная разреженная матрица.
load west0479;
A = west0479;Задайте b так, чтобы истинное решение было вектором из всех единиц.
b = full(sum(A,2));
Установите допуск и максимальное количество итераций.
tol = 1e-12; maxit = 20;
Используйте cgs, чтобы найти решение в требуемом допуске и количестве итераций.
[x0,fl0,rr0,it0,rv0] = cgs(A,b,tol,maxit);
fl0 равняется 1, потому что cgs не сходится к требуемому допуску 1e-12 в требуемых 20 итерациях. На самом деле поведение cgs так плохо, что исходное предположение (x0 = zeros(size(A,2),1) является лучшим решением и возвращен, как обозначено it0 = 0. MATLAB хранит остаточную историю в rv0.
Постройте график поведения cgs.
semilogy(0:maxit,rv0/norm(b),'-o'); xlabel('Iteration number'); ylabel('Relative residual');
График показывает, что решение не сходится. Можно использовать предварительный формирователь, чтобы улучшить результат.
Создайте предварительный формирователь с ilu, поскольку A несимметричен.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));
Error using ilu There is a pivot equal to zero. Consider decreasing the drop tolerance or consider using the 'udiag' option.
MATLAB не может создать неполный LU, когда это привело бы к сингулярному фактору, который бесполезен как предварительный формирователь.
Можно попробовать еще раз с уменьшенным допуском отбрасывания, как обозначено сообщением об ошибке.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6)); [x1,fl1,rr1,it1,rv1] = cgs(A,b,tol,maxit,L,U);
fl1 0, потому что cgs управляет относительной невязкой к 4.3851e-014 (значение rr1). Относительная невязка является меньше, чем предписанный допуск 1e-12 в третьей итерации (значение it1), когда предобусловленный неполной LU-факторизацией с допуском отбрасывания 1e-6. Выводом rv1(1) является norm(b), и выводом rv1(14) является norm(b-A*x2).
Можно следовать, прогресс cgs путем графического изображения относительных невязок в каждой итерации, начинающей с первоначальной сметы (выполните итерации номера 0).
semilogy(0:it1,rv1/norm(b),'-o'); xlabel('Iteration number'); ylabel('Relative residual');
[1] Барретт, Ягода R., M., Т. F. Чан, и др., Шаблоны для Решения Линейных систем: Стандартные блоки для Итеративных Методов, SIAM, Филадельфия, 1994.
[2] Sonneveld, Питер, “CGS: быстрый решатель Lanczos-типа для несимметричных линейных систем”, SIAM J. Научный Закон Comput., январь 1989, Издание 10, № 1, стр 36–52.