CPS615—Base Course

for the Simulation Track

of Computational Science

Fall Semester 1995 Finite Element Methods

and Conjugate Gradient Methods

Geoffrey Fox

NPAC Room 3-131 CST 111 College Place

Syracuse, New York 13244-4100

CPS615—Abstract

Finite Element/Conjugate Gradient Presentation

• This derives the finite element method for a simple two-dimensional Laplacian with triangular elements

• We use this to motivate conjugate gradient as a variant of steepest descent for variational principle underlying FEM

• We discuss preconditioning, parallelism and convergence of general conjugate gradient method

Integral Formulation of Finite Element Method

• Consider

F = \frac{1}{2} ∫ _{Ω} (∇φ)^{2} d^{3}\underline{x}

where φ is defined over a region Ω in d dimensional space — d = 2 in text as one can draw.

• φ is fixed = f on boundary ∂Ω of Ω (Figure FEM 1)

Figure: FEM 1

Variation in Integral

• Consider variation Let φ → φ + η with η = 0 on ∂Ω and η a function of \underline{x} over d dimensional Ω

• Then this change in φ induces a change F → F + δF where, to first order in η,

δF ≈ ∫_{Ω} (∇φ · ∇η) d^{3}\underline{x}

• We have Gauss’s theorem ∫_{Ω} ∇\underline{u}d^{3}x = ∫_{∂Ω} \underline{u} · d\underline{s}

where d\underline{s} is differential on surface ∂Ω

• apply to \underline{u} = η∇φ ∇\underline{u} = ∇η · ∇φ + η∇^{2}φ

where η = 0 on ∂Ω

• Then

∫_{Ω} ∇φ ∇ηd^{3}\underline{x} = − ∫ η∇^{2}φd^{3}\underline{x}

relates δF to ∇^{2}φ where

∇φ = (\frac{∂φ}{∂x}, \frac{∂φ}{∂y}, \frac{∂φ}{∂z})

and

(∇φ)^{2} = (\frac{∂φ}{∂x})^{2} + (\frac{∂φ}{∂y})^{2} + (\frac{∂φ}{∂z})^{2}

Equivalence of Integral and Differential Formulation of Laplace’s Equation

• Therefore, if F is stationary, i.e.

\frac{∂F}{∂ (η(\underline{x}))} = 0

for all η(\underline{x}) consistent with boundary conditions, then ∇^{2}φ = 0

• Thus, Laplace’s equation is entirely equivalent to VARIATIONAL PRINCIPLE

δF = 0
F = ∫_{Ω} \frac{1}{2} (∇φ)^{2} d^{3}\underline{x}
φ(\underline{x}) = known f(x) on ∂Ω

We have converted differential equation to an integral …

Discretization of Integral

• Convert this to a minimization principle over a finite dimensional space defining φ

• Consider those functions φ which are piecewise linear on a set of triangles (the mesh) which covers Ω (Figure FEM 2)

Figure: FEM 2

• i.e., we will use the generalized trapezoidal rule to integrate F

Triangular Elements in Two Dimensions

• First consider evaluating F on a triangle with a linear interpolation function (Figure FEM 3)

Figure: FEM 3

• so that

φ(\underline{r}) = φ_{0} + (φ_{1} − φ_{0}) \frac{(\underline{r} − \underline{r}_{0}) ∧ (\underline{r}_{2} − \underline{r}_{0}) |_{z}}{(\underline{r}_{1} − \underline{r}_{0}) ∧ (\underline{r}_{2} − \underline{r}_{0}) |_{z}}

+ (φ_{2} − φ_{0}) \frac{(\underline{r} − \underline{r}_{0}) ∧ (\underline{r}_{1} − \underline{r}_{0}) |_{z}}{(\underline{r}_{2} − \underline{r}_{0}) ∧ (\underline{r}_{1} − \underline{r}_{0}) |_{z}}

• It can be seen that φ is linear in \underline{r} and has value φ_{0} at \underline{r}_{0}, φ_{1} at \underline{r}_{1}, and φ_{2} at \underline{r}_{2}

Example for Two-Dimensional Triangular Elements

• Show that integrating over one triangle gives

F = ∫_{Δ} \frac{1}{2} (∇φ)^{2} d^{2}x
= \frac{1}{4} [φ_{0} φ_{1} φ_{2}] [ \begin{matrix}c_{1} + c_{2} & −c_{2} & −c_{1} \\ −c_{2} & c_{2} + c_{0} & −c_{0} \\ −c_{1} & −c_{0} & c_{0} + c_{1} \end{matrix} ] [ \begin{matrix}φ_{0} \\ φ_{1} \\ φ_{2} \end{matrix} ]

