Eigenvalue problem

Standard eigenvalue problem

In analogy to linear algebra, many differential equations may be posed as eigenvalue problems. That is, for some differential operator $\mathop{L}$, there are a family of functions $\mathop{u}_i(x)$ such that

\[\mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),\]

where $λ_i$ is the $i^{th}$ eigenvalue of the $L$ and has a corresponding eigenfunction $\mathop{u}_i(x)$.

Quantum harmonic oscillator

A classic eigenvalue problem is known as the quantum harmonic oscillator where

\[\mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,\]

and one demands that $\mathop{u}(∞) = \mathop{u}(-∞) = 0$. Because we expect the solutions to be exponentially suppressed for large $x$, we can approximate this with Dirichlet boundary conditions at a 'reasonably large' $x$ without much difference.

We can express this in ApproxFun as the following:

using ApproxFun

x = Fun(-8..8)
V = x^2/2
L = -𝒟^2/2 + V
S = space(x)
B = Dirichlet(S)
λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10);

Note that boundary conditions must be specified in the call to eigs. Plotting the first $20$ eigenfunctions offset vertically by their eigenvalue, we see

import Plots
using LinearAlgebra: norm
p = Plots.plot(V, legend=false, ylim=(-Inf, λ[22]))
for k=1:20
    Plots.plot!(real(v[k]) + real(λ[k]))
end
p

If the solutions are not relatively constant near the boundary then one should push the boundaries further out.

For problems with different constraints or boundary conditions, B can be any zero functional constraint, e.g., DefiniteIntegral().

Self-adjoint Eigenvalue Problem

Ref: J. L. Aurentz & R. M. Slevinsky (2019), arXiv:1903.08538

We solve the confined anharmonic oscillator

\[\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 operators that we diagonalize

n = 3000
λ = eigvals(SEg, n);

We retain a fraction of the eigenvalues with the least magnitude

λ = λ[1:round(Int, 3n/5)];

We plot the eigenvalues

import Plots
Plots.plot(λ, title = "Eigenvalues", legend=false)

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

\[\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 The eigenfunctions vodd are normalized on the half-domain, so we normalize the sum by dividing it by √2

voddfull = (voddimage .+ vodd)./√2;

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)./√2;

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;

We plot the first few eigenfunctions offset by their eigenvalues. The eigenfunctions appear in odd-even pairs by construction.

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]) + λ[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]) + λ[1], linewidth=2)

Plots.plot(p1, p2, layout = 2)

This page was generated using Literate.jl.