Example 1: Tutorial

This example should serve as a tutorial of the ALFA framework. We analyze several components of a two-grid method to solve the linear system of equations $Ax=b$, where $A$ corresponds to the two-dimensional Laplacian discretized on an equidistant rectangular lattice.

This example is structured as follows.

  1. We define the underyling linear operator: The second order approximation of the two-dimensional Laplacian $L$ obtained via finite central differences (5-point stencil). We further compute its spectrum.
  2. We show the connection of the operator $L$ to the system matrix $A$.
  3. We introduce and analyze the Jacobi-method.
  4. We define the coarse grid correction.
  5. We analyze the two-grid method using the Jacobi-method as a smoother.
  6. We introduce the lexicographic Gauss-Seidel and red-black Gauss-Seidel smoother and analyze the corresponding two-grid method.
  7. We use this framework to prototype an actual two-grid method.

Importing required packages

using ALFA
using LinearAlgebra
using Plots

The discretized Laplacian $L$ in 2D

We are going to define the 2D discretized Laplacian $L:\mathcal{L}(\mathbb{L}^s(\mathcal{A})) \rightarrow \mathcal{L}(\mathbb{L}^s(\mathcal{A}))$ on an equidistant rectangular lattice. Thus, $\mathbb{L}^s(\mathcal{A})$ describes the underlying structure of the domain and codomain of the operator $L$.

Definition of the underlying lattice structure

First, we need to define a basis of the underlying translational-invariance: an equidistant rectangular lattice is given by

\[\mathcal{A} = \frac{1}{h}\left[\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}\right] = \left[\begin{matrix} \mathcal{a}_1 & \mathcal{a}_2\end{matrix}\right].\]

The lattice generated by $\mathcal{A}$ is the set

\[\mathbb{L}(\mathcal{A}) = \mathcal{A}\mathbb{Z}^2 = \{x = \alpha_1\mathcal{a}_1 + \alpha_2\mathcal{a}_2 \ : \ \alpha_1,\alpha_2 \in \mathbb{Z}\}\]

The class ALFA.Lattice corresponds to $\mathbb{L}(\mathcal{A})$ and is initialized with the matrix $\mathcal{A}$.

h = .1  # define h
A = h*[1 0; 0 1]  # 2x2 identity matrix scaled with h

# initialize the lattice;
#In the constructor we need to specify the size/dimensionality (N=2)
# and the datatype of the lattice basis (T=Float64)
#  (the alternative to Float is Rational{BigInt})
LA = ALFA.Lattice{2, Float64}(A)
ALFA.Lattice{2, Float64}([0.1 0.0; 0.0 0.1])

We can plot a section of the lattice. Black arrows correspond to the primitive vectors $a_1$ and $a_2$

plot(LA)

Definition of the domain and codomain of the operator $L$

The structure element $s$ of $\mathbb{L}^s(\mathcal{A})$ corresponds of the location of the unknowns. We can simply choose $s=(s_1)$, where $s_1=(0,0)$, such that the crystal points coincide with the lattice structure.

The struct ALFA.Crystal is used to represent both crystals corresponding to the domain and codomain of an operator. We initialize it with ALFA.Crystal{N,T}(LA,s_domain, s_codomain), where sdomain corresponds to the structure element of the domain, and scodomain corresponds to the structure element of the codomain. (In our case we have sdomain = scodomain = $s$):

Domain = [[0,0]]
Codomain = [[0,0]]
CA = ALFA.Crystal{2,Float64}(LA,Domain,Codomain)
Lattice Basis: ALFA.Lattice{2, Float64}([0.1 0.0; 0.0 0.1])
Domain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]
Codomain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]

We can have a plot function to plot a section of a crystal:

plot(CA)

Initializing the operator $L$

As we have defined the underlying domain and codomain of our operator, we can initialize the multiplication operator $L$.

This class represents a multiplication operator corresponding to

\[(Lf)(x) = \sum_{y \in \mathbb{Z}^\text{d} } m_L^{y} \cdot f(x+\mathcal{A}y),\]

for all $x \in \mathbb{L}(\mathcal{A})$, where $d=2$ is the dimensionality. Note, that the position of the operator is given in fractional coordinates, not in cartesian coordinates $\mathcal{A}y$, such that $y$ is always integral.

