Die Konvergenz klassischer iterativer Löser für lineare Systeme wird durch den Spektralradius der Iterationsmatrix . Für ein allgemeines lineares System ist es schwierig, einen optimalen (oder sogar guten) SOR-Parameter zu bestimmen, da es schwierig ist, den Spektralradius der Iterationsmatrix zu bestimmen. Im Folgenden habe ich viele zusätzliche Details aufgeführt, einschließlich eines Beispiels für ein echtes Problem, bei dem das optimale SOR-Gewicht bekannt ist.ρ ( G )
Spektralradius und Konvergenz
Der Spektralradius ist definiert als der Absolutwert des Eigenwerts der größten Größe. Eine Methode konvergiert, wenn und ein kleinerer Spektralradius eine schnellere Konvergenz bedeutet. SOR ändert die Matrixaufteilung, die zum Ableiten der Iterationsmatrix verwendet wird, basierend auf der Wahl eines Gewichtungsparameters ω , wodurch hoffentlich der Spektralradius der resultierenden Iterationsmatrix verringert wird.ρ < 1ω
Matrixaufteilung
Für die folgende Diskussion gehe ich davon aus, dass das zu lösende System gegeben ist durch
A x = b ,
mit einer Iteration des Formulars
x( k + 1 )= v + G x( k ),
wobei ein Vektor ist und die Iterationszahl k mit x ( k ) bezeichnet wird .vkx( k)
SOR nimmt einen gewichteten Durchschnitt der alten Iteration und eine Gauß-Seidel-Iteration. Die Gauß-Seidel-Methode beruht auf einer Matrixaufteilung der Form
A = D + L + U.
wobei die Diagonale von A ist , L eine untere Dreiecksmatrix ist, die alle Elemente von A streng unterhalb der Diagonale enthält, und R eine obere Dreiecksmatrix ist, die alle Elemente von A genau oberhalb der Diagonale enthält. Die Gauß-Seidel-Iteration ist dann gegeben durchD.EINL.EINR.EIN
x( k + 1 )= ( D + L )- 1b + G.G - S.x( k )
und die Iterationsmatrix ist
GG - S.= - ( D + L )- 1U .
SOR kann dann geschrieben werden als
x( k + 1 )= ω ( D + ω L )- 1b + G.S O R.x( k )
wo
GS O R.= ( D + ω L )- 1( ( 1 - ω ) D - ω U ) .
ω
Optimale SOR
Ein realistisches Beispiel, bei dem der optimale Gewichtungskoeffizient bekannt ist, ergibt sich im Zusammenhang mit der Lösung einer Poisson-Gleichung:
∇2u = f i n Ω u = g o n ∂ Ω
Die Diskretisierung dieses Systems auf einer quadratischen Domäne in 2D unter Verwendung endlicher Differenzen zweiter Ordnung mit gleichmäßigem Gitterabstand führt zu einer symmetrischen Bandmatrix mit 4 auf der Diagonale, -1 unmittelbar über und unter der Diagonale und zwei weiteren Bändern von -1 in einiger Entfernung von der diagonal. Aufgrund der Randbedingungen gibt es einige Unterschiede, aber das ist die Grundstruktur. In Anbetracht dieser Matrix ist die nachweislich optimale Wahl für den SOR-Koeffizienten gegeben durch
ω = 21 + Sünde( πΔ x / L )
Δ xL.
Wie Sie sehen können, erreicht SOR die Maschinengenauigkeit in etwa 100 Iterationen, wobei Gauß-Seidel um etwa 25 Größenordnungen schlechter ist. Wenn Sie mit diesem Beispiel herumspielen möchten, habe ich den unten verwendeten MATLAB-Code eingefügt.
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')
Diese Seite der Dinge ist nicht wirklich meine Spezialität, aber ich denke nicht, dass dies ein super fairer Test für viele realistische Anwendungen ist.
Ich bin nicht sicher, welche Werte Sie für c und r verwendet haben , aber ich vermute, Sie haben mit extrem schlecht konditionierten Matrizen gearbeitet. (Nachfolgend finden Sie Python-Code, der zeigt, dass dies möglicherweise nicht die invertierbarsten Matrizen sind.)
Wenn Sie tatsächlich Matrizen invertieren müssten, die so schlecht konditioniert sind, würden Sie a) eine spezielle Methode verwenden und b) wahrscheinlich nur ein neues Feld suchen 😉
Bei gut konditionierten Matrizen jeder Größe ist SOR wahrscheinlich schneller. Für echte Probleme, bei denen es auf Geschwindigkeit ankommt, ist es selten, SOR zu verwenden - auf der anspruchsvollen Seite gibt es heutzutage viel bessere; Auf der langsamen, aber zuverlässigen Seite ist SOR nicht das Beste, was Sie tun können.
quelle
OK, also für symetrische Matrizen dieses Königs:
(Dies ist nur eine kaiserliche Beobachtung, nichts Strenges)
quelle