where c_{j} = cot θ_{j}

Bilinear Form of Integral with Triangular Elements

• Now consider a triangulation of Ω (Figure FEM 4)

Figure: FEM 4

• We can write

F = \frac{1}{4}\underline{φ}^{T}S\underline{φ}

where \underline{φ} is the vector of vertex values of φ and S is called a stiffness matrix

CPS615FEM, G.C.Fox, Oct 95  9  gcf@npac.syr.edu

Formula for Stiffness Matrix Element—I

Figure: FEM 5

• S_{i}i = ∑_{angles θ} cot θ where angles θ are opposite nodal point i in a triangle with a vertex at i.

CPS615FEM, G.C.Fox, Oct 95  10  gcf@npac.syr.edu

Formula for Stiffness Matrix Element—II

Figure: FEM 6

• If i and j are not connected in the triangulation, then S_{ij} = 0. If they are connected, written i ↔ j and i ≠ j then S_{ij} = − (cot θ_{1} + cot θ_{2}) (Figure FEM 5)

• Example: Show that

_{i}^{∑} S_{ii} + _{i≠j}^{∑} S_{ij} = 0

Hint: φ(\underline{x}) = 1 is a solution to Laplace’s equation for which approximations are exact

Finite Element Equations

• We now have

F = \frac{1}{4}φ^{T}Sφ
= \frac{1}{4} _{ij}^{∑} φ_{j}S_{ji}φ_{i}

• Let

φ_{i} → φ_{i} + η_{i} be variation
F = \frac{1}{4} _{ij}^{∑} φ_{j}S_{ji}φ_{i}
+ \frac{1}{4} (_{ij}^{∑} (η_{j}S_{ji}φ_{i}) + (φ_{j}S_{ji}η_{i}))

• Thus δF = \frac{1}{2} ∑_{ij} η_{j} S_{ji} φ_{i}

• This δF form is zero if and only if _{i}^{∑} S_{ji} φ_{i} = 0

• This condition can be split into internal and boundary terms

\underbrace{Σ_{(j not∈ ∂Ω)}S_{ij}φ_{j}}_{Unknown} + \underbrace{Σ_{(j∈∂Ω)}S_{ij} (φ = f_{j})}_{Known at Boundary} = 0

CPS615FEM, G.C.Fox, Oct 95  12  gcf@npac.syr.edu

Structure of Stiffness Matrix and Its Assembly

• N has contributions from five triangles in picture (Figure FEM 7)

• This is called “assembly of the stiffness matrix S”

• Note by using variational formalism, we have avoided need for regular grid

Figure: FEM 7

Conditions on Triangulation

• You need nicely shaped triangles to ensure that

_{j≠i}^{∑} |S_{ij}| ≤ |S_{ii}|

for each i (diagonal dominance)

• This is satisfied if triangulation is such that θ_{1} + θ_{2} ≤ 180^{∘} for every internal edge e (Figure FEM 8a).

Figure: FEM 8a

• Skinny triangles violate this!

Figure: FEM 8b

Introduction to Poor Person’s Conjugate Gradient

• Finite Element Method leads to minimization of a Quadratic Form

F (\underline{x}) = \frac{1}{2}\underline{x}^{T}A\underline{x} − \underline{x}^{T}b  (1)

with customary solution

A\underline{x} = \underline{b}  (2)

derived by differentiating (1) with respect to \underline{x} or by completing the square

F (\underline{x}) = \frac{1}{2} (∇ − A^{−1}\underline{b})^{T} \underline{A} (\underline{x} − A^{−1}\underline{b}) − \frac{1}{2}\underline{b}^{T}A^{−1}\underline{b}

• We will derive an iterative approach (PPCG) to

(2) by an intuitive minimization method applied to (1). This is not the true conjugate method, but has the basic idea. As usual, iterations are labelled by (k) and we start with \underline{x}^{(k)}

Conjugate Gradient Iteration for Quadratic Form

• \underline{x}^{(k)} is a point in an n dimension space (n is number of equations in A\underline{x} = \underline{b})

−\frac{∂F}{∂x} = −∇F = \underline{b} − A\underline{x}^{(k)} = \underline{r}^{(k)}

