improve accuracy for ODE solvers

Author: ytdHuang

Project.toml

@@ -5,7 +5,6 @@ authors = ["Alberto Mercurio", "Yi-Te Huang"]
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -16,7 +15,8 @@ LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -45,7 +45,6 @@ QuantumToolboxMakieExt = "Makie"
 ArrayInterface = "6, 7"
 CUDA = "5.0 - 5.8"
 ChainRulesCore = "1"
-DiffEqBase = "6"
 DiffEqCallbacks = "4.2.1 - 4"
 DiffEqNoiseProcess = "5"
 Distributed = "1"
@@ -59,7 +58,8 @@ LinearAlgebra = "1"
 LinearSolve = "2, 3"
 Makie = "0.24"
 OrdinaryDiffEqCore = "1"
-OrdinaryDiffEqTsit5 = "1"
+OrdinaryDiffEqLowOrderRK = "1"
+OrdinaryDiffEqVerner = "1"
 Pkg = "1"
 ProgressMeter = "1.11.0"
 Random = "1"

src/QuantumToolbox.jl

@@ -1,12 +1,16 @@
 module QuantumToolbox
 
-# Standard Julia libraries
+## Standard Julia libraries
 using LinearAlgebra
-import LinearAlgebra: checksquare
 using SparseArrays
+
+import Distributed: RemoteChannel
+import LinearAlgebra: checksquare
+import Pkg
+import Random: AbstractRNG, default_rng, seed!
 import Statistics: mean, std
 
-# SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
+## SciML packages (for QobjEvo, OrdinaryDiffEq, and LinearSolve)
 import SciMLBase:
     solve,
     solve!,
@@ -32,8 +36,10 @@ import SciMLBase:
     DiscreteCallback,
     AbstractSciMLProblem,
     AbstractODEIntegrator,
-    AbstractODESolution
-import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA2, SRIW1
+    AbstractODEAlgorithm,
+    AbstractODESolution,
+    AbstractSDEAlgorithm
+import StochasticDiffEq: SRA2, SRIW1
 import SciMLOperators:
     cache_operator,
     iscached,
@@ -49,43 +55,41 @@ import SciMLOperators:
     concretize
 import LinearSolve:
     SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES, UMFPACKFactorization, OperatorAssumptions
-import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, FunctionCallingCallback, FunctionCallingAffect, TerminateSteadyState
-import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
-import OrdinaryDiffEqTsit5: Tsit5
+import OrdinaryDiffEqVerner: Vern7
+import OrdinaryDiffEqLowOrderRK: DP5
 import DiffEqNoiseProcess: RealWienerProcess!, RealWienerProcess
 
-# other dependencies (in alphabetical order)
+## other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
-import Distributed: RemoteChannel
 import FFTW: fft, ifft, fftfreq, fftshift
 import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
 import LaTeXStrings: @L_str
-import Pkg
 import ProgressMeter: Progress, next!
-import Random: AbstractRNG, default_rng, seed!
 import SpecialFunctions: loggamma
 import StaticArraysCore: SVector, MVector
 
 # Export functions from the other modules
 
-# LinearAlgebra
+## LinearAlgebra
 export ishermitian, issymmetric, isposdef, dot, tr, svdvals, norm, normalize, normalize!, diag, Hermitian, Symmetric
 
-# SparseArrays
+## SparseArrays
 export permute
 
-# SciMLOperators
+## SciMLOperators
 export cache_operator, iscached, isconstant
 
-# Utility
+# Source files
+
+## Utility
 include("settings.jl")
 include("utilities.jl")
 include("versioninfo.jl")
 include("linear_maps.jl")
 
-# Quantum Object
+## Quantum Object
 include("qobj/space.jl")
 include("qobj/energy_restricted.jl")
 include("qobj/dimensions.jl")
@@ -102,7 +106,7 @@ include("qobj/superoperators.jl")
 include("qobj/synonyms.jl")
 include("qobj/block_diagonal_form.jl")
 
-# time evolution
+## time evolution
 include("time_evolution/time_evolution.jl")
 include("time_evolution/callback_helpers/callback_helpers.jl")
 include("time_evolution/callback_helpers/sesolve_callback_helpers.jl")
