The Convergence of SOR Relaxation for the Poisson Equation

Allan Struthers, Math Sciences, Michigan Tech

Large sparse linear system that arises from finite difference discretizations of Poisson' s equation are standard test problems for iterative linear solvers with SOR (Successive Over Relaxation) - using optimal (or near optimal) values of the relaxation parameter - frequently used as a benchmark.This presentation will describe some surprising behavior (and the linear algebra explanation for it) observed when developing such a benchmark.

Introduction

I generated the following pretty picture.  In some sense it should look like a bunch of straight lines with one slightly lower than the others. It does not!

This is the story of the picture and the explanation.

Poisson's Equation

Our 1D problem on  ax≤x≤bx with either periodic u(ax)=u(bx)  or homogeneous Dirichlet conditions u(ax)=u(bx)=0 is

(∂_xx)u =f(x)

For the periodic conditions we need the compatibility condition that (∫_ax^bx)f(x)dx=0.

Our  2D problem on  ax≤x≤bx and  ay ≤ y ≤ by with either periodic u(ax,y)=u(bx,y)  and u(x, ay) = u(x,by) or homogeneos Dirichlet conditions u(ax,y)=u(bx,y)=u(x,ay)=u(x,by)=0 is

(∂_xx)u + (∂_yy)u =f(x,y)

Finite Differences 1D

Centered differences gives a standard approximation for the second partial derivative as follows

(∂_xx)u ≈((u[x-Δx)-2U(x)+U(x+Δx))/((Δx)^2))

Choosing a uniform gird with nx interior points gives Δx = (bx - ax)/(nx + 1) and the x values (for i=0,1,…,nx,nx+1)

(Subscript[x, i])=a+i Δx

We get the following system of equations ( i=1,2,…,nx)

(Subscript[u, i-1])-2(Subscript[u, i])+(Subscript[u, i+1])=((Δx)^2)f((Subscript[x, i]))

which almost define the standard finitie difference approximation (Subscript[u, i])≈u((Subscript[x, i])).  Of course this is nx equations in nx+2 unknowns. The boundary conditions supply the other two equations: 
a) For dirichlet conditions eliminate (Subscript[u, 0])=0 and (Subscript[u, nx+1])=0 giving an nx by nx system.
b) For periodic conditions eliminate (Subscript[u, xn+1]) (it is equal to (Subscript[u, 0])) and append the equation 
	(Subscript[u, nx])-2(Subscript[u, 0])+(Subscript[u, o+1])=0.5*((Δx)^2)*(f((Subscript[x, 0]))+f((Subscript[x, nx])))

Finite Difference Matrices 1D

It is a good idea to write these as standard linear algebra A.(OverVector[u])=(OverVector[b]) equations for the vector if unknowns (Subscript[u, i]). For Dirichlet conditions the matrix (with 6 interior points) is

A=(({{-2, 1, 0, 0, 0, 0}, {1, -2, 1, 0, 0, 0}, {0, 1, -2, 1, 0, 0}, {0, 0, 1, -2, 1, 0}, {0, 0, 0, 1, -2, 1}, {0, 0, 0, 0, 1, -2}}))

For periodic conditions the singular matrix (with 6 interior points) is

A=(({{-2, 1, 0, 0, 0, 0, 1}, {1, -2, 1, 0, 0, 0, 0}, {0, 1, -2, 1, 0, 0, 0}, {0, 0, 1, -2, 1, 0, 0}, {0, 0, 0, 1, -2, 1, 0}, {0, 0, 0, 0, 1, -2, 1}, {1, 0, 0, 0, 0, 0, 1}}))

Right Hand Sides 1D

For Dirichlet conditions the RHS vector (OverVector[b]) is

(OverVector[b])=((Δx)^2)(({{f((Subscript[x, 1]))}, {f((Subscript[x, 2]))}, {f((Subscript[x, 3]))}, {f((Subscript[x, 4]))}, {f((Subscript[x, 5]))}, {f((Subscript[x, 6]))}}))

For periodic conditions the singular matrix (with 6 interior points) is

(OverVector[b])=((Δx)^2)(({{0.5*(f((Subscript[x, 0]))+f((Subscript[x, 7])))}, {f((Subscript[x, 1]))}, {f((Subscript[x, 2]))}, {f((Subscript[x, 3]))}, {f((Subscript[x, 4]))}, {f((Subscript[x, 5]))}, {f((Subscript[x, 6]))}}))

The 1D null space of the periodic A is the span of (OverVector[1])={1,1,1,1,1,1,1}.  The consistency condition is that (OverVector[b])(OverVector[.1])=0.

Jacobi Iteration 1D

Solving (5) for (Subscript[u, i]) gives

(Subscript[u, i])=0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i]))

and interpreting the result as an iteration gives a Jacobi iteration for the linear system which can be written as

((OverVector[u])^(n+1))=G.((OverVector[u])^n)-0.5((Δx)^2)(OverVector[f])

For Dirichlet conditions (with 6 interior points) the iteration matrix is