• \underline{r}^{(k)} is residual

• Let \underline{x}^{(k+1)} = \underline{x}^{(k)} + α\underline{r}^{(k)} (α a scalar, i.e., a number).

• F becomes a function of α for this special choice of \underline{x}^{(k+1)}. This is useful as it is much easier to minimize functions of one variable than those that depend on n.

F (α) = α^{2} [\frac{1}{2} \underline{r}^{(k)T}A\underline{r}^{(k)}]
− α\underline{r}^{(k)T} · \underline{r}^{(k)}
+ Terms Independent of α

Conjugate Gradient and Method of Steepest Descent

• We can complete the square

F (α) = \frac{1}{2}\underline{r}^{(k)T}A\underline{r}^{(k)} (α − α_{k})^{2}
+ Terms Independent of α

• where α_{k} = \underline{r}^{(k)T} · \underline{r}^{(k)}/\underline{r}^{(k)T}A\underline{r}^{(k)}

• Note that if

\underline{x}^{(k+1)} = \underline{x}^{(k)} + α_{k}\underline{r}^{(k)}
F (\underline{x}^{(k+1)}) ≤ F (\underline{x}^{(k)}) rigorously

• This is Poor Person’s Conjugate Gradient or when applied to minimization, the method of steepest descent. It is called this because we move from \underline{x}^{(k)} along direction −∇F —the direction along which F varies most rapidly.

Conjugate Gradient for Finite Element Problems

• Conjugate Gradient can be applied to any matrix problem.

• It is associated naturally with finite element problems because both finite element and conjugate gradient use similar minimization procedures.

Poor Person’s Conjugate Gradient and Eigenvalues of Matrix

• What is wrong with this version (PPCG) of conjugate gradient?

• More generally, what is wrong with steepest descent?

• The matrix A is symmetric and can be diagonalized as usual

A = P^{−1}ΛP
Λ = diag. (λ_{1} , λ_{2} , … λ_{n})
= [ \begin{matrix}λ_{1} & 0 & 0 & … \\ 0 & λ_{2} & 0 & … \\ 0 & 0 & λ_{3} & … \end{matrix} ]
λ_{1} ≥ λ_{2} ≥ …

• P is orthogonal matrix P^{T}P = I

• Let

\underline{x̃}^{(k)} = P \underline{x}^{(k)}
\underline{b̃} = P \underline{b}
\underline{r̃}^{(k)} = P \underline{r}^{(k)}

Diagonalization of Quadratic Form

• Note that for any vector \underline{x}

\frac{1}{2}\underline{x}^{T}A\underline{x} − \underline{b}^{T} · \underline{x}

= \frac{1}{2} \underbrace{\underline{x}^{T}P^{−1}}_{x^{T} P^{T}} ΛP \underline{x} − \underbrace{\underline{b}^{T}P^{−1}}_{(P\underline{b})^{T}} P \underline{x}

= \frac{1}{2}(P \underline{x})^{T}ΛP \underline{x} − (P \underline{b})^{T}P \underline{x}
= \frac{1}{2}\underline{x̃}^{T}Λ\underline{x̃} − b̃^{T} \underline{x̃}

• Note this transformation is correct but not useful as a solution technique as it takes longer to find Λ and P than to solve equations A\underline{x} = \underline{b}.

• In transformed space denoted by superscript

x̃|_{i compt} = b̃_{i}/λ_{i}

CPS615FEM, G.C.Fox, Oct 95  20  gcf@npac.syr.edu

Diagonalization of Conjugate Gradient Equations

• We can take all the PPCG equations and multiply them by the transformation matrix P and convert them into “ ” equations.

• Scale factor α

α^{(k)} = \frac{\underline{r̃}^{(k−1)T} · \underline{r̃}^{(k−1)}}{\underline{r̃}^{(k−1)T}Λ\underline{r̃}^{(k−1)}}

• Iteration

\underline{x̃}^{(k)} = \underline{x̃}^{(k−1)} + α_{k}\underline{r̃}^{(k−1)}

• Residual

\underline{r̃}^{(k)} = \underline{b̃} − Λ \underline{x̃}^{(k)}
= \underline{b̃} − Λ \underline{x̃}^{(k−1)} − α_{k}Λ\underline{r̃}^{(k−1)}
= (I − α_{k}Λ) r̃^{(k−1)}

• Residual components in diagonalized coordinate system

r̃^{(k)}|_{i} = (1 − α_{k}λ_{i}) r̃^{(k−1)}|_{i}

