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

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

rng default
A = sprand(600,600,.5);
A = A'*A + speye(size(A));
b = rand(600,1);

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

x = cgs(A,b);

cgs 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.068.

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

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

x = cgs(A,b,1e-4,40);

cgs stopped at iteration 40 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 39) has relative residual 0.0035.

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

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

L = ichol(A);
x = cgs(A,b,1e-4,40,L);

cgs converged at iteration 13 to a solution with relative residual 7.1e-05.

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

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

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

fl0 = 1

rr0

rr0 = 1

it0

it0 = 0

fl0 является 1 потому что cgs не сходится к запрошенному допуску 1e-12 в пределах запрошенных 20 итераций. На самом деле, поведение cgs настолько плох, что начальное предположение 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 в качестве входов в cgs.

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

fl1 = 0

rr1

rr1 = 4.3851e-14

it1

it1 = 3

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

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

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

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

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

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

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

cgs converged at iteration 17 to a solution with relative residual 8.8e-07.

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

cgs converged at iteration 4 to a solution with relative residual 8e-07.

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

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

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

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

cgs converged at iteration 11 to a solution with relative residual 1.3e-14.

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