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

**URL:** <https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766>\
**Category:** Quantum\
**Created:** [August 4, 2021, 12:47am UTC](https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766 "2021-08-04T00:47:26Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Noel\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noel_araujo/32/4862_2.png) [@Noel\_Araujo](https://discourse.julialang.org/u/Noel_Araujo)\
**Post date:** [August 4, 2021, 12:47am UTC](https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766/1 "2021-08-04T00:47:26Z")

</div>

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](https://en.wikipedia.org/wiki/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.

```julia
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)`

```julia
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
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 😫?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 5, 2021, 12:31am UTC](https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766/2 "2021-08-05T00:31:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jagot](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jagot/32/12217_2.png) [@jagot](https://discourse.julialang.org/u/jagot)\
**Post date:** [August 5, 2021, 8:22am UTC](https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766/3 "2021-08-05T08:22:46Z")

</div>

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?](https://discourse.julialang.org/t/how-to-efficiently-exponentiate-a-complex-symmetric-tridiagonal-matrix/56982)

Are you sure you want this Hamiltonian matrix though?

> [@Noel\_Araujo](#):
>
> `H = [1 1;-1 1]./sqrt(2)`

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

---

<div class="post-metadata">

**Author:** ![Noel\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/noel_araujo/32/4862_2.png) [@Noel\_Araujo](https://discourse.julialang.org/u/Noel_Araujo)\
**Post date:** [August 5, 2021, 12:20pm UTC](https://discourse.julialang.org/t/how-do-i-create-a-quantum-time-operator-with-differentialequations/65766/4 "2021-08-05T12:20:21Z")

</div>

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.

```julia
(...)
β₀ = [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)

```

![toDiscourse](https://global.discourse-cdn.com/julialang/original/3X/0/b/0bb1567ffee038db1c5dd7c46191a230a42c5f4a.png)

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