9  Relaxation Methods

Lecture 9

Overview

1. Convergence of Jacobi & Gauss-Seidel
2. Sucessive over-relaxed (SOR) method

Note

I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.

Warning

These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Tip

This lecture is mostly based on Burden, Faires, and Burden (2015) sections 7.4.

using Plots, LaTeXStrings, LinearAlgebra

9.1 Recap

Burden, Faires, and Burden (2015) 7.3

We will let \|\cdot\| to be the induced matrix norm (also called natural or an operator norm) from some vector norm |\cdot|: that is,

\begin{align} \|A\| := \max_{|x| = 1} |Ax| % = \max_{ x \not= 0} \frac{|Ax|}{|x|}. \end{align}

Notice:

• |Ax| \leq \|A\| |x|,
  
• \|AB\| = \max_{|x| = 1} |ABx| \leq \|A\| \max_{|x| = 1} |Bx| = \|A\| \|B\|,
  
• In particular, \|A^k\| \leq \|A\|^k,
  
• Recall \rho(A) := \max_{\lambda\in\sigma(A)}|\lambda| is the spectral radius,
  
• Recall that for every A and \varepsilon>0, there exists a induced matrix norm \| \cdot \|_\varepsilon such that

\begin{align} \rho(A) \leq \|A\|_{\varepsilon} \leq \rho(A) + \varepsilon \end{align}

• The left-hand side of this inequality holds for every submultiplicative matrix norm.

Therefore, if \rho(T) < 1, then there exists a an induced matrix norm \| \cdot \| such that \|T\| < 1. As a result, we have | T^{k+1} x | \leq \|T\|^{k+1} |x| \to 0 as k \to \infty.

Exercise 9.1 Let p be a polynomial. Show that \sigma( p(A) ) = \{ p(\lambda) : \lambda \in \sigma(A) \}.

In particular, we have \rho(A^k) = \rho(A)^k.

Last time, we saw Jacobi and Gauss Seidel methods:

function Jacobi(A,b,x; Niter=100, tol=1e-10) 
    
    X = [x]

    L = -tril(A, -1)
    U = -triu(A, 1)
    D = Diagonal(diag(A))

    if (any(iszero, diag(D)))
        @warn "diagonal of A is not invertible"
        return X
    end

    for i ∈ 1:Niter
        push!( X, inv(D)*( b + (L+U)*X[end] ) ) 
        r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)
        if (r < tol)
            return X
        end 
    end

    @warn "max iterations exceeded"
    return X
end;

function GaussSeidel(A,b,x; Niter=100, tol=1e-10) 
    
    X = [x]

    L = -tril(A, -1)
    U = -triu(A, 1)
    D = Diagonal(diag(A))

    if (any(iszero, diag(D)))
        @warn "diagonal of A is not invertible"
        return X
    end

    for i ∈ 1:Niter
        push!( X, inv(D-L)*( b + U*X[end] ) )
        r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)
        if (r < tol)
            return X
        end 
    end

    @warn "max iterations exceeded"
    return X
end;

Example 9.1 Applying the Jacobi and Gauss-Seidel iterations to Ax = b with

b = [6, 25, -11, 15];
A = [10 -1 2 0; 
    -1 11 -1 3; 
    2 -1 10 -1; 
    0 3 -1 8]

4×4 Matrix{Int64}:
 10  -1   2   0
 -1  11  -1   3
  2  -1  10  -1
  0   3  -1   8

gives:

X, Y = Jacobi(A,b,[0. ,0., 0., 0.]), GaussSeidel(A,b,[0. ,0., 0., 0.])
display( X[end] )
display( Y[end] )

4-element Vector{Float64}:
  1.0000000000277243
  1.9999999999534517
 -0.9999999999639979
  0.9999999999487015

4-element Vector{Float64}:
  0.9999999999812351
  1.9999999999852736
 -0.9999999999921585
  1.0000000000065026

Exercise 9.2 Verify that

\begin{align} x = \begin{pmatrix} 1 \\ 2 \\ -1 \\ 1 \end{pmatrix} \end{align}

solves Ax = b with A and b given by Example 10.1.

In fact, we have the following error curves:

Exercise 9.3 What is the rate of approximation in the error curves of Example 10.1?

9.2 Convergence

Both Jacobi and Gauss-Siedel can be written as

\begin{align} x^{(k+1)} = T x^{(k)} + c \tag{T} \end{align}

Exercise 9.4 Write down what T is for (i) Jacobi, (ii) Gauss-Siedel.

Theorem 9.1 The iteration \text{(T)} converges to the unique solution of x = Tx + c if and only if \rho(T) < 1.

Proof.

We will use the following geometric series result

Lemma 9.1 If \rho(T) < 1, then I-T is invertible with

\begin{align} (I - T)^{-1} = \sum_{n=0}^\infty T^n. \end{align}

Since \rho(T) < 1, there exists an induced matrix norm \| \cdot\| such that \|T\|<1 and thus \|T^{k+1}\| \leq \|T\|^{k+1} \to 0 as k \to \infty. That is, |T^{k+1}x| \to 0 as k\to\infty for all x. Using this and

\begin{align} (I-T)\big[ I + T + T^2 + \cdots + T^k] = I - T^{k+1},\nonumber \end{align}

we find that \sum_{n=0}^\infty T^n x = (I-T)^{-1}x, as required.

Let’s return to the proof of the theorem. If \rho(T) < 1, then

