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

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

x = tfqmr(A,b);

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

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

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

x = tfqmr(A,b,1e-4,100);

tfqmr stopped at iteration 200 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 13) has relative residual 0.3.

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

Вычислите неполную факторизацию Холесского A, и использовать L' коэффициент как вход предкондиционера для tfqmr.

L = ichol(A);
x = tfqmr(A,b,1e-4,100,L');

tfqmr converged at iteration 32 to a solution with relative residual 5.2e-05.

Использование предкондиционера улучшает числовые свойства задачи достаточно, чтобы tfqmr способен сходиться.

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

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

load west0479
A = west0479;

Определите b чтобы истинное решение $\mathrm{Ax}=\mathit{b}$ является вектором всех таковых.

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

[x,fl0,rr0,it0,rv0] = tfqmr(A,b,tol,maxit);
fl0

fl0 = 1

rr0

rr0 = 0.9846

it0

it0 = 4

fl0 равен 1, потому что tfqmr не сходится к запрошенному допуску 1e-12 в пределах запрошенных 20 итераций. Десятая итерация является лучшим аппроксимационным решением и возвращается как обозначено it0 = 10.

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

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

fl1 = 0

rr1

rr1 = 4.3298e-14

it1

it1 = 3

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

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

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

Исследуйте эффект подачи tfqmr с начальным угадыванием решения.

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

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

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

tfqmr converged at iteration 19 to a solution with relative residual 9.6e-07.

x0 = 0.99*e;
x2 = tfqmr(A,b,[],maxit,[],[],x0);

tfqmr converged at iteration 4 to a solution with relative residual 7.9e-07.

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

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

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

Выражение $\mathrm{Ax}$ становится:

$\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]$.

Получившийся вектор может быть записан как сумма трех векторов:

$\mathrm{Ax}=\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}$ путем предоставления tfqmr с указателем на функцию, который вычисляет A*x. Используйте допуск 1e-12 и 50 итераций.

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

tfqmr converged at iteration 10 to a solution with relative residual 6.7e-15.

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