Документация

x = pcg(A,b) попытки решить систему линейных уравнений A*x = b для x использование Предобусловленного Метода сопряженных градиентов. Когда попытка успешна, pcg отображает сообщение, чтобы подтвердить сходимость. Если pcg сбои, чтобы сходиться после максимального количества итераций или остановов по любой причине, это отображает диагностическое сообщение, которое включает относительный остаточный norm(b-A*x)/norm(b) и номер итерации, в который остановленный метод.

x = pcg(A,b,tol) задает допуск к методу. Допуском по умолчанию является 1e-6.

x = pcg(A,b,tol,maxit) задает максимальное количество итераций, чтобы использовать. pcg отображает диагностическое сообщение, если ему не удается сходиться в maxit итерации.

x = pcg(A,b,tol,maxit,M) задает матрицу перед формирователем M и вычисляет x путем эффективного решения системы $H^{- 1} A H^{- T} y = H^{- 1} b$ для y, где $y = H^{T} x$ и $H = M^{1 / 2} = {(M_{1} M_{2})}^{1 / 2}$ . Используя предварительный формирователь матрица может улучшить числовые свойства проблемы и КПД вычисления.

x = pcg(A,b,tol,maxit,M1,M2) задает факторы матрицы перед формирователем M таким образом, что M = M1*M2.

x = pcg(A,b,tol,maxit,M1,M2,x0) задает исходное предположение для вектора решения x. Значением по умолчанию является нулевой вектор.

[x,flag] = pcg(___) возвращает флаг, который задает, сходился ли алгоритм успешно. Когда flag = 0, сходимость была успешна. Можно использовать этот выходной синтаксис с любой из предыдущих комбинаций входных аргументов. Когда вы задаете flag выведите, pcg не отображает диагностических сообщений.

[x,flag,relres] = pcg(___) также возвращает относительный остаточный norm(b-A*x)/norm(b). Если flag 0, затем relres <= tol.

[x,flag,relres,iter] = pcg(___) также возвращает номер итерации iter в котором x был вычислен.

[x,flag,relres,iter,resvec] = pcg(___) также возвращает вектор нормы невязки в каждой итерации, включая первый остаточный norm(b-A*x0).

Примеры

Итеративное решение линейной системы

Решите квадратную линейную систему с помощью pcg с настройками по умолчанию, и затем настраивают допуск и количество итераций, используемых в процессе решения.

Создайте случайную разреженную матрицу A с 50%-й плотностью. Также создайте векторный b из сумм строки A для правой стороны $Ax = b$ так, чтобы истинное решение $x$ вектор из единиц.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = sum(A,2);

Решить $Ax = b$ использование pcg. Выходное отображение включает значение относительной остаточной ошибки $\frac{‖ Ax - b ‖}{‖ b ‖}$ .

x = pcg(A,b);

pcg stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 3.6e-06.

pcg по умолчанию использование 20 итераций и допуск 1e-6, и алгоритм не может сходиться в тех 20 итерациях для этой матрицы. Однако невязка близко к допуску, таким образом, для алгоритма, вероятно, только нужно больше итераций, чтобы сходиться.

Решите систему снова с помощью допуска 1e-7 и 150 итераций.

x = pcg(A,b,1e-7,150);

pcg converged at iteration 130 to a solution with relative residual 9.9e-08.

Используя `pcg` с Предварительным формирователем

Исследуйте эффект использования матрицы перед формирователем с pcg решить линейную систему.

Создайте симметричную положительную определенную, полосную матрицу коэффициентов.

A = delsq(numgrid('S',102));

Задайте b для правой стороны линейного уравнения.

b = ones(size(A,1),1);

Установите погрешность и максимальное количество итераций.

tol = 1e-8;
maxit = 100;

Используйте pcg найти решение в требуемом допуске и количестве итераций. Задайте пять выходных параметров, чтобы возвратить информацию о процессе решения:

x0 вычисленное решение A*x0 = b.
fl0 флаг, указывающий, сходился ли алгоритм.
rr0 невязка вычисленного ответа x0.
it0 номер итерации когда x0 был вычислен.
rv0 вектор остаточной истории для $‖ Ax - b ‖$ .

[x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);
fl0

fl0 = 1

rr0

rr0 = 0.0131

it0

it0 = 100

