Решите квадратную линейную систему с помощью 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 \mathrm{Ax}-\mathit{b}\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 так, чтобы истинное решение было вектором из всех единиц.

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

[x0,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*ones(size(A,2),1);
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 умножает элементы в 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}$ путем обеспечения 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