where i = 1 … n labels components

Convergence of Conjugate Gradient in Diagonalized Form

• We want r̃_{k}|_{i} = 0 for all i

• But

α_{k} = \frac{∑_{i} [(r̃^{(k−1)}) |_{i}]^{2}}{∑_{i} λ_{i} [r̃^{(k−1)}|_{i}]^{2}}

• We order eigenvalues λ_{1} ≥ λ_{2} ≥ …

• If r̃^{(k−1)}|_{1} non zero and reasonably big

• Then

α_{k} ∼ \frac{1}{λ_{1}} and
r̃^{(k)}|_{1} = (1 − \frac{λ_{1}}{λ_{1}}) r̃^{(k−1)}|_{1}
∼ 0

• But other r̃_{k}|_{i} i > 1 are largely unchanged if eigenvalues λ_{i} ≪ λ_{1} because as r̃^{(k−1)}|_{i} is multiplied by (1 − λ_{i}/λ_{1}) ∼ 1 (when λ_{i} ≪ λ_{1}) on going from iteration (k − 1) to (k)

Clarification of Eigenvalue Analysis for Conjugate Gradient and Jacobi Iteration

• Note in eigenvalues: In Jacobi/SOR, “large” eigenvalues converged slowly, “small” eigenvalues converged fast.

• But the matrices aren’t the same.

• Conjugate gradient (CG) uses complete matrix A

• SOR uses decomposed matrix A = M + N

• M^{−1}N eigenvalues determine SOR A eigenvalues determine CG

• Large eigenvalues of A → small eigenvalues of SOR and vice versa.

Intuitive Description of Poor Person’s Conjugate Gradient Algorithm

• So Poor Person’s Conjugate Gradient is very good at large eigenvalues, but converges slowly for small eigenvalues of A.

• Note A is positive definite for simple pde’s ⇒ λ_{i} ≥ 0.

• What we would like to do is to arrange first that r̃^{(k)}|_{1} = 0, then we see that equations give r̃^{(k+1)}|_{2} = 0, and then the next iteration will give r̃^{(k+2)}|_{3} = 0, and so on.

• So qualitative goal of real conjugate gradient is to successively minimize along the “principal directions of matrix A.”

• As described in my optimization notes, this is easily understood in case n = 2 in terms moving in hills and valleys.

Improvement of Poor Person’s Conjugate Gradient with Orthonormal Iteration

• In practice |λ_{1} = λ_{max}| ≫ |λ_{n} = λ_{min}| and in simple cases \frac{λ_{n}}{λ_{1}} ∼ h^{2} where h is grid spacing

• We must modify PPCG to keep search directions properly orthogonal so that eigenvectors with large eigenvalues do not dominate

• We set

\underline{x}^{(k)} = \underline{x}^{(k−1)} + α_{k}\underline{p}^{(k)}

where we choose

\underline{p}^{(k)T}A\underline{p}^{(k−1)} = 0

Full Conjugate Gradient Algorithm

• Initialize in true conjugate gradient

\underline{r}^{(0)} = \underline{b} − A\underline{x}^{(0)}
k = 1  p^{(1)} = \underline{r}^{(0)}
k > 1  p^{(k)} = \underline{r}^{(k−1)} + β_{k}\underline{p}^{(k−1)}

• This subtracts off \underline{p}^{(k−1)} component in \underline{r}^{(k−1)}

• as long as:

β_{k} = −\frac{\underline{p}^{(k−1)T}A\underline{r}^{(k−1)}}{\underline{p}^{(k−1)T}A\underline{p}^{(k−1)}}

• Then as before, solve for α_{k} so that \underline{x}^{(k−1)} + α_{k}\underline{p}^{(k)} minimizes the quadratic form −1/2\underline{x}^{T}A\underline{x} − \underline{b}^{T} · \underline{x} • This gives

α_{k} = \underline{p}^{(k)T} \underline{p}^{(k)}/\underline{p}^{(k)T}A\underline{p}^{(k)}
\underline{x}^{(k)} = \underline{x}^{(k−1)} + α_{k}\underline{p}^{k}
\underline{r}^{(k)} = b − Ax^{(\underline{k})}

Overview of Parallelism in Conjugate Gradient

• One must 1. Calculate matrix elements A_{ij}

2. Calculate various:

• Matrix-vector ∑_{j} A_{ij}x_{j} products

