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

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

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

x = minres(A,b);

minres 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.017.

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

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

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

minres converged at iteration 200 to a solution with relative residual 7e-13.

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

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

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

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 100;

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

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

fl0 = 1

rr0

rr0 = 0.0013

it0

it0 = 100

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

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

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

fl1 = 0

rr1

rr1 = 2.6691e-15

it1

it1 = 4

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

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

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','ICHOL preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

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

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

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

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

minres converged at iteration 27 to a solution with relative residual 9.5e-07.

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

minres converged at iteration 7 to a solution with relative residual 6.7e-07.

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

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

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

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

minres converged at iteration 11 to a solution with relative residual 4.1e-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