Уравновесьте матрицу с большим числом обусловленности, чтобы улучшить эффективность и стабильность решения линейной системы с помощью итерационного решателя gmres.

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

load west0479
A = west0479;
c1 = condest(A)

c1 = 1.4244e+12

Попытайтесь решить линейную систему $\mathrm{Ax}=\mathit{b}$ использование gmres с 450 итерациями и допуском 1e-11. Задайте пять выходов так, чтобы gmres возвращает нормы невязки решения при каждой итерации (использование ~ для подавления ненужных выходов). Постройте график норм невязки на полулогарифмический график. График показывает, что gmres не в состоянии достичь разумной нормы невязки, и поэтому расчетное решение для $\mathit{x}$ не надежно.

b = ones(size(A,1),1);
tol = 1e-11;
maxit = 450;
[x,flx,~,~,rvx] = gmres(A,b,[],tol,maxit);
semilogy(rvx)
title('Residual Norm at Each Iteration')

Использование equilibrate для транспозиции и пересмотра A. Создайте новую матрицу B = R*P*A*C, который имеет лучшее число обусловленности и диагональные элементы только 1 и -1.

[P,R,C] = equilibrate(A);
B = R*P*A*C;
c2 = condest(B)

c2 = 5.1036e+04

Использование выходов, возвращенной equilibrate, можно переформулировать задачу $\mathrm{Ax}=\mathit{b}$ в $\mathrm{By}=\mathit{d}$, где $\mathit{B}=\mathrm{RPAC}$ и $\mathit{d}=\mathrm{RPb}$. В этой форме вы можете восстановить решение в исходной системе с $\mathit{x}=\mathrm{Cy}$.

Использование gmres решить $\mathrm{By}=\mathit{d}$ для $\mathit{y}$, и затем повторите нормы невязки при каждой итерации. График показывает, что уравновешивание матрицы улучшает стабильность задачи с gmres сходится к желаемому допуску 1e-11 менее чем за 200 итераций.

d = R*P*b;
[y,fly,~,~,rvy] = gmres(B,d,[],tol,maxit);
hold on
semilogy(rvy)
legend('Original', 'Equilibrated', 'Location', 'southeast')
title('Relative Residual Norms (No Preconditioner)')
hold off

semilogy(rvy)
hold on

[L1,U1] = ilu(B,struct('type','ilutp','droptol',1e-1,'thresh',0));
[yp1,flyp1,~,~,rvyp1] = gmres(B,d,[],tol,maxit,L1,U1);
semilogy(rvyp1)

[L2,U2] = ilu(B,struct('type','ilutp','droptol',1e-2,'thresh',0));
[yp2,flyp2,~,~,rvyp2] = gmres(B,d,[],tol,maxit,L2,U2);
semilogy(rvyp2)

legend('No preconditioner', 'ILUTP(1e-1)', 'ILUTP(1e-2)')
title('Relative Residual Norms with ILU Preconditioner (Equilibrated)')
hold off