Complete Example: The Harmonic Oscillator

We will now turn away from general relativity and use QuasinormalModes.jl to compute to compute the energy eigenvalues of the quantum harmonic oscillator following this paper.

Mathematical preliminaries

If we measure the energy of the system in units of $\hbar\omega$ and distance in units of $\sqrt{\hbar/(m\omega)}$ the time independent Schrödinger equation for the quantum harmonic oscillator is written as

\[-\psi^{\prime\prime}(x) + x^2\psi(x) = \epsilon\psi(x),\]

where we defined $\epsilon \equiv 2 E$ and $E$ is the quantum state's energy. Imposing that $\psi(x)$ decays like a Gaussian distribution asymptotically, we apply the ansatz

\[\psi(x) = e^{-x^2/2}f(x)\]

which, substituting in the original equation, yields

\[f^{\prime\prime}(x) = 2 x f^\prime(x) + (1-\epsilon)f(x).\]

This allows us to easily identify $\lambda_0 = 2x$ and $s_0 = 1 - \epsilon$. In all our implementations we shall refer the sought eigenvalue $\epsilon$ using the variable ω in order to maintain consistency with the previous example.

Implementing the master equation as an analytic problem

The first step is to load the required packages to run this example: QuasinormalModes and SymEngine:

using QuasinormalModes
using SymEngine

Next, we create a parametric type that sub-types AnalyticAIMProblem. As the eigenvalue in the master equation is a quadratic polynomial, we will sub-type QuadraticEigenvalueProblem with the following structure:

struct HarmonicOscilatorData{N,T} <: QuadraticEigenvalueProblem{N,T}
    nIter::N
    x0::T

    vars::Tuple{Basic, Basic}
    exprs::Tuple{Basic, Basic}
end

Now we implement the constructor and extend the default implementations:

function HarmonicOscilatorData(nIter::N, x0::T) where {N,T}
	
    vars = @vars x ω

    λ0 = 2*x
    S0 = 1 - ω

    return HarmonicOscilatorData{N,T}(nIter, x0, vars, (λ0, S0))
end

QuasinormalModes.λ0(d::HarmonicOscilatorData{N,T}) where {N,T} = d.exprs[1]
QuasinormalModes.S0(d::HarmonicOscilatorData{N,T}) where {N,T}  = d.exprs[2]

QuasinormalModes.get_niter(d::HarmonicOscilatorData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::HarmonicOscilatorData{N,T}) where {N,T} = d.x0

QuasinormalModes.get_ODEvar(d::HarmonicOscilatorData{N,T}) where {N,T} = d.vars[1]
QuasinormalModes.get_ODEeigen(d::HarmonicOscilatorData{N,T}) where {N,T} = d.vars[2]

Implementing the master equation as a numeric problem

The structure, constructor and extensions are

struct NHarmonicOscilatorData{N,T} <: NumericAIMProblem{N,T}
    nIter::N
    x0::T
end

function NHarmonicOscilatorData(nIter::N, x0::T) where {N,T}
    return NHarmonicOscilatorData{N,T}(nIter, x0)
end

QuasinormalModes.λ0(::NHarmonicOscilatorData{N,T}) where {N,T} = (x,ω) -> 2*x
QuasinormalModes.S0(::NHarmonicOscilatorData{N,T}) where {N,T} = (x,ω) -> 1 - ω + x - x

QuasinormalModes.get_niter(d::NHarmonicOscilatorData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::NHarmonicOscilatorData{N,T}) where {N,T} = d.x0

Constructing problems and initializing the cache

Once again, we create our problems and cache objects by calling the constructors:

p_ana = HarmonicOscilatorData(0x0000A, 0.5);
p_num = NHarmonicOscilatorData(0x0000A, 0.5);

c_ana = AIMCache(p_ana)
c_num = AIMCache(p_num)

Here we are setting up problems to be solved using 10 iterations with x0 = 0.5

Computing the eigenvalues

Once again we compute the eigenvalues by calling

ev_ana = computeEigenvalues(Serial(), p_ana, c_ana)
ev_num = eigenvaluesInGrid(Serial(), p_num, c_num, (0.0, 21.0))

The results are two arrays, containing the eigenvalues. As before, we define a function to print the results to stdout

function printEigen(eigenvalues)
    println("--------------------------------------")

    for i in eachindex(eigenvalues)
        println("n = $i, ω = $(eigenvalues[i])")
    end
    
    println("--------------------------------------")

    return nothing
end

println("Analytic results")
printEigen(reverse!(ev_ana))

println("Numeric results")
printEigen(ev_num)

The complete source file for this example can be found in harmonic_oscillator.jl. The output is agreement with the expected result for the eigenenergies of the harmonic oscillator, that is, $E_n = n + 1/2$

Analytic results
--------------------------------------
n = 1, ω = 0.9999999999999999 + 0.0im
n = 2, ω = 2.9999999999999964 + 0.0im
n = 3, ω = 4.999999999999426 + 0.0im
n = 4, ω = 7.000000000006788 + 0.0im
n = 5, ω = 8.999999999980533 + 0.0im
n = 6, ω = 10.999999804542819 + 0.0im
n = 7, ω = 13.000000959453153 + 0.0im
n = 8, ω = 14.999998108295404 + 0.0im
n = 9, ω = 17.00000187312756 + 0.0im
n = 10, ω = 18.999999068409203 - 0.0im
n = 11, ω = 21.000000186185098 + 0.0im
--------------------------------------
Numeric results
--------------------------------------
n = 1, ω = 1.0
n = 2, ω = 3.000000000000006
n = 3, ω = 5.000000000000002
n = 4, ω = 7.000000000000006
n = 5, ω = 8.999999999999988
n = 6, ω = 11.0
n = 7, ω = 12.999999999999977
n = 8, ω = 15.000000000000014
n = 9, ω = 16.999999999999908
n = 10, ω = 19.000000000000025
n = 11, ω = 21.0
--------------------------------------