fl0 1 потому что pcg не сходится к требуемому допуску 1e-8 в требуемых 100 итерациях.

Чтобы помочь с медленной сходимостью, можно задать матрицу перед формирователем. Начиная с A симметрично, используйте ichol сгенерировать предварительный формирователь $M = L L^{T}$ . Решите предобусловленную систему путем определения L и L' как вводит к pcg.

L = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L,L');
fl1

fl1 = 0

rr1

rr1 = 8.0992e-09

it1

it1 = 79

Использование ichol предварительный формирователь производит относительную невязку меньше, чем предписанный допуск 1e-8 в 79-й итерации. Выход rv1(1) norm(b) и rv1(end) norm(b-A*x1).

Теперь используйте michol опция, чтобы создать модифицированный неполный предварительный формирователь Холесского.

L = ichol(A,struct('michol','on'));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L,L');
fl2

fl2 = 0

rr2

rr2 = 9.9605e-09

it2

it2 = 47

Этот предварительный формирователь лучше, чем тот, произведенный неполной факторизацией Холесского с нулевой заливкой для матрицы коэффициентов в этом примере, таким образом, pcg может сходиться еще более быстрый.

Вы видите, как предварительные формирователи влияют на уровень сходимости pcg путем графического вывода каждой из остаточных историй, начинающих с первоначальной оценки (выполняют итерации номера 0). Добавьте линию для заданного допуска.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),'-o')
semilogy(0:length(rv2)-1,rv2/norm(b),'-o')
yline(tol,'r--');
legend('No Preconditioner','Default ICHOL','Modified ICHOL','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Предоставление исходного предположения

Исследуйте эффект предоставления pcg с исходным предположением решения.

Создайте трехдиагональную разреженную матрицу. Используйте сумму каждой строки как вектор для правой стороны $Ax = b$ так, чтобы ожидаемое решение для $x$ вектор из единиц.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

Используйте pcg решить $Ax = b$ дважды: одно время с исходным предположением по умолчанию, и одно время с хорошим исходным предположением решения. Используйте 200 итераций в обоих решениях и задайте исходное предположение как вектор со всеми элементами, равными 0.99.

maxit = 200;
x1 = pcg(A,b,[],maxit);

pcg converged at iteration 35 to a solution with relative residual 9.5e-07.

x0 = 0.99*ones(size(A,2),1);
x2 = pcg(A,b,[],maxit,[],[],x0);

pcg converged at iteration 7 to a solution with relative residual 8.7e-07.

В этом случае, предоставляющем исходное предположение, включает pcg сходиться более быстро.

Возвращение промежуточных результатов

Также можно использовать исходное предположение, чтобы получить промежуточные результаты путем вызова pcg в цикле for. Каждый вызов решателя выполняет несколько итераций и хранит расчетное решение. Затем вы используете то решение в качестве начального вектора для следующего пакета итераций.

Например, этот код выполняет 100 итераций четыре раза и хранит вектор решения после каждой передачи в цикле for:

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = pcg(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

X(:,k) вектор решения, вычисленный в итерации k из цикла for и R(k) относительная невязка того решения.

Используя указатель на функцию вместо числовой матрицы

Решите линейную систему путем обеспечения pcg с указателем на функцию, который вычисляет A*x вместо матрицы коэффициентов A.

Используйте gallery сгенерировать 20 20 положительную определенную трехдиагональную матрицу. Супер - и поддиагонали имеют единицы, в то время как основные диагональные элементы рассчитывают по сравнению с от 20 до 1. Предварительно просмотрите матрицу.

n = 20;
A = gallery('tridiag',ones(n-1,1),n:-1:1,ones(n-1,1));
full(A)

ans = 20×20

    20     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1    19     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1    18     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1    17     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1    16     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1    15     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1    14     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1    13     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1    12     1     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1    11     1     0     0     0     0     0     0     0     0     0
      ⋮

Поскольку эта трехдиагональная матрица имеет специальную структуру, можно представлять операцию A*x с указателем на функцию. Как каждая строка A умножает элементы в x, только несколько результатов являются ненулевыми (соответствие ненулям на tridiagonals). Кроме того, только основная диагональ имеет ненули, которые не равны 1. Выражение $Ax$ становится:

$[\begin{array}{c} 20 & 1 & 0 & \dots & \dots & \dots & 0 & 0 \\ 1 & 19 & 1 & 0 & 0 \\ 0 & 1 & 18 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 17 & 1 & 0 \\ 0 & 1 & 16 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 15 & 1 & 0 \\ 0 & 1 & 14 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 13 & ⋱ & 0 \\ 0 & 0 & ⋱ & ⋱ & 1 \\ 0 & 0 & \dots & \dots & \dots & 0 & 1 & 1 \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ ⋮ \\ ⋮ \\ x_{20} \end{array}] = [\begin{array}{c} 2 {0 x}_{1} + x_{2} \\ x_{1} + 19 x_{2} + x_{3} \\ x_{2} + 18 x_{3} + x_{4} \\ ⋮ \\ x_{18} + 2 x_{19} + x_{20} \\ x_{19} + x_{20} \end{array}]$ .

Итоговый вектор может быть записан как сумма трех векторов:

$[\begin{array}{c} 2 {0 x}_{1} + x_{2} \\ x_{1} + 19 x_{2} + x_{3} \\ x_{2} + 18 x_{3} + x_{4} \\ ⋮ \\ x_{18} + 2 x_{19} + x_{20} \\ x_{19} + x_{20} \end{array}]$ = $[\begin{array}{c} 0 \\ x_{1} \\ ⋮ \\ x_{19} \end{array}] + [\begin{array}{c} 20 x_{1} \\ 19 x_{2} \\ ⋮ \\ x_{20} \end{array}] + [\begin{array}{c} x_{2} \\ ⋮ \\ x_{20} \\ 0 \end{array}]$ .

В MATLAB® запишите функцию, которая создает эти векторы и добавляет их вместе, таким образом давая значение A*x:

function y = afun(x)
y = [0; x(1:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 0];
end

(Эта функция сохранена как локальная функция в конце примера.)

Теперь решите линейную систему $Ax = b$ путем обеспечения pcg с указателем на функцию, который вычисляет A*x. Используйте допуск 1e-12 и 50 итераций.

b = ones(20,1);
tol = 1e-12;  
maxit = 50;
x1 = pcg(@afun,b,tol,maxit)

pcg converged at iteration 20 to a solution with relative residual 4.4e-16.

Проверяйте тот afun(x1) дает вектор из единиц.

afun(x1)

Локальные функции

function y = afun(x)
y = [0; x(1:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 0];
end

Входные параметры

`A` — Матрица коэффициентов
матрица | указатель на функцию

Матрица коэффициентов в виде симметричной положительной определенной матрицы или указателя на функцию. Эта матрица является матрицей коэффициентов в линейной системе A*x = b. Обычно A большая разреженная матрица или указатель на функцию, который возвращает продукт большой разреженной матрицы и вектор-столбца. Смотрите Определяют, Является ли Матрица Симметричной Положительный Определенный для получения информации о том, как подтвердить тот A симметричен положительный определенный.

Определение `A` как указатель на функцию

Можно задать матрицу коэффициентов как указатель на функцию вместо матрицы, чтобы сохранить память в вычислении. Указатель на функцию возвращает матричные векторные произведения вместо того, чтобы формировать целую матрицу коэффициентов, делая вычисление более эффективным.

Чтобы использовать указатель на функцию, используйте функциональную подпись function y = afun(x). Параметризация Функций объясняет, как предоставить дополнительные параметры функциональному afun, при необходимости. Вызов функции afun(x) должен возвратить значение A*x.

Типы данных: double | function_handle
Поддержка комплексного числа: Да

`b` — Правая сторона линейного уравнения
вектор

Правая сторона линейного уравнения в виде вектор-столбца. b должен быть вектор-столбец с длиной, равной size(A,1).

Типы данных: double
Поддержка комплексного числа: Да

`tol` — Допуск метода
`[]` или `1e-6` (значение по умолчанию) | положительная скалярная величина

Допуск метода в виде положительной скалярной величины. Используйте этот вход для точности компромисса и времени выполнения в вычислении. pcg должен соответствовать допуску в количестве позволенных итераций, чтобы быть успешным. Меньшее значение tol означает, что ответ должен быть более точным для вычисления, чтобы быть успешным.

Типы данных: double

`maxit` — Максимальное количество итераций
`[]` или `min(size(A,1),20)` (значение по умолчанию) | положительное скалярное целое число

Максимальное количество итераций в виде положительного скалярного целого числа. Увеличьте значение maxit позволить больше итераций для pcg соответствовать допуску tol. Обычно меньшее значение tol средние значения больше итераций требуются, чтобы успешно завершать вычисление.

`M`, `M1`, `M2` — Матрицы перед формирователем (в качестве отдельных аргументов)
`eye(size(A))` (значение по умолчанию) | матрицы | указатели на функцию

Матрицы перед формирователем в виде отдельных аргументов матриц или указателей на функцию. Можно задать матрицу перед формирователем M или его матричные факторы M = M1*M2 улучшить числовые аспекты линейной системы и облегчить для pcg сходиться быстро. Можно использовать неполные матричные функции факторизации ilu и ichol сгенерировать матрицы перед формирователем. Также можно использовать equilibrate до факторизации, чтобы улучшить число обусловленности матрицы коэффициентов. Для получения дополнительной информации о предварительных формирователях смотрите Итерационные методы для Линейных систем.

pcg обрабатывает незаданные предварительные формирователи как единичные матрицы.

Определение `M` как указатель на функцию

Можно задать любой M, M1, или M2 как указатели на функцию вместо матриц, чтобы сохранить память в вычислении. Указатель на функцию выполняет матрично-векторные операции вместо того, чтобы формировать целую матрицу перед формирователем, делая вычисление более эффективным.

Чтобы использовать указатель на функцию, используйте функциональную подпись function y = mfun(x). Параметризация Функций объясняет, как предоставить дополнительные параметры функциональному mfun, при необходимости. Вызов функции mfun(x) должен возвратить значение M\x или M2\(M1\x).

Типы данных: double | function_handle
Поддержка комплексного числа: Да

`x0` — Исходное предположение
`[]` или вектор-столбец из нулей (значение по умолчанию) | вектор-столбец

Исходное предположение в виде вектор-столбца с длиной равняется size(A,2). Если можно обеспечить pcg с более разумным исходным предположением x0 чем нулевой вектор по умолчанию затем это может сохранить время вычисления и помочь алгоритму сходиться быстрее.

Типы данных: double
Поддержка комплексного числа: Да

Выходные аргументы

`x` — Решение для линейной системы
вектор

Решение для линейной системы, возвращенное как вектор. Этот выход дает приближенное решение линейной системы A*x = b. Если вычисление успешно (flag = 0), затем relres меньше чем или равно tol.

Каждый раз, когда вычисление не успешно (flag ~= 0), решение x возвращенный pcg тот с минимальной нормой невязки, вычисленной по всем итерациям.

`flag` — Флаг Convergence
скаляр

Флаг Convergence, возвращенный как одно из скалярных значений в этой таблице. Флаг сходимости указывает, было ли вычисление успешно и дифференцируется между несколькими различными формами отказа.

Флаговое значение	Сходимость
0	Успех — `pcg` сходившийся к желаемому допуску `tol` в `maxit` итерации.
1	Отказ — `pcg` выполненный с помощью итераций `maxit` итерации, но не сходились.
2	Отказ — матрица перед формирователем `M` или `M = M1*M2` isIllConditioned.
3	Отказ — `pcg` застоявшийся после того, как две последовательных итерации были тем же самым.
4	Отказ — Один из скаляров, вычисленных `pcg` алгоритм стал слишком маленьким или слишком большим, чтобы продолжить вычислять.

`relres` — Относительная остаточная ошибка
скаляр

Относительная остаточная ошибка, возвращенная как скаляр. Относительная остаточная ошибка relres = norm(b-A*x)/norm(b) индикация относительно того, насколько точный ответ. Если вычисление сходится к допуску tol в maxit итерации, затем relres <= tol.

Типы данных: double

`iter` — Номер итерации
скаляр

Номер итерации, возвращенный как скаляр. Этот выход указывает на номер итерации в который вычисленный ответ для x был вычислен.

Типы данных: double

`resvec` — Остаточная ошибка
вектор

Остаточная ошибка, возвращенная как вектор. Остаточная ошибка norm(b-A*x) показывает, как близко алгоритм к схождению для данного значения x. Число элементов в resvec равно количеству итераций. Можно исследовать содержимое resvec помочь решить, изменить ли значения tol или maxit.

Типы данных: double

Больше о