@@ -119,7 +123,7 @@ include("time_evolution/ssesolve.jl")
 include("time_evolution/smesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
 
-# Others
+## Others
 include("correlations.jl")
 include("wigner.jl")
 include("spin_lattice.jl")
@@ -131,7 +135,7 @@ include("steadystate.jl")
 include("spectrum.jl")
 include("visualization.jl")
 
-# deprecated functions
+## deprecated functions
 include("deprecated.jl")
 
 end

src/qobj/eigsolve.jl

@@ -422,7 +422,7 @@ end
         H::Union{AbstractQuantumObject{HOpType},Tuple},
         T::Real,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = DP5(),
         params::NamedTuple = NamedTuple(),
         ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
         eigvals::Int = 1,
@@ -440,7 +440,7 @@ Solve the eigenvalue problem for a Liouvillian superoperator `L` using the Arnol
 - `H`: The Hamiltonian (or directly the Liouvillian) of the system. It can be a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a tuple of the form supported by [`mesolve`](@ref).
 - `T`: The time at which to evaluate the time evolution.
 - `c_ops`: A vector of collapse operators. Default is `nothing` meaning the system is closed.
-- `alg`: The differential equation solver algorithm. Default is `Tsit5()`.
+- `alg`: The differential equation solver algorithm. Default is `DP5()`.
 - `params`: A `NamedTuple` containing the parameters of the system.
 - `ρ0`: The initial density matrix. If not specified, a random density matrix is used.
 - `eigvals`: The number of eigenvalues to compute.
@@ -465,7 +465,7 @@ function eigsolve_al(
     H::Union{AbstractQuantumObject{HOpType},Tuple},
     T::Real,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     params::NamedTuple = NamedTuple(),
     ρ0::AbstractMatrix = rand_dm(prod(H.dimensions)).data,
     eigvals::Int = 1,

src/steadystate.jl

@@ -44,7 +44,7 @@ end
 
 @doc raw"""
     SteadyStateODESolver(
-        alg = Tsit5(),
+        alg = DP5(),
         ψ0 = nothing,
         tmax = Inf,
         terminate_reltol = 1e-4,
@@ -66,7 +66,7 @@ or
 ```
 
 # Arguments
-- `alg::OrdinaryDiffEqAlgorithm=Tsit5()`: The algorithm to solve the ODE.
+- `alg::AbstractODEAlgorithm=DP5()`: The algorithm to solve the ODE.
 - `ψ0::Union{Nothing,QuantumObject}=nothing`: The initial state of the system. If not specified, a random pure state will be generated.
 - `tmax::Real=Inf`: The final time step for the steady state problem.
 - `terminate_reltol` = The relative tolerance for stationary state terminate condition. Default to `1e-4`.
@@ -78,13 +78,13 @@ or
 For more details about the solving `alg`orithms, please refer to [`OrdinaryDiffEq.jl`](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 """
 Base.@kwdef struct SteadyStateODESolver{
-    MT<:OrdinaryDiffEqAlgorithm,
+    MT<:AbstractODEAlgorithm,
     ST<:Union{Nothing,QuantumObject},
     TT<:Real,
     RT<:Real,
     AT<:Real,
 } <: SteadyStateSolver
-    alg::MT = Tsit5()
+    alg::MT = DP5()
     ψ0::ST = nothing
     tmax::TT = Inf
     terminate_reltol::RT = 1e-4

src/time_evolution/lr_mesolve.jl

@@ -25,7 +25,7 @@ struct TimeEvolutionLRSol{
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     RetT<:Enum,
-    AlgT<:OrdinaryDiffEqAlgorithm,
+    AlgT<:AbstractODEAlgorithm,
     TolT<:Real,
     TSZB<:AbstractVector,
     TM<:Vector{<:Integer},
@@ -45,7 +45,7 @@ struct TimeEvolutionLRSol{
 end
 
 lr_mesolve_options_default = (
-    alg = Tsit5(),
+    alg = DP5(),
     progress = true,
     err_max = 0.0,
     p0 = 0.0,

src/time_evolution/mcsolve.jl

@@ -274,7 +274,7 @@ end
         ψ0::QuantumObject{Ket},
         tlist::AbstractVector,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = DP5(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         params = NullParameters(),
         rng::AbstractRNG = default_rng(),
@@ -329,7 +329,7 @@ If the environmental measurements register a quantum jump, the wave function und
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
-- `alg`: The algorithm to use for the ODE solver. Default to `Tsit5()`.
+- `alg`: The algorithm to use for the ODE solver. Default to `DP5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
 - `rng`: Random number generator for reproducibility.
@@ -361,7 +361,7 @@ function mcsolve(
     ψ0::QuantumObject{Ket},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params = NullParameters(),
     rng::AbstractRNG = default_rng(),
@@ -398,7 +398,7 @@ end
 
 function mcsolve(
     ens_prob_mc::TimeEvolutionProblem,
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     ntraj::Int = 500,
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
     keep_runs_results = Val(false),

src/time_evolution/mesolve.jl

@@ -129,7 +129,7 @@ end
         ψ0::QuantumObject,
         tlist::AbstractVector,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = DP5(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         params = NullParameters(),
         progress_bar::Union{Val,Bool} = Val(true),
@@ -155,7 +155,7 @@ where
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
-- `alg`: The algorithm for the ODE solver. The default value is `Tsit5()`.
+- `alg`: The algorithm for the ODE solver. The default value is `DP5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -180,7 +180,7 @@ function mesolve(
     ψ0::QuantumObject{StateOpType},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -225,7 +225,7 @@ function mesolve(
     end
 end
 
-function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); kwargs...)
+function mesolve(prob::TimeEvolutionProblem, alg::AbstractODEAlgorithm = DP5(); kwargs...)
     sol = solve(prob.prob, alg; kwargs...)
 
     return _gen_mesolve_solution(sol, prob.times, prob.dimensions, prob.kwargs.isoperket)
@@ -237,7 +237,7 @@ end
         ψ0::Union{QuantumObject,AbstractVector{<:QuantumObject}},
         tlist::AbstractVector,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = DP5(),
         ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         params::Union{NullParameters,Tuple} = NullParameters(),
@@ -267,7 +267,7 @@ for each combination in the ensemble.
 - `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
-- `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
+- `alg`: The algorithm for the ODE solver. The default is `DP5()`.
 - `ensemblealg`: Ensemble algorithm to use for parallel computation. Default is `EnsembleThreads()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: A `Tuple` of parameter sets. Each element should be an `AbstractVector` representing the sweep range for that parameter. The function will solve for all combinations of initial states and parameter sets.
@@ -290,7 +290,7 @@ function mesolve_map(
     ψ0::AbstractVector{<:QuantumObject{StateOpType}},
     tlist::AbstractVector,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params::Union{NullParameters,Tuple} = NullParameters(),
@@ -357,7 +357,7 @@ mesolve_map(
 function mesolve_map(
     prob::TimeEvolutionProblem{<:ODEProblem},
     iter::AbstractArray,
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = DP5(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,