Harmonics of the Sphere

This page shows how to use Decapodes tooling to explore the harmonics of a discrete manifold. This isn't using any Decapodes specific code, but it is emblematic of a more advanced analysis you might want to do on your Decapode.

In this case we are trying to visualize the roots of the Laplacian on a discrete manifold.

Load the dependencies

# Meshing:
using CombinatorialSpaces
using CoordRefSystems
using GeometryBasics: Point3
const Point3D = Point3{Float64};

# Visualization:
using CairoMakie

# Simulation:
using LinearAlgebra

Load the mesh

const RADIUS = 1.0
s = loadmesh(Icosphere(3, RADIUS));
sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s);
subdivide_duals!(sd, Barycenter());

Compute the Laplacian eigenvectors using LinearAlgebra.eigen. This requires making the sparse Laplacian matrix dense with collect. Alternatively, use Arpack.jl.

Δ0 = -Δ(0,sd)
λ = eigen(collect(Δ0))

LinearAlgebra.Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
162-element Vector{Float64}:
  9.977942897419346e-16
  1.850195339821431
  1.9229341081041116
  1.9658810495474286
  5.146220910776434
  5.528877811211735
  5.639164681656218
  5.737050666884109
  5.85994369579021
  9.3172725539988
  ⋮
 77.17507216665005
 86.50765980625566
 86.51040566694232
 86.58168500302766
 86.59363600121723
 86.67718662766426
 86.68365853045188
 86.68546749816456
 86.68738445511194
vectors:
162×162 Matrix{Float64}:
 0.0785674  -0.00139719  -0.104545   …  -0.000848792  -0.00294379
 0.0785674   0.0696797    0.0192111     -0.153257     -0.00715503
 0.0785674   0.107855    -0.0713353      0.488914      0.487552
 0.0785674   0.00417104  -0.128525      -0.00110173   -0.00296474
 0.0785674  -0.112362    -0.0657565     -0.575396      0.387291
 0.0785674  -0.0834091    0.0226547  …   0.117539      0.00955895
 0.0785674   0.114972     0.0665464     -0.216674     -0.476139
 0.0785674   0.0820901   -0.0259179     -0.115597     -0.00385017
 0.0785674  -0.0704939   -0.0222209      0.155239      0.0135661
 0.0785674  -0.110218     0.0716917      0.301402     -0.409684
 ⋮                                   ⋱   ⋮            
 0.0785674   0.0323531   -0.130654      -0.00544982   -0.00626171
 0.0785674   0.0306632   -0.122617      -0.00539986   -0.00625665
 0.0785674   0.0654118   -0.113619   …   0.0338602     0.0338033
 0.0785674  -0.0284136   -0.0111406      0.0102456     0.000514537
 0.0785674   0.02021     -0.0121854     -0.0124337    -0.000887686
 0.0785674  -0.0018298   -0.0418315      0.000199872  -0.000553123
 0.0785674   0.0856251   -0.068554       0.0357922     0.0341249
 0.0785674   0.0508548   -0.0772451  …  -0.00421754   -0.00604669
 0.0785674   0.0734731   -0.040202      -0.0163659    -0.00595885

Let's check that our eigenvalues satisfy the right equation. The first eigenvector should be the kernel of the laplacian. So the following norm should be close to 0.

q1 = λ.vectors[:,1]
norm(Δ0 *q1)

5.5079433242025446e-14

The first eigenvector is boring to visualize, because it is constant. So we will make some plots of the second eigenvector. If you run this on the desktop, you can use GLMakie and get an interactive plot to explore. We will just draw two angles.

q = λ.vectors[:,2]
fig = Figure()
Label(fig[1, 1, Top()], "Default Angle", padding = (0, 0, 5, 0))
ax = LScene(fig[1,1], scenekw=(lights=[],))
msh = CairoMakie.mesh!(ax, s, color=q)
Colorbar(fig[1,2], msh, size=32)
# Second Angle
Label(fig[2, 1, Top()], "Bottom Angle", padding = (0, 0, 5, 0))
ax = LScene(fig[2,1], scenekw=(lights=[],))
update_cam!(ax.scene, Vec3f(-1/2,-1/2,1.0/2), Vec3f(1,1,1), Vec3f(0, 0, 1))
msh = CairoMakie.mesh!(ax, s, color=q)
Colorbar(fig[2,2], msh, size=32)
fig

q = λ.vectors[:,12]
fig = Figure()
Label(fig[1, 1, Top()], "Default Angle", padding = (0, 0, 5, 0))
ax = LScene(fig[1,1], scenekw=(lights=[],))
msh = CairoMakie.mesh!(ax, s, color=q)
Colorbar(fig[1,2], msh, size=32)
# Second Angle
Label(fig[2, 1, Top()], "Bottom Angle", padding = (0, 0, 5, 0))
ax = LScene(fig[2,1], scenekw=(lights=[],))
update_cam!(ax.scene, Vec3f(-1/2,-1/2,1.0/2), Vec3f(1,1,1), Vec3f(0, 0, 1))
msh = CairoMakie.mesh!(ax, s, color=q)
Colorbar(fig[2,2], msh, size=32)
fig

Exploring solutions with Krylov methods

We can also use the information about the eigenvectors for spectral techniques in solving the equations. Krylov methods are a bridge between linear solvers and spectral information.

using Krylov
b = zeros(nv(sd))
b[1] = 1
b[end] = -1
x, stats = Krylov.gmres(Δ0, b, randn(nv(sd)), restart=true, memory=20, atol = 1e-10, rtol=1e-8, history=true, itmax=10000)
x̂ = x .- sum(x)./length(x)
norm(x̂)
stats
norm(Δ0*(x) - b)