# How do I create a Quantum Time Operator with DifferentialEquations?

For the special case of unitary matrices, my quantum mechanics books tells me to use Time Evolution Operators.
`β(t) = exp(-im*H*t/h)*β(0)`

As far I know, I just need to use Matrix Exponential.
First, decompose `H = ψ*λ*inv(ψ)`

Second, we create a time evolution operator `U = exp(-im*H*t) = ψ*exp(-im*λ*t)*inv(ψ)`

If one sets a small enough time interval `δt` : `δU = exp(-im*H*δt) = ψ*exp(-im*λ*δt)*inv(ψ)`

Then, I can evolve iteratively the solution `βₜ₊₁ = δU*βₜ` of my problem.

``````using LinearAlgebra
using Random
Random.seed!(1256)
N = 2

H = [1 1;-1 1]./sqrt(2)
β₀= rand(ComplexF64, N);

maxTimeSteps, t_max = 1000, 10
t = range(1, stop=t_max, length=maxTimeSteps);
δₜ = t[2]-t[1];

@time begin
λ, ψ = eigen(H);

δU = ψ*exp(-im*Diagonal(λ*δₜ))*inv(ψ)

βₜ = []
push!(βₜ,  deepcopy(β₀));
for idx_t = 2:maxTimeSteps
push!(βₜ, δU*βₜ[idx_t-1]  )
end
end
``````

I would like to solve this problem with `DifferentialEquations`, but I don’t know how to properly write down. It should be simple, but I’m taking time.
My bet was to use `du/dt = a*u` as input, because its solution is `u(t)=exp(a*t)`

``````using DifferentialEquations
function ff(du, u, p,t)
du[:] .= -im*H*u
end
tspan = (0.0,t_max);
prob = ODEProblem(ff,β₀,tspan);
@time sol = solve(prob, Tsit5(), reltol=δₜ, abstol=δₜ);
``````

I can check the difference of the last results are far away of each other.

``````julia> abs.(βₜ[end] .- sol.u[end])
2-element Vector{Float64}:
424.06945981919637
424.06825305239903
``````

Is there any one with physics background that already knows where my silly mistake is ?

1 Like

See the examples which use the exponential integrators and Magnus methods.

2 Likes

Chris is right about Magnus propagators being the method of choice for QM. See also this earlier post: How to efficiently exponentiate a complex (symmetric) tridiagonal matrix?

Are you sure you want this Hamiltonian matrix though?

This is non-Hermitian, which is very untypical in Markovian (=normal) QM.

2 Likes

I knew it. I did some silly mistake, the wrong Hamiltonian was the first.
The second was to initialize the first method at `t_0=1`, and the other with `t_0=0`.

BUT I found major mistake :

I setup `reltol=δₜ, abstol=δₜ`, and in this example, `δₜ=0.10010`, which is a terrible tolerance. Just using default values, or not adaptive methods, everything worked as expected.

``````(...)
β₀ = [1.0+0im,1.0+0im] ./sqrt(2)

(...)
@time sol1 = solve(prob, Tsit5(), adaptive=true, dt=δₜ, dense=true, reltol=δₜ, abstol=δₜ);
@time sol2 = solve(prob, Tsit5(), adaptive=true, dt=δₜ, dense=true);
@time sol3 = solve(prob, Tsit5(), adaptive=false, dt=δₜ);

using Plots
fig1 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol1(t[i])[1]) for i=1:length(t)], label="sol1", lw=4, linestyle=:dot)

fig2 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol2(t[i])[1]) for i=1:length(t)], label="sol2", lw=4, linestyle=:dot)

fig3 = plot(t,  [abs2(βₜ[i][1]) for i=1:length(t)], label="exp")
plot!(t,  [abs2(sol3(t[i])[1]) for i=1:length(t)], label="sol3", lw=4, linestyle=:dot)

plot(fig1, fig2, fig3)
``````

Anyway, using the matrix exponential still faster approach if I don’t need to guarantee stability.

2 Likes