• and vector-vector ∑_{j} x_{j}y_{j} scalar products • for several different vectors

(A is always stiffness matrix).

Parallel Issues in Calculation of Matrix Elements

• Divide world up into geometric domains as usual (Figure FEM 9)

Figure: FEM 9

• Information (triangles with i, j as an edge) to calculate A_{ij} can involve ≥ 2 processors (four is most in two-dimensional example), but communication is local.

• This is called assembly.

• Each processor should know A_{ii} and A_{ij} for all j connected to each i “owned” by processor.

Scalar Products in Parallel Conjugate Gradient

• When one calculates ∑_{j} A_{ij}x_{j}, send x_{j} to all processors “owning” x_{i} where A_{ij} ≠ 0

• This gives local communication

• When one calculates

n ; ∑ ; j=1 y_{j}x_{j}

First calculate sums in a given processor, summing in Processor P over all nodal points j in P. Then one must form

• a Global “combine”.

_{Processors P}^{∑} [_{jϵP}^{∑} y_{j}x_{j}]

Typically, this same sum must be sent to all processors.

Preconditioning in Conjugate Gradient

• This is a general technique, where one finds M such that M\underline{x} = \underline{y} can be easily solved and M^{−1}A is numerically faster (easier or more stable) to solve than A i.e., M is an approximation to A whose inverse can be found easily—more precisely, for which linear equations M\underline{x} = \underline{y} can be solved as one never inverts matrices in this type of problem.

• The simplest preconditioner is M = Diagonal part of A.

• The best solution methods use sophisticated preconditioners such as those given by multigrid.

Convergence of Conjugate Gradient

• The number of conjugate gradient iterations needed for given accuracy for simple problems (positive definite matrix A) is given by

k = \sqrt{\frac{λ_{max}}{λ_{min}}}
∝ \frac{1}{h}

where h is grid size

• Remember SOR was 1/h \underline{if} overrelaxation parameter w chosen correctly, which is not easy.

• Jacobi and Gauss Seidel have a number of iterations k ∝ 1/h^{2} which is much more

• Conjugate gradient iterations are more complicated and time consuming than those for Jacobi or Gauss Seidel

Mathematical and Pseudo Code Form of Gauss Seidel Iteration Method

Choose an initial guess x^{(0)} to the solution x. for k = 1,2,…

for i = 1,2,…,n σ = 0 for j = 1,2,…,i − 1

σ = σ + a_{i,j} x_{j}^{(k)} end for j = i + 1,…,n

σ = σ + a_{i,j} x_{j}^{(k−1)} end x_{i}^{(k)} = (b_{i} − σ) /a_{i,i}

end check convergence; continue if necessary

end

Pseudo code for Gauss-Seidel Method

Mathematical (Matrix) Form of Gauss Seidel

• Mathematically, Gauss Seidel is

x_{i}^{(k)} = \frac{[b_{i} − _{j<i}^{∑} a_{i,j} x_{j}^{(k)} − _{j>i}^{∑} a_{i,j} x_{j}^{(k−1)}]}{a_{i,i}}

• In matrix language \underline{x}^{(k)} = (D − L)^{−1} (U\underline{x}^{(k−1)} + \underline{b})

• where D=Diagonal, L=lower, and U=Upper parts of matrix A

• A= D − L − U

Parallelism in Gauss-Seidel Iteration

• Conventional Gauss-Seidel canNOT be easily parallelized because it has sequential update order

x_{n}^{(k)} requires x_{n−1}^{(k)} …
x_{n−1}^{(k)} requires x_{n−2}^{(k)} …
⋮
x_{2}^{(k)} requires x_{1}^{(k)}, etc.

• There are two parallelization strategies

1. Keep same algorithm—pipelined wavefront method

2. Change algorithm—red/black two-phase

6 × 6 Matrix Example Stencil

• We can exploit special structure of matrix of A (e.g., x_{16}^{(k)} requires x_{15}^{(k)}, x_{17}^{(k−1)}, x_{10}^{(k)}, x_{22}^{(k−1)}, etc.).

CPS615FEM, G.C.Fox, Oct 95  35  gcf@npac.syr.edu

6 × 6 Matrix—Wavefront Parallelism for Gauss Seidel

• Here number l labels l = 1 … 11 phases

• Each phase l ONLY depends on \underline{x}^{(k)} components from (l − 1) th phase

• Can execute each update in a given phase concurrently

• Some Load Imbalance

• Speedup of 2.5 in example of three processors

