Demonstrate the use of SymmetricEigensystem in examples (#881)

* Update eigen examples

* remove eigs

* Eigenfunctions in example

* Add comments

* remove checks in readmetests if there are no actual tests

* revert removing tests

* rename tunneling file

* fix filename

* make functions multiline

* separate function for interlace

docs/make.jl

@@ -33,7 +33,7 @@ end
 makedocs(
             format = Documenter.HTML(),
             sitename = "ApproxFun.jl",
-            authors = "Sheehan Olver, Jishnu Bhattacharya and contributors",
+            authors = "Sheehan Olver, Mikael Slevinsky, Jishnu Bhattacharya and contributors",
             pages = Any[
                     "Home" => "index.md",
                     "Usage" => Any[

examples/Eigenvalue.jl

@@ -16,26 +16,33 @@ p
 # If the solutions are not relatively constant near the boundary then one should push
 # the boundaries further out.
 
-include("Eigenvalue_tunnelling.jl")
-
-# We plot the first few eigenfunctions
-p = Plots.plot(V, legend=false)
-Plots.vline!([-Lx/2, Lx/2], linecolor=:black)
-p_twin = Plots.twinx(p)
-for k=1:6
-    Plots.plot!(p_twin, real(v[k]/norm(v[k]) + λ[k]), label="$k")
-end
-p
-
-# Note that the parity symmetry isn't preserved exactly at finite matrix sizes.
-# In general, it's better to preserve the symmetry of the operator matrices (see section below),
-# and projecting them on to the appropriate subspaces.
-
 # For problems with different contraints or boundary conditions,
 # `B` can be any zero functional constraint, e.g., `DefiniteIntegral()`.
 
-include("Eigenvalue_symmetric.jl")
+include("Eigenvalue_anharmonic.jl")
 
 # We plot the eigenvalues
 import Plots
-Plots.plot(λ; title = "Eigenvalues", legend=false)
+Plots.plot(λ, title = "Eigenvalues", legend=false)
+
+include("Eigenvalue_well_barrier.jl")
+
+# We plot the first few eigenfunctions offset by their eigenvalues.
+# The eigenfunctions appear in odd-even pairs as expected.
+
+import Plots
+using LinearAlgebra: norm
+p1 = Plots.plot(Vfull, legend=false, ylim=(-Inf, λ[14]), title="Eigenfunctions",
+    seriestype=:path, linestyle=:dash, linewidth=2)
+Plots.vline!([-Lx/2, Lx/2], color=:black)
+for k=1:12
+    Plots.plot!(real(v[k]/norm(v[k]) + λ[k]))
+end
+
+# Zoom into the ground state:
+p2 = Plots.plot(Vfull, legend=false, ylim=(-Inf, λ[3]), title="Ground state",
+    seriestype=:path, linestyle=:dash, linewidth=2)
+Plots.vline!([-Lx/2, Lx/2], color=:black)
+Plots.plot!(real(v[1]/norm(v[1]) + λ[1]), linewidth=2)
+
+Plots.plot(p1, p2, layout = 2)

examples/Eigenvalue_anharmonic.jl

@@ -0,0 +1,38 @@
+# ## Self-adjoint Eigenvalue Problem
+# Ref:
+# [J. L. Aurentz & R. M. Slevinsky (2019), arXiv:1903.08538](https://arxiv.org/abs/1903.08538)
+
+# We solve the confined anharmonic oscillator
+# ```math
+# \left[-\frac{d^2}{dx^2} + V(x)\right] u = λu,
+# ```
+# where ``u(\pm 10) = 0``, ``V(x) = ωx^2 + x^4``, and ``ω = 25``.
+
+using ApproxFun
+using LinearAlgebra
+using BandedMatrices
+
+# Define parameters
+ω = 25.0
+d = -10..10;
+S = Legendre(d); # Equivalently, Ultraspherical(0.5, d)
+
+# Construct the differential operator
+V = Fun(x -> ω * x^2 + x^4, S)
+L = -Derivative(S, 2) + V;
+
+# Boundary conditions that are used in the basis recombination
+B = Dirichlet(S);
+
+# The system may be recast as a generalized eigenvalue problem
+# ``A_S\,v = λ\, B_S\, v``, where ``A_S`` and ``B_S`` are symmetric band matrices.
+# We wrap the operators in `SymmetricEigensystem` to implicitly perform the basis recombination
+
+SEg = ApproxFun.SymmetricEigensystem(L, B);
+
+# We construct `n × n` matrix representations of the opertors that we diagonalize
+n = 3000
+λ = eigvals(SEg, n);
+
+# We retain a fraction of the eigenvalues with the least magnitude
+λ = λ[1:round(Int, 3n/5)];

examples/Eigenvalue_standard.jl

@@ -9,7 +9,7 @@
 # where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding
 # *eigenfunction* ``\mathop{u}_i(x)``.
 
-### Quantum harmonic oscillator
+# ### Quantum harmonic oscillator
 
 # A classic eigenvalue problem is known as the quantum harmonic oscillator where
 # ```math

examples/Eigenvalue_well_barrier.jl

@@ -0,0 +1,69 @@
+# ### Infinite well with a barrier
+
+# We solve the Schrodinger equation in an infinite square well in the domain `-Lx/2..Lx/2`, with a finite barrier in the middle from `-d/2..d/2`
+# ```math
+# \mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + V(x),
+# ```
+# where ``V(x)`` is given by an infinite well with a smoothed barrier at the center.
+
+# We note that the system has parity symmetry, so the solutions may be separated into odd and even functions.
+# We may therefore solve the problem only for half the domain, with Dirichlet boundary conditions at the midpoint
+# for odd functions, and Neumann conditions for even functions.
+# This projection to subspaces allows us to halve the sizes of matrices that we need to diagonalize
+
+Lx = 4 # size of the domain
+S = Legendre(0..Lx/2)
+x = Fun(S)
+d = 1 # size of the barrier
+Δ = 0.1 # smoothing scale of the barrier
+V = 50 * (1 - tanh((x - d/2)/Δ))/2; # right half of the potential barrier
+H = -Derivative(S)^2/2 + V;
+
+# Odd solutions, with a zero Dirichlet condition at `0` representing a node
+B = Dirichlet(S);
+Seig = ApproxFun.SymmetricEigensystem(H, B);
+
+# Diagonalize `n × n` matrix representations of the basis-recombined operators
+# We specify a tolerance to reject spurious solutions arising from the discretization
+n = 100
+λodd, vodd = ApproxFun.eigs(Seig, n, tolerance=1e-8);
+
+# To extend the solutions to the full domain, we construct the left-half space.
+Scomplement = Legendre(-Lx/2..0);
+
+# We use the fact that Legendre polynomials of odd orders are odd functions,
+# and those of even orders are even functions
+# Using this, for the odd solutions, we negate the even-order coefficients to construct the odd image in `-Lx/2..0`
+function oddimage(f, Scomplement)
+	coeffs = [(-1)^isodd(m) * c for (m,c) in enumerate(coefficients(f))]
+	Fun(Scomplement, coeffs)
+end;
+voddimage = oddimage.(vodd, Scomplement);
+
+# Construct the functions over the entire domain `-Lx/2..Lx/2` as piecewise sums over the two half domains `-Lx/2..0` and `0..Lx/2`
+voddfull = voddimage .+ vodd;
+
+# Even solutions, with a Neumann condition at `0` representing the symmetry of the function
+B = [lneumann(S); rdirichlet(S)];
+
+Seig = ApproxFun.SymmetricEigensystem(H, B);
+λeven, veven = ApproxFun.eigs(Seig, n, tolerance=1e-8);
+
+# For the even solutions, we negate the odd-order coefficients to construct the even image in `-Lx/2..0`
+function evenimage(f, Scomplement)
+	coeffs = [(-1)^iseven(m) * c for (m,c) in enumerate(coefficients(f))]
+	Fun(Scomplement, coeffs)
+end;
+vevenimage = evenimage.(veven, Scomplement);
+vevenfull = vevenimage .+ veven;
+
+# We interlace the eigenvalues and eigenvectors to obtain the entire spectrum
+function interlace(a::AbstractVector, b::AbstractVector)
+	vec(permutedims([a b]))
+end;
+λ = interlace(λeven, λodd);
+v = interlace(vevenfull, voddfull);
+
+# Symmetrize the potential using an even image (this is mainly for plotting/post-processing)
+Vevenimage = evenimage(V, Scomplement);
+Vfull = Vevenimage + V;

test/ReadmeTest.jl

@@ -9,6 +9,11 @@ const EXAMPLES_DIR = joinpath(dirname(dirname(pathof(ApproxFun))), "examples")
 
 include(joinpath(@__DIR__, "testutils.jl"))
 
+function get_included_files(filename)
+    v = [l for l in eachline(joinpath(EXAMPLES_DIR, filename)) if contains(l, "include")]
+    strip.(getindex.(split.(getindex.(split.(v, "include("), 2), ")"), 1), '\"')
+end
+
 @verbose @testset "Readme" begin
     @testset "Calculus and algebra" begin
         x = Fun(identity,0..10)
@@ -58,34 +63,44 @@ include(joinpath(@__DIR__, "testutils.jl"))
     end
 
     @testset "NonlinearBVP" begin
-        include(joinpath(EXAMPLES_DIR, "NonlinearBVP1.jl"))
-        include(joinpath(EXAMPLES_DIR, "NonlinearBVP2.jl"))
+        for f in get_included_files("NonlinearBVP.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
 
     @testset "ODE" begin
-        include(joinpath(EXAMPLES_DIR, "ODE_BVP.jl"))
-        include(joinpath(EXAMPLES_DIR, "ODE_increaseprec.jl"))
+        for f in get_included_files("ODE.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
     @testset "System of Equations" begin
-        include(joinpath(EXAMPLES_DIR, "System1.jl"))
+        for f in get_included_files("system_of_eqn.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
 
     @testset "Periodic" begin
-        include(joinpath(EXAMPLES_DIR, "Periodic1.jl"))
+        for f in get_included_files("Periodic.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
     @testset "Sampling" begin
-        include(joinpath(EXAMPLES_DIR, "Sampling1.jl"))
+        for f in get_included_files("Sampling.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
     @testset "Eigenvalue" begin
-        include(joinpath(EXAMPLES_DIR, "Eigenvalue_standard.jl"))
-        include(joinpath(EXAMPLES_DIR, "Eigenvalue_symmetric.jl"))
+        for f in get_included_files("Eigenvalue.jl")
+            include(joinpath(EXAMPLES_DIR, f))
+        end
     end
 
     # this is broken on v1.6 (on BandedArrays < v0.16.16), so we only test on higher versions
     if VERSION >= v"1.7"
         @testset "PDE" begin
-            include(joinpath(EXAMPLES_DIR, "PDE_Poisson.jl"))
-            include(joinpath(EXAMPLES_DIR, "PDE_Helmholtz.jl"))
+            for f in get_included_files("PDE.jl")
+                include(joinpath(EXAMPLES_DIR, f))
+            end
         end
     end