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

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{M}=\mathit{L}{\mathit{L}}^{\mathit{T}}$. Задайте '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{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}$ путем обеспечения 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