Zbieżność klasycznych solwerów iteracyjnych dla układów liniowych jest określona przez promień widmowy macierzy iteracji, . W przypadku ogólnego układu liniowego trudno jest określić optymalny (lub nawet dobry) parametr SOR ze względu na trudność w określeniu promienia widmowego macierzy iteracji. Poniżej zamieściłem wiele dodatkowych szczegółów, w tym przykład prawdziwego problemu, w którym znana jest optymalna waga SOR.ρ(G)
Promień widmowy i zbieżność
Promień spektralny jest zdefiniowany jako wartość bezwzględna wartości własnej największej wielkości. Metoda zbiegnie się, jeśli a mniejszy promień widmowy oznacza szybszą zbieżność. SOR działa poprzez zmianę podziału macierzy stosowanego do uzyskania macierzy iteracji w oparciu o wybór parametru ważenia , miejmy nadzieję, że zmniejsza promień widmowy powstałej macierzy iteracji.ρ<1ω
Podział macierzy
W poniższej dyskusji założę, że system do rozwiązania jest podany przez
Ax=b,
z iteracją formularza
x(k+1)=v+Gx(k),
gdzie jest wektorem, a liczba iteracji jest oznaczona .vkx(k)
SOR przyjmuje średnią ważoną starej iteracji i iteracji Gaussa-Seidela. Metoda Gaussa-Seidela polega na podziale macierzy na formę
A=D+L+U
gdzie to przekątna , to dolna trójkątna matryca zawierająca wszystkie elementy ściśle poniżej przekątnej, a to górna trójkątna matryca zawierający wszystkie elementy ściśle nad przekątną. Iteracja Gaussa-Seidela jest następnie podawana przezDALARA
x(k+1)=(D+L)−1b+GG−Sx(k)
a macierz iteracji to
GG−S=−(D+L)−1U.
SOR można następnie zapisać jako
x(k+1)=ω(D+ωL)−1b+GSORx(k)
gdzie
GSOR=(D+ωL)−1((1−ω)D−ωU).
Określenie współczynnika zbieżności schematu iteracyjnego sprowadza się naprawdę do ustalenia promienia widmowego tych matryc iteracyjnych. Zasadniczo jest to trudny problem, chyba że wiesz coś konkretnego o strukturze matrycy. Jest bardzo niewiele przykładów, o których wiem, gdzie można obliczyć optymalny współczynnik ważenia. W praktyce musi być ustalana w locie w oparciu o zaobserwowaną (przypuszczalnie) zbieżność działającego algorytmu. To działa w niektórych przypadkach, ale w innych nie.ω
Optymalny SOR
Jeden realistyczny przykład, w którym znany jest optymalny współczynnik wagowy, pojawia się w kontekście rozwiązywania równania Poissona:
∇2u=f in Ωu=g on ∂Ω
Dyskretyzacja tego systemu w dziedzinie kwadratu w 2D przy użyciu różnic skończonych drugiego rzędu z równomiernymi odstępami siatki daje w wyniku symetryczną macierz pasmową z 4 na przekątnej, -1 bezpośrednio powyżej i poniżej przekątnej oraz dwa kolejne pasma -1 pewnej odległości od przekątna. Istnieją pewne różnice wynikające z warunków brzegowych, ale taka jest podstawowa struktura. Biorąc pod uwagę tę macierz, możliwy do udowodnienia optymalny wybór współczynnika SOR daje
ω=21+sin(πΔx/L)
gdzie to odstępy między osiami, a to rozmiar domeny. Robiąc to dla prostego przypadku ze znanym rozwiązaniem daje następujący błąd w stosunku do liczby iteracji dla tych dwóch metod:ΔxL
Jak widać, SOR osiąga precyzję maszyny w około 100 iteracjach, w których Gauss-Seidel jest o około 25 rzędów wielkości gorszy. Jeśli chcesz się pobawić tym przykładem, zamieściłem kod MATLAB, którego użyłem poniżej.
clear all
close all
%number of iterations:
niter = 150;
%number of grid points in each direction
N = 16;
% [x y] = ndgrid(linspace(0,1,N),linspace(0,1,N));
[x y] = ndgrid(linspace(-pi,pi,N),linspace(-pi,pi,N));
dx = x(2,1)-x(1,1);
L = x(N,1)-x(1,1);
%desired solution:
U = sin(x/2).*cos(y);
% Right hand side for the Poisson equation (computed from U to produce the
% desired known solution)
Ix = 2:N-1;
Iy = 2:N-1;
f = zeros(size(U));
f(Ix,Iy) = (-4*U(Ix,Iy)+U(Ix-1,Iy)+U(Ix+1,Iy)+U(Ix,Iy-1)+U(Ix,Iy+1));
figure(1)
clf
contourf(x,y,U,50,'linestyle','none')
title('True solution')
%initial guess (must match boundary conditions)
U0 = U;
U0(Ix,Iy) = rand(N-2);
%Gauss-Seidel iteration:
UGS = U0; EGS = zeros(1,niter);
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
UGS(ix,iy) = -1/4*(f(ix,iy)-UGS(ix-1,iy)-UGS(ix+1,iy)-UGS(ix,iy-1)-UGS(ix,iy+1));
end
end
%error:
EGS(iter) = sum(sum((U-UGS).^2))/sum(sum(U.^2));
end
figure(2)
clf
contourf(x,y,UGS,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow
%SOR iteration:
USOR = U0; ESOR = zeros(1,niter);
w = 2/(1+sin(pi*dx/L));
for iter=1:niter
for iy=2:N-1
for ix=2:N-1
USOR(ix,iy) = (1-w)*USOR(ix,iy)-w/4*(f(ix,iy)-USOR(ix-1,iy)-USOR(ix+1,iy)-USOR(ix,iy-1)-USOR(ix,iy+1));
end
end
%error:
ESOR(iter) = sum(sum((U-USOR).^2))/sum(sum(U.^2));
end
figure(4)
clf
contourf(x,y,USOR,50,'linestyle','none')
title(sprintf('Gauss-Seidel approximate solution, iteration %d', iter))
drawnow
figure(5)
clf
semilogy(EGS,'b')
hold on
semilogy(ESOR,'r')
title('L2 relative error')
xlabel('Iteration number')
legend('Gauss-Seidel','SOR','location','southwest')