Change default ODE solver to Vern7

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumToolbox"
 uuid = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
-version = "0.38.1"
 authors = ["Alberto Mercurio", "Yi-Te Huang"]
+version = "0.38.1"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -16,7 +16,7 @@ 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"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -59,7 +59,7 @@ LinearAlgebra = "1"
 LinearSolve = "2, 3"
 Makie = "0.24"
 OrdinaryDiffEqCore = "1"
-OrdinaryDiffEqTsit5 = "1"
+OrdinaryDiffEqVerner = "1.5.0"
 Pkg = "1"
 ProgressMeter = "1.11.0"
 Random = "1"

src/QuantumToolbox.jl

@@ -32,7 +32,8 @@ import SciMLBase:
     DiscreteCallback,
     AbstractSciMLProblem,
     AbstractODEIntegrator,
-    AbstractODESolution
+    AbstractODESolution,
+    AbstractODEAlgorithm
 import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA2, SRIW1
 import SciMLOperators:
     cache_operator,
@@ -51,8 +52,7 @@ 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 DiffEqNoiseProcess: RealWienerProcess!, RealWienerProcess
 
 # other dependencies (in alphabetical order)

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 = Vern7(),
         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 `Vern7()`.
 - `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 = Vern7(),
     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 = Vern7(),
         ψ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=Vern7()`: 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 = Vern7()
     ψ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 = Vern7(),
     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 = Vern7(),
         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 `Vern7()`.
 - `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 = Vern7(),
     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 = Vern7(),
     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 = Vern7(),
         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 `Vern7()`.
 - `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 = Vern7(),
     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 = Vern7(); 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 = Vern7(),
         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 `Vern7()`.
 - `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 = Vern7(),
     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 = Vern7(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,

src/time_evolution/sesolve.jl

@@ -97,7 +97,7 @@ end
         H::Union{AbstractQuantumObject{Operator},Tuple},
         ψ0::QuantumObject{Ket},
         tlist::AbstractVector;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = Vern7(lazy=false),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         params = NullParameters(),
         progress_bar::Union{Val,Bool} = Val(true),
@@ -116,7 +116,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ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.
-- `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
+- `alg`: The algorithm for the ODE solver. The default is `Vern7(lazy=false)`.
 - `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.
@@ -139,7 +139,7 @@ function sesolve(
     H::Union{AbstractQuantumObject{Operator},Tuple},
     ψ0::QuantumObject{Ket},
     tlist::AbstractVector;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = Vern7(lazy=false),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -171,7 +171,7 @@ function sesolve(
     end
 end
 
-function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(); kwargs...)
+function sesolve(prob::TimeEvolutionProblem, alg::AbstractODEAlgorithm = Vern7(lazy=false); kwargs...)
     sol = solve(prob.prob, alg; kwargs...)
 
     return _gen_sesolve_solution(sol, prob.times, prob.dimensions)
@@ -182,7 +182,7 @@ end
         H::Union{AbstractQuantumObject{Operator},Tuple},
         ψ0::Union{QuantumObject{Ket},AbstractVector{<:QuantumObject{Ket}}},
         tlist::AbstractVector;
-        alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+        alg::AbstractODEAlgorithm = Vern7(lazy=false),
         ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
         e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
         params::Union{NullParameters,Tuple} = NullParameters(),
@@ -205,7 +205,7 @@ for each combination in the ensemble.
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state(s) of the system. Can be a single [`QuantumObject`](@ref) or a `Vector` of initial states.
 - `tlist`: List of time points at which to save either the state or the expectation values of the system.
-- `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
+- `alg`: The algorithm for the ODE solver. The default is `Vern7(lazy=false)`.
 - `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.
@@ -226,7 +226,7 @@ function sesolve_map(
     H::Union{AbstractQuantumObject{Operator},Tuple},
     ψ0::AbstractVector{<:QuantumObject{Ket}},
     tlist::AbstractVector;
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = Vern7(lazy=false),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
     e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
     params::Union{NullParameters,Tuple} = NullParameters(),
@@ -266,7 +266,7 @@ sesolve_map(H::Union{AbstractQuantumObject{Operator},Tuple}, ψ0::QuantumObject{
 function sesolve_map(
     prob::TimeEvolutionProblem{<:ODEProblem},
     iter::AbstractArray,
-    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
+    alg::AbstractODEAlgorithm = Vern7(lazy=false),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,

src/time_evolution/time_evolution.jl

@@ -65,7 +65,7 @@ struct TimeEvolutionSol{
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     RETT<:Enum,
-    AlgT<:OrdinaryDiffEqAlgorithm,
+    AlgT<:AbstractODEAlgorithm,
     TolT<:Real,
 }
     times::TT1
@@ -138,7 +138,7 @@ struct TimeEvolutionMCSol{
     TE<:Union{AbstractArray,Nothing},
     TJT<:Vector{<:Vector{<:Real}},
     TJW<:Vector{<:Vector{<:Integer}},
-    AlgT<:OrdinaryDiffEqAlgorithm,
+    AlgT<:AbstractODEAlgorithm,
     TolT<:Real,
 } <: TimeEvolutionMultiTrajSol{TS,TE}
     ntraj::Int
@@ -406,7 +406,7 @@ end
 
 function _ensemble_dispatch_solve(
     ens_prob::TimeEvolutionProblem,
-    alg::Union{<:OrdinaryDiffEqAlgorithm,<:StochasticDiffEqAlgorithm},
+    alg::Union{<:AbstractODEAlgorithm,<:StochasticDiffEqAlgorithm},
     ensemblealg::ET,
     ntraj::Int,
 ) where {ET<:Union{EnsembleSplitThreads,EnsembleDistributed}}
@@ -427,7 +427,7 @@ function _ensemble_dispatch_solve(
 end
 function _ensemble_dispatch_solve(
     ens_prob::TimeEvolutionProblem,
-    alg::Union{<:OrdinaryDiffEqAlgorithm,<:StochasticDiffEqAlgorithm},
+    alg::Union{<:AbstractODEAlgorithm,<:StochasticDiffEqAlgorithm},
     ensemblealg,
     ntraj::Int,
 )