G=0.5*(({{0, 1, 0, 0, 0, 0}, {1, 0, 1, 0, 0, 0}, {0, 1, 0, 1, 0, 0}, {0, 0, 1, 0, 1, 0}, {0, 0, 0, 1, 0, 1}, {0, 0, 0, 0, 1, 0}}))

For periodic conditions (with 6 interior points) the iteration matrix is

G=0.5*(({{0, 1, 0, 0, 0, 0, 1}, {1, 0, 1, 0, 0, 0, 0}, {0, 1, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 1, 0, 0}, {0, 0, 0, 1, 0, 1, 0}, {0, 0, 0, 0, 1, 0, 1}, {1, 0, 0, 0, 0, 1, 0}}))

Convergence 1D

The iteration converges iff the eigenvalues of  G are less than one in magnitude.  Writing the {(Subscript[λ, i]),(Subscript[(OverVector[v]), i])} for the eigenvalues-eigenvector pairs and ((OverVector[u])^∞) for the solution of the linear system we have the generic asymptotic convergence theorem

((OverVector[u])^n)-((OverVector[u])^∞)=∑(Subscript[c, i])(Subscript[λ, i]^n)(Subscript[(OverVector[v]), i])

for some constants (Subscript[c, i]) determined by the initial guess of ((OverVector[u])^0).  The standard asymptotic theory use only the largest single eigenvalue.

Gauss-Seidel Iteration: Red/Black

The Jacobi iteration is naturally implemented storing both an old and a new (OverVector[u]) vector.  A Gauss-Seidel or relaxation iteration overwrites the old (OverVector[u]) vector with the new values and uses new values as they become avaliable.  For example,

Do[
	(Subscript[u, i])=0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i])),
	i…2,3,…,dim-1]

(with appropriate updates for the boundary points) is a Gauss-Seidel iteration.  Such iterations clearly depends on the order in which the variables are processed.  However, our equation partitions neatly into even and odd values (called red/black in 2 and higher dimensions) which allows us to continue the analysis provided we do all the odds then all the evens. The following is a step of the red/black ordered GS iteration

Do[
	(Subscript[u, i])=0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i])),
	i:odd]
	Do[
	(Subscript[u, i])=0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i])),
	i:even]

Our final iteration matrix is G=(Subscript[G, b])(ω).(Subscript[G, e])(ω)

Succesive Over Relaxation

The idea behind SOR is "If a little bit of the Gauss-Seidel iteration is good more might be better."   In other words inflate the step by some factor ω and weight it backoff with a 1-ω.  Here it is.

Do[
	(Subscript[u, i])=ω ⟨0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i]))⟩ + (1-ω) (Subscript[u, i]),
	i:odd]
	Do[
	(Subscript[u, i])=ω ⟨0.5*((Subscript[u, i+1])+(Subscript[u, i-1]))-0.5((Δx)^2)f((Subscript[x, i]))⟩ + (1-ω) (Subscript[u, i]),
	i:even]

Now, we can analyze this because the update matrix is just the product of two matrices.   For

For the Dirichlet problem (in our dimension) we have

(Subscript[G, r])(ω)=(({{1-ω, (ω/2), 0, 0, 0}, {0, 1, 0, 0, 0}, {0, (ω/2), 1-ω, (ω/2), 0}, {0, 0, 0, 1, 0}, {0, 0, 0, (ω/2), 1-ω}}))
	(Subscript[G, b])(ω)=(({{1, 0, 0, 0, 0}, {(ω/2), 1-ω, (ω/2), 0, 0}, {0, 0, 1, 0, 0}, {0, 0, (ω/2), 1-ω, (ω/2)}, {0, 0, 0, 0, 1}}))

Our final iteration matrix is G=(Subscript[G, b])(ω).(Subscript[G, e])(ω) and if we want GS we just put ω=1.

Optimal Asymptotic Convergence

The eigenvalues and eigenvectors of both of our iteration matrices are available analytically. For simplicity we will illustrate the results numerically for the periodic conditions.

GE[ω_,n_]:=SparseArray[{
Band[{1,1},{n,n}]->{1,(1-ω)},
Band[{1,2},{n,n}]->{0,ω/2},
Band[{2,1},{n,n}]->{ω/2,0}
(* {n,1}->ω/2 *)
},
{n,n}]
GO[ω_,n_]:=SparseArray[{
(* {1,n}->ω/2, *)
Band[{1,1},{n,n}]->{(1-ω),1},
Band[{1,2},{n,n}]->{ω/2,0},
Band[{2,1},{n,n}]->{0,ω/2}
},
{n,n}]

Here is a plot of the largest eigen-value in absolute value with a clear minimum.

n=16;
Manipulate[
Plot[
Abs[Eigenvalues[GE[ω,n].GO[ω,n],1][[1]]],
{ω,0,2},
PlotRange->{0.6,1},
PlotLabel->StringForm["n=``",n]
],
{n,12,24,1}]

Relaxation Iteration 1D: Results

Code for Red/Black 1D SOR

