Setting `abstol` and `reltol` when solving the Schrödinger equation with OrdinaryDiffEq

Specific Domains Quantum
,
4

Indeed, I found while playing around with this that the tstops here actually hurt the accuracy, and asking for an abstol also hurt accuray (rather than a reltol).

In addition, Vern8() tends to be a more accurate solver for this sort of problem since it’s symplectic and thus will keep the Hamiltonian on the symplectic manifold, conserving various properties. (edit: you can of course do better with the sorts of solvers Chris mentioned above, but Vern8 does better than Tsit5 here and is easier to switch out).

Here’s what I see for the total norm of the wave-function over time after making some adjustments:

using OrdinaryDiffEq: OrdinaryDiffEq as ODE, ODEProblem, solve as ode_solve
using Plots

function f_ode(u, p, t)
	Δ = 1.0
	α = 4.6
	H = Hermitian([
        0.5 * α * t   Δ
        Δ            -0.5 * α * t
    ])
    -im * (H * u)
end

function tst(;tlist = range(-100, 100; length=1001), alg=ODE.Vern8())

    Ψ₀ = [-0.999998+0.0im, 0.0021739+0.0im]  # eigenstate at t = -100

    ode_problem_direct = ODEProblem(f_ode, Ψ₀, (tlist[begin], tlist[end]))

    Ψ = ode_solve(ode_problem_direct, alg;
                  reltol=1e-10)
    norms = map(t -> norm(Ψ(t)), tlist)
    plot(tlist, norms;
         label="⟨Ψ⟩²",
         legend=:right, title="Projection onto canonical basis states", xlabel="time", ylabel="population")

end

image
image599×399 21.2 KB

which I think is a pretty impressive level of norm-conservation.

If it becomes a problem beyond that though, you can do tricks like using callbacks to re-normalize your solution.