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

Создание разреженной тридиагональной матрицы A в качестве матрицы коэффициентов. Используйте суммы плотных строк $\mathit{A}$ в качестве вектора для правой стороны $\mathrm{Ax}\mathit{=}$ b, так что $\mathit{решение}$ x должно быть вектором единиц.

n = 400; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = full(sum(A,2));

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

x = symmlq(A,b);

symmlq 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.045.

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

Повторное решение системы с использованием допуска 1e-4 и 250 итераций.

x = symmlq(A,b,1e-4,250);

symmlq converged at iteration 199 to a solution with relative residual 1.5e-14.

Анализ влияния использования матрицы предварительного кондиционирования с помощью symmlq для решения линейной системы.

Создайте симметричную матрицу положительных определенных, ограниченных коэффициентов.

A = delsq(numgrid('S',102));

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 100;

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

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

fl0 = 1

rr0

rr0 = 0.0031

it0

it0 = 100

fl0 равно 1, потому что symmlq не сходится с требуемым допуском 1e-12 в пределах запрошенных 100 итераций.

Для облегчения сходимости можно указать матрицу предварительного условия. С тех пор A симметричен, используется ichol для создания устройства предварительного кондиционирования $\mathit{М}\mathit{=}{\mathit{L}}^{\mathit{}}$LT. Укажите 'ict' использование неполной факторизации Холеского с понижением порога и указание значения диагонального сдвига 1e-6 во избежание непозволительных поворотов. Решение предварительно кондиционированной системы путем указания L и L' в качестве входных данных для symmlq.

setup = struct('type','ict','diagcomp',1e-6,'droptol',1e-14);
L = ichol(A,setup);
[x1,fl1,rr1,it1,rv1,rvcg1] = symmlq(A,b,tol,maxit,L,L');
fl1

fl1 = 0

rr1

rr1 = 2.7782e-15

it1

it1 = 3

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

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

semilogy(0:length(rvcg0)-1,rvcg0/norm(b),'-o')
hold on
semilogy(0:length(rvcg1)-1,rvcg1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ICHOL preconditioner','Tolerance','Location','SouthEast')
xlabel('Iteration number')
ylabel('Relative residual')

Проверить влияние подачи symmlq с первоначальным предположением решения.

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

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

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

symmlq converged at iteration 34 to a solution with relative residual 9.5e-07.

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

symmlq converged at iteration 6 to a solution with relative residual 8.7e-07.

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

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

Решение линейной системы путем предоставления symmlq с дескриптором функции, который вычисляет 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{}\begin{array}{cccccccccc}\mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & \\ & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & \\ & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \\ & & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ & & & & & \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}\\ \mathrm{}& & & & & & \mathrm{}& & & \mathrm{}\\ \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}& \mathrm{Ax=[1010\cdots \cdots \cdots 001910001810⋮⋮0171001610⋮⋮0151001410⋮⋮013\ddots 000\ddots \ddots 100\cdots \cdots \cdots 0110}\end{array}]\begin{array}{c}{\mathit{[}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ \\ \\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}x1x2x3x4x5⋮⋮x21]=[10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21$].

Результирующий вектор может быть записан как сумма трех векторов:

$\mathrm{}\begin{array}{c}\mathrm{}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}\end{array}$Ax=[0+10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21+0]=[0x1⋮x20]+[10x19x2⋮10x21]+[x2⋮x210$\begin{array}{c}\mathrm{}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}\mathrm{}{\mathit{}}_{\mathrm{}}\\ \mathrm{}{\mathit{}}_{\mathrm{}}\\ \\ \mathrm{}{\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\\ \mathrm{}\end{array}$].

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

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

symmlq converged at iteration 10 to a solution with relative residual 4.5e-16.

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