n=125;
{a,b}={0.0,1.0};
h=((b-a)/n);
𝒜= SparseArray[{
Band[{1,1}]->2.0,
Band[{2,1}]->-1.0,
Band[{1,2}]->-1.0
},
{n-1,n-1}
];
f[x_]= 3 Sin[x] + 3 Cos[x]- 4 Cos[2 x]-Sin[6 x];
𝒻=Table[f[x],{x,a+h,b-h,h}];
𝓊= LinearSolve[𝒜,(h^2)𝒻];

SetAttributes[SORUpDate1DRedBlack, HoldFirst]
SORUpDate1DRedBlack[u_,𝒻_,h_,ω_]:=Module[
{i,n=Length[𝒻]},
u[[1]]=(1-ω)u[[1]]+(ω/2)((h^2)𝒻[[1]]+u[[2]]);
u[[n]]=(1-ω)u[[n]]+(ω/2)((h^2)𝒻[[n]]+u[[n-1]]);
Do[
u[[i]]=(1-ω)u[[i]]+(ω/2)((h^2)𝒻[[i]]+u[[i-1]]+u[[i+1]]),
{i,3, n-1,2}];
Do[
u[[i]]=(1-ω)u[[i]]+(ω/2)((h^2)𝒻[[i]]+u[[i-1]]+u[[i+1]]),
{i,2,n-1,2}];
]

Results

If you code it and run it you get the following.

Clear[DataSORrb]
{ωMin,ωMax,Δω}={1.94, 2.0,0.00125};
NMax=400;
DataSORrb=Table[
u=0*𝒻;
ScaleNorm=Norm[𝓊];
Table[
SORUpDate1DRedBlack[u,𝒻,h,ω];
Norm[u-𝓊,2]/ScaleNorm,
{NMax}],
{ω,ωMin, ωMax-Δω,Δω}];

ListContourPlot[Log[10,DataSORrb],
Contours->20,
DataRange->{{1,NMax},{ωMin,ωMax-Δω}}]

Relaxation Iteration 1D: More Results

Focusing in on the interesting bit at the top!

Clear[DataSORrb]
{ωMin,ωMax,Δω}={1.94, 2.0,0.00125};
NMax=100;
DataSORrb=Table[
u=0*𝒻;
ScaleNorm=Norm[𝓊];
Table[
SORUpDate1DRedBlack[u,𝒻,h,ω];
Norm[u-𝓊,2]/ScaleNorm,
{NMax}],
{ω,ωMin, ωMax-Δω,Δω}];

ListLogPlot[ Take[DataSORrb, -29],
PlotStyle->Table[Hue[i/25],{i, 0, 24}],
GridLines->Automatic,
Frame->True,
FrameLabel->{"n","relative residual"}]

Analytical Explanation

The exlanation is simple.  The asymptotic results are correct in (ℂ^n) but we are running the computation (and everybody does) in (ℝ^n).  For ω near 2 the eigenvalues are complex conjugate pairs and since the computation is in reals there is an exponential decay and a periodic portion.  That periodic portion gives the spikes.  The amplification associated with a complex conjugate eigenpair {ρ (e^(±i θ)),(OverVector[vr])±i (OverVector[vi])} over n iterations is

(ρ^n)||(OverVector[vr]) Cos[n θ]- (OverVector[vi])Sin[n θ]||

This is easy enough to implement and plot.

AmpFactors[ω_,{nd_,ni_}]:= Module[
{λs,Qs,i},
{λs,Qs}=Eigensystem[Normal[GO[ω,nd].GE[ω,nd]]];
Table[
λ=λs[[i]]; v=Qs[[i]];
ρ=Abs[λ];θ=Arg[λ];
vr=Re[v]; vi = Im[v];
(ρ^ni) Norm[ vr Cos[θ ni]-vi Sin[θ ni]],{i,1,nd}]
];
MaxAmp[ω_,{nd_,ni_}]:=Max[ AmpFactors[ω,{nd,ni}]];

{dim,nIter}={14,14};
LogPlot[{
MaxAmp[ω,{dim,nIter}],
((Abs[Eigenvalues[GO[ω,dim].GE[ω,dim],1]])^nIter)
},
{ω,1.5, 1.8},
PlotLabel->StringForm["dimension=`` and iteration=``",dim,nIter]
]

Optimality Revisted

The optimal value of ω for a fixed number of iterations is slightly greater than the spectral radius estimate.

This provides an additional explanation for the accepted wisdom when talking about the .  "It is better to over estimate the ω than under estimate it!"

Further Work

1) The same features persist in 2 and 3D dimensions.  The analysis needs to be extended and the magnitude of the effect quantified. plotted.

2) It is unknown if similar features exist for other discretizations and other linear systems.

References

[1] Applied Numerical Linear Algebra, James W. Demmel, SIAM, Philadelphia, 1997.

[2] Iterative Methods for Solving Linear Systems, Anne Greenbaum, SIAM, Philadelphia, 1997.

[2]The Optimal rleaxation Parameter for the SOR Method Applied to the Poisson Equation in Any Space Dimension, Shiming Yang and Mathiass Gobbert, Applied Mathematics Letters, Volume 22, pp325-331, 2009