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

Создание случайной симметричной разреженной матрицы A. Также создайте вектор b строковых сумм A для правой стороны $\mathrm{Ax}\mathit{=}$ b, так что истинное $\mathit{решение}$ x является вектором единиц.

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

Решить $\mathrm{Ax}\mathit{=}$ b с помощьюpcg. Выходной дисплей включает в себя значение относительной остаточной ошибки $\frac{\Vert \mathit{}\mathrm{}b-Ax}{\mathit{\Vert }}$‖ 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 для решения линейной системы.

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

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

Определить b для правой стороны линейного уравнения $\mathrm{Ax}\mathit{=}$b.

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

Задайте допуск и максимальное количество итераций.

tol = 1e-8;
maxit = 100;

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

[x,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 для создания устройства предварительного кондиционирования $\mathit{М}\mathit{=}{\mathit{L}}^{\mathit{}}$LT. Решение предварительно кондиционированной системы путем указания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 для создания измененного неполного предварительного условия Cholesky.

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 с первоначальным предположением решения.

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

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

Использовать pcg для решения $\mathrm{Ax}\mathit{=}$ 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*e;
x2 = pcg(A,b,[],maxit,[],[],x0);

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

В этом случае предоставление начального предположения позволяет pcg для более быстрого сближения.

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

Решение линейной системы путем предоставления 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 умножает вектор, большинство элементов в полученном векторе равны нулю. Ненулевые элементы в результате соответствуют ненулевым тридиагональным элементам A. При этом только основная диагональ имеет ненулевые значения, не равные 1.

Выражение $\mathrm{Ax}$ становится следующим:

$\begin{array}{cccccccccc}\mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & \\ & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & \\ & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \\ & & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ & & & & & \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}\\ \mathrm{}& & & & & & \mathrm{}& & & \mathrm{}\\ \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ \\ \\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}\mathrm{}{\mathrm{}\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{[2010\cdots \cdots \cdots 00119100011810⋮⋮011710011610⋮⋮011510011410⋮⋮0113\ddots 000\ddots \ddots 100\cdots \cdots \cdots 011][x1x2x3x4x5⋮⋮x20]=[20x1+x2x1+19x2+x3x2+18x3+x4⋮x18+2x19+x20x19+x20}}\end{array}]$.

Результирующий вектор может быть записан как сумма трех векторов:

$\begin{array}{c}\mathrm{}{\mathrm{}\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\end{array}[20x1+x2x1+19x2+x3x2+18x3+x4⋮x18+2x19+x20x19+x20]$= $\begin{array}{c}\mathrm{}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}[0x1⋮x19]+\begin{array}{c}\mathrm{}{\mathit{}}_{\mathrm{}}\\ \mathrm{}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}[20x119x2⋮x20]+\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\\ \mathrm{}\end{array}[x2⋮x200]$.

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

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

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

Теперь решите линейную систему $\mathrm{Ax}\mathit{=}$ 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.

x1 = 20×1

    0.0476
    0.0475
    0.0500
    0.0526
    0.0555
    0.0588
    0.0625
    0.0666
    0.0714
    0.0769
      ⋮

Проверьте, что afun(x1) производит вектор из единиц.

afun(x1)

ans = 20×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

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