L = ALFA.CrystalOperator{2,Float64}(CA)
Lattice Basis: ALFA.Lattice{2, Float64}([0.1 0.0; 0.0 0.1])
Domain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]
Codomain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]
Multiplier: ALFA.Multiplier[]

Adding the multipliers of $L$

We use push!(L,m) to add multipliers $m_L^{y}$ to the operator in order to define the discretized Laplacian. As the structure elements of the domain and codomain are both $1$-dimensional, the matrices m_L^{y} are of size $1\times 1$.

The multipliers are then saved within a SortedSet L.m which is lexicographically ordered with respect to the position y.

push!(L, ALFA.Multiplier([0 0], [-4/h^2]))
push!(L, ALFA.Multiplier([0 -1], [1/h^2]))
push!(L, ALFA.Multiplier([0 1], [1/h^2]))
push!(L, ALFA.Multiplier([1 0], [1/h^2]))
push!(L, ALFA.Multiplier([-1 0], [1/h^2]))

@show L
Lattice Basis: ALFA.Lattice{2, Float64}([0.1 0.0; 0.0 0.1])
Domain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]
Codomain: 1-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [0.0, 0.0]
Multiplier: 5-element Vector{ALFA.Multiplier}:
 ALFA.Multiplier{2}([-1, 0], [99.99999999999999;;])
 ALFA.Multiplier{2}([0, -1], [99.99999999999999;;])
 ALFA.Multiplier{2}([0, 0], [-399.99999999999994;;])
 ALFA.Multiplier{2}([0, 1], [99.99999999999999;;])
 ALFA.Multiplier{2}([1, 0], [99.99999999999999;;])

A schematic representation can be created with the plot command.

plot(L)

Computing the spectrum of $L$

The spectrum of a CrystalOperator is computed via the $ALFA.symbol(L::CrystalOperator,k::Vector)$ method, where $k$ is the frequency/wavevector.

ALFA.symbol(L,[0, 0])
1×1 Matrix{ComplexF64}:
 0.0 + 0.0im

In order to compute or approximate the "complete" spectrum, we divide $[0,1)^\text{d}$ into $N^d$ equidistant points to get a discretization of the primitive cell of the dual lattice $\mathcal{A}^{-T}[0,1)^\text{d}$. Next, we compute the eigenvalues for each symbol(L,k), where k is sampled on these $N^d$ points.

ALFA.eigvals(L,N=10)
100-element Vector{ComplexF64}:
                 0.0 + 0.0im
  -38.19660112501053 + 0.0im
  -38.19660112501053 + 0.0im
  -38.19660112501053 + 0.0im
  -38.19660112501053 + 0.0im
  -76.39320225002103 + 0.0im
  -76.39320225002103 + 0.0im
  -76.39320225002103 + 0.0im
  -76.39320225002103 + 0.0im
 -138.19660112501043 + 0.0im
                     ⋮
  -723.6067977499788 + 0.0im
  -723.6067977499788 + 0.0im
  -723.6067977499788 + 0.0im
  -723.6067977499788 + 0.0im
  -761.8033988749894 - 1.964386723728473e-15im
  -761.8033988749894 - 1.964386723728473e-15im
  -761.8033988749894 + 0.0im
  -761.8033988749894 + 0.0im
  -799.9999999999999 + 0.0im

We may also save everything in a dataframe via ALFA.eigvals_df(L,N=N) or directly produce a plot of the absolut part of the spectrum as follows

plotspectrum(L,N=20)

Obtaining a system matrix $A$ from a multiplication operator $L$

We can obtain the system matrix $A$ corresponding to the a discretization of the Laplacian on the unit square [0,1]^2 with periodic boundary conditions in two steps.

Step 1

We rewrite $L$ with respect to the sublattice $\mathbb{L}(\mathcal{Z}) \subset \mathbb{L}(\mathcal{A})$ with $\mathcal{Z}=\left(\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}\right)$.

Due to the fact that we find $1/h^2 = 100$ lattice points in $Z[0,1)^2 \cap \mathbb{L}(\mathcal{A}) =: \tilde{s}$, we obtain an operator

