# Operators and Matrices

In this tutorial we explain how the hyperbolic graph structure generated by hypertiling can be used to construct lattice operators. You will learn how the standard (listed) neighbour format can be easily transferred into adjacency and degree matrices, convientinly encoded in scipy’s sparse matrix format. This allows to construct discrete (finite-difference) representation of several second order differential operators, such as a Laplacian or Helmholtz operator, which even support variable (position dependent) weights and Dirichlet boundaries.

:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt


## Basic Usage

Import tiling and graph conversion functions

:

from hypertiling import HyperbolicTiling
from hypertiling.operators import adjacency_matrix_sparse, degree_matrix_sparse, identity_matrix_sparse

:

T = HyperbolicTiling(7,3,3)


Note that almost everything in this notebook can also be done using the HyperbolicGraph class instead of HyperbolicTiling

Extract neighbours using the kernel-default method and convert to adjacency matrix

:

nbrs = T.get_nbrs_list()


The sparse representation can be put back to a full matrix with the todense() function. We use it for plotting. As can be seen, the adjaceny matrix is indeed very sparse, and furthermore symmetric (due to the graph being undirected). Every connection is represented by yellow color.

:

plt.matshow(A.todense()); plt.colorbar(); Similarly, we can construct a degree matrix, which contains the number of connections attached to each vertex. As we can see for central notes this is 7 (the parameter p of the tiling), whereas for boundary points it can be fewer than p.

:

G = degree_matrix_sparse(nbrs)

plt.matshow(G.todense()); plt.colorbar(); We are now prepared to build a Laplacian for our lattice …

:

A = adjacency_matrix_sparse(nbrs)
G = degree_matrix_sparse(nbrs)

D = A - G

plt.matshow(D.todense()); plt.colorbar(); … as well as a Helmholtz operator containting a mass term. We also added weights for the Laplacian part. In case we want to use constant weights, it is enought to multiply them as a factor. For vertex-dependend weights the weight keyword provided by adjacency_matrix_sparse and degree_matrix_sparse can be used.

:

w = 2 # weight
m = 5 # mass

:

A = adjacency_matrix_sparse(nbrs)
G = degree_matrix_sparse(nbrs)
I = identity_matrix_sparse(nbrs)  # identical: sp.sparse.identity(len(T))

H = w*(A - G) - m*I

plt.matshow(D.todense()); plt.colorbar(); ## Dirichlet boundary

:

from hypertiling.graphics.svg import make_svg, draw_svg


The operator functions provide the boundary keyword which takes a boolean array, encoding whether a vertex is considered a boundary vertex. For boundary points, the corresponding rows in the matrix are left empty, it therefore becomes non-symmetric. This way it is ensured that the boundary point act on their interiour neighbours, but stay themselves constant (Dirichlet boundary)

We use a helper function to determine at which maximal radial distance the lattice can be divded into a closed boundary and interiour bulk nodes

:

from hypertiling.experimental.boundary import maximal_radial_cutoff

/home/schrauth/code/hypertiling/hypertiling/experimental/__init__.py:3: UserWarning: Experimental functions. Use with caution!
warnings.warn("Experimental functions. Use with caution!")

:

T = HyperbolicTiling(5,4,3)               # generate tiling
boundary, rc = maximal_radial_cutoff(T)   # locate boundary


This is an example; boundary cells are plotted in black, bulk cells in white color

:

svg = make_svg(T, boundary.astype("int"), cmap="Greys", edgecolor="grey", lw=0.1)
draw_svg(svg) ## Helmholtz equation

Now we are ready to solve a Helmholtz equation in the bulk region of a hyperbolic lattice

:

T = HyperbolicTiling(7,3,6)

# locate boundary

# extract neighbours
nbrs = T.get_nbrs_list()

# cut away boundary in one quadrant; those cells belong to the bulk region now
for i in range(len(T)):
boundary[i] = False

bulk = np.invert(boundary)

# mass parameter
m = 1

# weights
weights = [list(np.random.rand(len(nb))+5) for nb in nbrs] # some random weights between [4,5]

:

# construct discrete differential operator
G = degree_matrix_sparse(nbrs, weights, boundary=boundary)     # Degree matrix
D = A-G                                                        # Laplacian
M = m*identity_matrix_sparse(nbrs, boundary=boundary)          # Mass matrix
B = sp.sparse.diags(boundary.astype("int"))                    # Boundary matrix

H = D - M + B  # Helmholtzian


Reminder: The boundary is included as trivial lines in the operator matrix; in some literature this approach is refered to as “ghost cells”; they are not updated but affect their bulk neighbours

The right hand side of the resulting linear system thus carries precisely the Dirichlet boundary values for those sites

Now are ready to solve the linear system $$Hx=y$$, where $$H$$ is the differential operator, $$x$$ is the solution vector and $$y$$ is the right hand side (rhs). We use the GMRES solver (Generalized minimal residual method).

:

from scipy.sparse.linalg import gmres

:

rhs = boundary.astype("int") # right hand side
x0 = np.ones(len(T))         # initial guess

sol, _ = gmres(H, rhs, x0=x0) # solve


Boundary cells at angles between 90 and 360 degrees have been fixed to a value of 1 (dark blue), whereas all remaining cells (including the boundary between 0 and 90 degress) have been left free. The colormap is roughly symmetric, hence white cells have a value of 0 and dark red cells a value of -1.

:

svg = make_svg(T, sol, cmap="RdBu", edgecolor="white", lw=0.1)
draw_svg(svg) ## Laplace equation

Now we model an electrostatic problem, where the boundary nodes are held to a fixed electric potential

:

T = HyperbolicTiling(3,7,6, center="vertex")

# locate boundary

# extract neighbours
nbrs = T.get_nbrs_list(method="RO")

weights = [list(np.random.rand(len(nb))+5) for nb in nbrs] # some random weights between [4,5]

:

# construct discrete differential operator
# no weights -> homogeneous surface
G = degree_matrix_sparse(nbrs, weights=weights, boundary=boundary)      # Degree matrix
D = A-G                                                # Laplacian
B = sp.sparse.diags(boundary.astype("int"))            # Boundary matrix

H = D + B  # Laplacian


Build Dirichlet boundary and solve Laplace equation in the bulk

:

# some regions are fixed to a value of 1, some to a value of -1
rhs = boundary.astype("int")
for i in range(len(T)):
if boundary[i]:
rhs[i] = -1
rhs[i] = -1

x0 = np.ones(len(T)) # initial guess

sol, _ = gmres(H, rhs, x0=x0, tol=1e-2)  # solve


Plot the result. This time we fixed the entire boundary to values of either -1 (red) or +1 (blue). In the interiour of the lattice, the Laplace equation provides a nice, smooth solution (electrostatic potential).

:

svg = make_svg(T, sol, cmap="RdBu", edgecolor="white", lw=0.1)
draw_svg(svg) 