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

Создайте случайную разреженную матрицу A с 50%-й плотностью. Также создайте случайный векторный b для правой стороны $\mathrm{Ax}=\mathit{b}$.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

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

x = gmres(A,b);

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

gmres по умолчанию использование 10 итераций и допуск 1e-6, и алгоритм не может сходиться в тех 10 итерациях для этой матрицы. Поскольку невязка является все еще большой, это - хороший индикатор, что необходимо больше итераций (или матрица перед формирователем). Также можно уменьшать допуск, чтобы облегчить для алгоритма сходиться.

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

tol = 1e-4;
maxit = 100;
x = gmres(A,b,[],tol,maxit);

gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.052.

Даже с более свободным допуском и большим количеством итераций остаточная ошибка не улучшается очень. Когда итеративный алгоритм останавливается этим способом, это - хорошая индикация, что матрица перед формирователем необходима. Однако gmres также имеет вход, который управляет количеством внутренних итераций. Путем определения значения для внутренних итераций, gmres действительно больше работает на внешнюю итерацию.

Решите систему снова с помощью restart значение 100 и maxit значение 250. Вместо того, чтобы делать 100 итераций однажды, gmres выполняет 100 итераций между перезапусками и повторяет это 250 раз.

restart = 100;
maxit = 250;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 133 (inner iteration 55) to a solution with relative residual 0.0001.

В этом случае, задающем большое значение перезапуска для gmres позволяет ему сходиться к решению в позволенном количестве итераций. Однако большие значения перезапуска могут использовать большую память когда A является также большим.

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

Загрузите west0479, действительное 479 479 несимметричная разреженная матрица.

load west0479
A = west0479;

Задайте b так, чтобы истинное решение было вектором из всех единиц.

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

fl0 = 1

rr0

rr0 = 0.7603

it0

it0 = 1×2

     1    20

fl0 1 потому что gmres не сходится к требуемому допуску 1e-12 в требуемых 20 итерациях. Лучшее приближенное решение, что gmres возвраты являются последним (как обозначено it0(2) = 20). MATLAB хранит остаточную историю в rv0.

Чтобы помочь с медленной сходимостью, можно задать матрицу перед формирователем. Начиная с A несимметрично, используйте ilu сгенерировать предварительный формирователь $\mathit{M}=\mathit{L}\mathit{U}$. Задайте допуск отбрасывания, чтобы проигнорировать недиагональные записи со значениями, меньшими, чем 1e-6. Решите предобусловленную систему ${\mathit{M}}^{-1}\mathit{A}\mathit{x}={\mathit{M}}^{-1}\mathit{b}$ путем определения L и U как вводит к gmres. Обратите внимание на то, что, когда вы задаете предварительный формирователь, gmres вычисляет норму невязки предобусловленной системы для выходных параметров rr1 и rv1.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1

fl1 = 0

rr1

rr1 = 7.2303e-14

it1

it1 = 1×2

     1     6

Использование ilu предварительный формирователь производит относительную невязку меньше, чем предписанный допуск 1e-12 в шестой внешней итерации. Первый остаточный rv1(1) norm(U\(L\b)), где M = L*U. Последний остаточный rv1(end) norm(U\(L\(b-A*x1))).

Можно следовать за прогрессом gmres путем графического вывода относительных остаточных значений в каждой итерации. Постройте остаточную историю каждого решения с линией для заданного допуска.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Используя предварительный формирователь с перезапущенным gmres.

Загрузите west0479, действительное 479 479 несимметричная разреженная матрица.

load west0479
A = west0479;

Задайте b так, чтобы истинное решение было вектором из всех единиц.

b = sum(A,2);

Создайте неполный предварительный формирователь LU с допуском отбрасывания 1e-6.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

Преимущество к использованию перезапущенного gmres должен ограничить объем памяти, требуемый выполнить метод. Без перезапуска, gmres требует maxit векторы устройства хранения данных, чтобы сохранить основание подпространства Крылова. Кроме того, gmres должен ортогонализировать против всех предыдущих векторов на каждом шаге. Перезапуск ограничивает количество используемой рабочей области и объем работы, сделанный на внешнюю итерацию.

Выполните gmres(3), gmres(4), и gmres(5) использование неполных факторов LU как предварительные формирователи. Используйте допуск 1e-12 и максимум 20 внешних итераций.

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3

fl3 = 0

fl4

fl4 = 0

fl5

fl5 = 0

fl3, fl4, и fl5 весь 0, потому что в каждом случае перезапустил gmres управляет относительной невязкой к меньше, чем предписанный допуск 1e-12.

Следующий график показывает, что история сходимости каждого перезапустила gmres метод. gmres(3) сходится во внешней итерации 5, внутренняя итерация 3 (it3 = [5, 3]) который совпал бы с внешней итерацией 6, внутренняя итерация 0, следовательно маркировка 6 на итоговой отметке деления.

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

В целом, чем больше количество внутренних итераций, тем больше работает gmres делает на внешнюю итерацию, и быстрее она может сходиться.

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

Создайте трехдиагональную разреженную матрицу. Используйте сумму каждой строки как вектор для правой стороны $\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);

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

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

gmres converged at iteration 27 to a solution with relative residual 9.5e-07.

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

gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

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

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

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

Один из Уилкинсона тестирует матрицы, сгенерированные gallery 21 21 трехдиагональная матрица. Предварительно просмотрите матрицу.

A = gallery('wilk',21)

A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

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

$\left[\begin{array}{cccccccccc}10& 1& 0& \cdots & & \cdots & & \cdots & 0& 0\\ 1& 9& 1& 0& & & & & & 0\\ 0& 1& 8& 1& 0& & & & & ⋮\\ ⋮& 0& 1& 7& 1& 0& & & & \\ & & 0& 1& 6& 1& 0& & & ⋮\\ ⋮& & & 0& 1& 5& 1& 0& & \\ & & & & 0& 1& 4& 1& 0& ⋮\\ ⋮& & & & & 0& 1& 3& \ddots & 0\\ 0& & & & & & 0& \ddots & \ddots & 1\\ 0& 0& \cdots & & \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ {\mathit{x}}_4\\ {\mathit{x}}_5\\ ⋮\\ \\ ⋮\\ \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

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

$\left[\begin{array}{c}0+10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}+0\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}10{\mathit{x}}_1\\ 9{\mathit{x}}_2\\ ⋮\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

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

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

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

Теперь решите линейную систему $\mathrm{Ax}=\mathit{b}$ путем обеспечения gmres с указателем на функцию, который вычисляет A*x. Используйте допуск 1e-12, 15 внешних итераций и 10 внутренних итераций перед перезапуском.

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)

gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.

x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

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

afun(x1)

ans = 21×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:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end