\begin{align} x^{(k+1)} &= Tx^{(k)} + c \nonumber\\ % &= T\big[ Tx^{(k-1)} + c \big] + c \nonumber\\ % &= T^3 x^{(k-2)} + \big[ I + T + T^2 \big] + c \nonumber\\ % &= \cdots = \nonumber\\ % &= T^{k+1} x^{(0)} + \sum_{n=0}^{k} T^{n} c \end{align}

The lemma implies that the limit of this as k\to\infty is given by

\begin{align} x := 0 + \sum_{n=0}^\infty T^n c = (I - T)^{-1} c. \end{align}

That is, x = Tx + c.

Suppose that \text{(T)} converges to x = Tx + c for any choice of x^{(0)}. Then,

\begin{align} x - x^{(k+1)} &= T(x - x^{(k)}) \nonumber\\ % &= \cdots = \nonumber\\ % &= T^{(k+1)}(x - x^{(0)}). \end{align}

Because we can choose any x^{(0)}, we have \|T^{k+1}\| \to 0 as k\to \infty. In particular, we have

\begin{align} \rho(T)^{k+1} &= \rho(T^{k+1}) \nonumber\\ % &\leq \|T^{k+1}\| \to 0 \end{align}

as k \to \infty. That is, \rho(T) < 1.

Exercise 9.5 Let

\begin{align} A = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}, % \,\, b = \begin{pmatrix} 1 \\ 0 \end{pmatrix}.\nonumber \end{align}

Do you expect Jacobi and Gauss-Seidel methods to converge? What rate of convergence?

A = [0 0; 1 0]; b = [1, 0]

X, Y = Jacobi(A,b,[0. ,0.]), GaussSeidel(A,b,[0.,0.])
display(X)

┌ Warning: diagonal of A is not invertible

└ @ Main In[3]:10

┌ Warning: diagonal of A is not invertible

└ @ Main In[3]:35

1-element Vector{Vector{Float64}}:
 [0.0, 0.0]

Exercise 9.6 (Exercise 7.3 #9 from Burden, Faires, and Burden (2015)) Let

\begin{align} A = \begin{pmatrix} 2 & -1 & 1 \\ 2 & 2 & 2 \\ -1 & -1 & 2 \end{pmatrix}, % \qquad b = \begin{pmatrix} -1 \\ 4\\ -5 \end{pmatrix}. \end{align}

Do you expect Jacobi and Gauss-Seidel methods to converge? What rate of convergence?

Here are the eigenvalues of D^{-1}(L + U) and (D-L)^{-1} U (using the same notations as before for D,L,U):

A = [2 -1 1 ; 2 2 2 ; -1 -1 2]; b = [-1, 4, -5]

L = -tril(A, -1)
U = -triu(A, 1)
D = Diagonal(diag(A))

display( eigen(inv(D)*(L+U)).values )
display( eigen(inv(D-L)*U).values )

3-element Vector{ComplexF64}:
 2.6009723803125757e-17 + 0.0im
  5.551115123125783e-17 - 1.1180339887498956im
  5.551115123125783e-17 + 1.1180339887498956im

3-element Vector{Float64}:
 -0.5
 -0.5
  0.0

Here, we apply Jacobi and Gauss-Seidel iterations:

A = [2 -1 1 ; 2 2 2 ; -1 -1 2]; b = [-1, 4, -5]
x0 = [0.,0.,0.]

X, Y = Jacobi(A,b,x0), GaussSeidel(A,b,x0)
display(X[end])
display(Y[end])

┌ Warning: max iterations exceeded

└ @ Main In[3]:22

3-element Vector{Float64}:
  -42037.95392974449
 -168153.81571897797
   42037.95392974449

3-element Vector{Float64}:
  1.0000000000636646
  1.9999999999326974
 -1.000000000001819

9.3 Relaxation Schemes

Mostly based on Burden, Faires, and Burden (2015) section 7.4

function SOR(A,b,x; ω=1/2, Niter=1000, tol=1e-8) 
    
    X = [x]

    L = -tril(A, -1)
    U = -triu(A, 1)
    D = Diagonal(diag(A))

    if (any(iszero, diag(D)))
        @warn "diagonal of A is not invertible"
        return X
    end

    for i ∈ 1:Niter
        push!( X, inv(D-ω*L)*( ω*b + ((1-ω)D+ω*U)*X[end] ) )
        r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)
        if (r < tol)
            return X
        end 
    end

    @warn "max iterations exceeded"
    return X
end;

Example 9.2 Here, we apply SOR with a few different choices of \omega:

A = [4 3 0; 3 4 -1; 0 -1 4]; b = [24,30,-24]
x0 = [0.,0.,0.]

X = Jacobi(A,b,x0)
Y = GaussSeidel(A,b,x0)
Z = []
Ω = 1.125:.25:2;
for ω ∈ Ω
    push!( Z, SOR(A,b,x0; ω=ω) )
end

┌ Warning: max iterations exceeded

└ @ Main In[3]:22

TipTO-DO

• Assignment 3: Due Wednesday 11:59 PM,
• (Sign-up to office hours if you would like)

Chapter 1 reading:

• Burden, Faires, and Burden (2015) sections 5.1 - 5.6,
  
• Driscoll and Braun (2018) zero stability

Chapter 2 reading:

• Recap matrices: sections 7.1 - 7.2 of Burden, Faires, and Burden (2015)
• Jacobi & Gauss-Seidel: section 7.3 Burden, Faires, and Burden (2015)
  
• Today: section 7.4 Burden, Faires, and Burden (2015)

Next:

• Read: section 7.6 Burden, Faires, and Burden (2015)

Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.

Driscoll, Tobin A, and Richard J Braun. 2018. Fundamentals of Numerical Computation. https://fncbook.com/.