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

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

x = bicgstab(A,b);

bicgstab 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 0.13.

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

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

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

bicgstab 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.044.

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

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

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

bicgstab converged at iteration 30.5 to a solution with relative residual 7.7e-05.

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

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

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

fl0 = 1

rr0

rr0 = 1

it0

it0 = 0

fl0 является 1 потому что bicgstab не сходится к запрошенному допуску 1e-12 в пределах запрошенных 20 итераций. На самом деле, поведение bicgstab настолько плох, что начальное предположение x0 = zeros(size(A,2),1) является лучшим решением и возвращается, как указано it0 = 0.

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

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

fl1 = 0

rr1

rr1 = 5.5615e-14

it1

it1 = 3

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

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

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')

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

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

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

maxit = 50;
x1 = bicgstab(A,b,[],maxit);

bicgstab converged at iteration 20.5 to a solution with relative residual 9.3e-07.

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

bicgstab converged at iteration 4 to a solution with relative residual 8.7e-07.

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

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

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

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

bicgstab converged at iteration 11.5 to a solution with relative residual 5.2e-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