The Red-Black Two Phase Parallel Gauss Seidel Iteration

Red ⊙ Black •

• Note calculation of ⊙ only depends on • points and vice versa

• So two phase algorithm. Calculate all ⊙—Phase I. Then all •—Phase II.

Analysis of Parallel Red Black Gauss Seidel

• This can be parallelized as easily as Jacobi! • Phase I follows by Phase II.

• Red black is a different ordering from original Gauss Seidel.

• For (some) elliptic problems, can show red black as good as and maybe better than original ordering.

— i.e., red black does not need more iterations than pure Gauss Seidel.

• For related hyperbolic problems—see NAS benchmark SSOR discussion in CPS713—not so clear.

• For higher order differencing, need more phases (more colors) but principle the same.

• Also not clear how to use if irregular matrix, as in sparse matrices from power systems.

Eigenvalues of Gauss Seidel Iteration Matrix

• It was relatively easy to find the eigenvalues and eigenfunctions of Jacobi’s method for a simple geometry.

• Gauss-Seidel it is not so easy as matrix M^{−1}N is not symmetric. There is a nice discussion in

“Numerical Computation of Internal and External Flows,” C. Hirsch, Wiley, 1988.

• Eigenvalues λ_{GS} of Gauss-Seidel Iteration Matrix are related to those of Jacobi Iteration λ_{J} by λ_{GS} = λ_{J}^{2}

• We showed that

1 − ϵ < ∼ λ_{J} < ∼ −1 + ϵ ϵ ∼ \frac{π^{2}}{2N^{2}}

• Therefore,

|λ_{GS}|_{max} = (|λ_{J}|_{max})^{2} ∼ 1 − 2ϵ

Comparison of Convergence of Gauss-Seidel and Jacobi Iteration

• Convergence can be characterized by a given fractional error (say .001).

• This correponds to

|λ|_{max}^{k} ∼ .001

where k number of required iterations

• exp [−k (1 − |λ|_{max})] = .001 i.e., value of k ∝ \frac{1}{(1−|λ|_{max})}

• We showed Gauss Seidel doubles value of

1/(1 − |λ|_{max}),

i.e., it converges in half as many iterations as Jacobi

Successive Overrelaxation Iteration Method (SOR)

• This is relaxation (acceleration) method applied to Gauss-Seidel.

• In Jacobi, the optimal case was ω = 1 (i.e., no gain in using relaxation).

• However, for Gauss-Seidel, we will find ω ∼ 2 is best and does improve the method.

• The spectral radius (|λ|_{max}) of Jacobi is given by

ρ_{J} = 1 − \frac{π^{2}}{2N^{2}}

and for conventional (unaccelerated) Gauss Seidel by

ρ_{GS} = 1 − \frac{π^{2}}{N^{2}}

• For choice

ω_{opt} ∼ 2 (1 − \frac{π}{N} )

one can show that

ρ_{SOR} ∼ 1 − \frac{π}{N}

Convergence of SOR Compared to Jacobi and Gauss Seidel

• Remember number of iterations is proportional to

\frac{1}{1 − ρ} = \frac{2N^{2}}{π^{2}} Jacobi
= \frac{N^{2}}{π^{2}} Gauss Seidel
= \frac{N}{π}  SOR

— i.e., number of iterations decreases in SOR by a factor N/π, which is large!

• Unfortunately, the above analysis is only valid in simple cases

• simple geometry

• some orderings in Gauss-Seidel

Estimate of Over Relaxation Parameter ω

• Thus, in general for practical problems, one cannot precalculate optimal overrelaxation parameter.

• Find ω by numerical experimentation.

• However, note that one can expect a choice 1 < ω < 2 to be better than natural ω = 1.

• One can use the estimate

ω_{opt} = 2 − \frac{constant}{N}
= 2 − constant · h

(h is grid spacing)

Pseudo Code for SOR—Successive Over Relaxation

Choose an initial guess x^{(0)} to the solution x. for k = 1,2,…

for i = 1,2,…,n σ = 0 for j = 1,2,…,i − 1

σ = σ + a_{i,j} x_{j}^{(k)} end for j = i + 1,…,n

σ = σ + a_{i,j} x_{j}^{(k−1)} end σ = (b_{i} − σ) /a_{i,i} x_{i}^{(k)} = x_{i}^{(k−1)} + ω (σ − x_{i}^{(k−1)})

end check convergence; continue if necessary

end

The SOR Method