\[\tilde{A}:\mathcal{L}(\mathbb{L}^{\tilde{s}}(\mathcal{Z})) \rightarrow \mathcal{L}(\mathbb{L}^{\tilde{s}}(\mathcal{A})) \text{ with } (\tilde{A} g)(x) = \sum_{y \in \mathbb{Z}^\text{d} } m_{\tilde{A}}^{y} \cdot g(x+\mathcal{Z}y), \quad m_{\tilde{A}} \in \mathbb{C}^{100\times 100}\]

using SparseArrays
tA = ALFA.wrtLattice(L,ALFA.Lattice{2, Float64}([1 0; 0 1]))
tA = ALFA.normalize(tA); # make sure all points of the structure elements are within [0,1)^2 and lexicographically ordered.

We create a spy plot of all five multipliers of tA. One multiplier corresponds to the interaction of the unknowns within the unit cell. The other four multipliers correspond to the connections beyond the boundaries:

using SparseArrays
parr = []
for x in tA.M
    p = spy(SparseMatrixCSC(x.mat), c=:blues, title="position: " * string(x.pos))
    push!(parr, p)
end
plot(parr...);

Step two

As we impose periodic boundary conditions on the unit square, all lattice points of $\mathbb{L}(Z)$ are identified. Thus, we obtain the system matrix $A$ by adding all multipliers of tA.

tAm = sum(x.mat for x in tA.M);
spy(SparseMatrixCSC(tAm), c=:blues);

We can compute the eigenvalues and eigenvectors of $A$ and compare them with the eigenvalues computed via LFA:

eig1 = abs.(sort(eigvals(tAm), by=abs))
eig2 = abs.(ALFA.eigvals(L,N=10, by=abs))

p = plot(eig1, label="eigenvalues of the system matrix")
plot!(eig2, label="eigenvalues computed via the symbol")
plot!(abs.(eig1.-eig2), label="pairwise difference of the eigenvalues.")

The eigenvalues are equal as long as N=1/h, where $h$ corresponds to the lattice spacing and $N$ corresponds to the discretization of the frequency space.

Analysis of stationary iterative methods

Given a linear system of equations (LSE) $Ax=b$, the Jacobi method produces iterates via

\[x_{k+1} \leftarrow (I-S^{-1}A)x_k + S^{-1}b,\]

where $I$ is the identity and $S$ is simply the diagonal of $A$.

Denote with $x^* = A^{-1}b$ the solution of the LSE and with $e_k = x^* - x_k$ the $k$th error. Since

\[e_{k+1} = (I-S^{-1}A)e_k\]

the operator $G:=(I-S^{-1}A)$ is called error propagator. We are interested in the spectral radius $\rho(G)=\max\{|\lambda | \ : \ \lambda \text{ Eigenvalue of }G \}$ of the error propagator, as we have

\[||e_{k+1}|| \approx \rho(G) ||e_{k}||\]

for large k (if the initial error has a component in the direction of the eigenvector corresponding to the largest absolute eigenvalue).

The Jacobi method

Definition of the Jacobi-method

We compute the spectrum via LFA, thus we need to define $I$ and $S$ as multiplication operators. In case of the Jacobi method, the multipliers of

\[S:\mathcal{L}(\mathbb{L}^s(\mathcal{A})) \rightarrow \mathcal{L}(\mathbb{L}^s(\mathcal{A})), \quad\quad \quad (Sf)(x) = \sum_{y \in \mathbb{Z}^\text{d} } m_S^{y} \cdot f(x+\mathcal{A}y),\]

are given by

\[m_S^{y} = \begin{cases} m_L^{y} & \text{ if } y = 0 \\ 0 & \text{ else.} \end{cases}\]

S_jac = ALFA.CrystalOperatorCopyWithMultipliers(L)

plot(S_jac)

Analysis of the Jacobi-method

We analyze the error propagator of the Jacobi method using underrelaxation of .8.

f_jac =:(I-0.8*inv($S_jac)*$L) # construct an expression holding the CrystalOperators S and L.
oc = ALFA.OperatorComposition(f_jac)

plotspectrum(oc, 40) # plot the absolute part of the spectrum using 40^2 points.