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')
Ta strona rzeczy nie jest tak naprawdę moją specjalnością, ale nie sądzę, że jest to bardzo uczciwy test dla wielu realistycznych aplikacji.
Nie jestem pewien, co ceni uzywasz dla C i R , ale podejrzewam, że pracujesz z bardzo źle uwarunkowanych matryc. (Poniżej znajduje się kod w języku Python pokazujący, że mogą to nie być najbardziej odwracalne macierze).
Jeśli faktycznie trzeba było odwracać macierze tak źle uwarunkowane, a) należy użyć specjalistycznej metody ib) prawdopodobnie należy po prostu znaleźć nowe pole 😉
W przypadku dobrze kondycjonowanych matryc o dowolnym rozmiarze SOR prawdopodobnie będzie szybszy. W przypadku prawdziwych problemów, w których liczy się prędkość, rzadko byłoby używać SOR - po wyrafinowanej stronie jest obecnie znacznie lepiej; powolna, ale niezawodna strona, SOR nie jest najlepszym rozwiązaniem.
źródło
OK, więc dla matryc symetrycznych tego króla:
SOR zbiega się szybciej niż Gaussa-Seidela, jeśli liczba w każdym rzędzie jest mała (znacznie mniejsza niż wymiar A) i jeśli wszystkie są podobne. Używałem wygenerowanych w ten sposób:t t t
Jeśli s różnią się znacznie i są wyśrodkowane wokół 0 ( ), to Gauss-Seidel jest szybszy. Gauss-Seidel jest także szybszy, jeśli każdy rząd jest więcej niż w połowie wypełniony przez . Oznacza to również, że SOR jest lepszy dla bardzo dużych i bardzo rzadkich matryc.t c=0,r=0.1 t
(To tylko obserwacja cesarska, nic rygorystycznego)
źródło