# Seeking Tips for Solving an ODE with DifferentialEquations

Our package QuantumAnnealing solves ODEs that arise in design of quantum annealing algorithms. The package includes two solution methods `simulate` (a hand crafted solution) and `simulate_de` (a generic solution approach using DifferentialEquations).

Source code here.

We have noticed that there can be correctness issues with our `simulate_de` implementation when running long time horizon problems, which makes me suspect we need to tweak the `solve` settings we are using. To that end we are seeking advise on the following points,

• Are there changes we can make to improve the solution accuracy for long time horizons?
• Are there changes we can make to improve runtime performance?
• Are there any ways to avoid the precompile time of DifferentialEquations when running in CI? Compilation of DifferentialEquations is currently about 80% of our CI runtime.

Here is a MWEs to highlight the correctness issue. Note that `sum(z_measure_probabilities(ρ)` should always be 1.0 in a perfect simulation.

``````using QuantumAnnealing

ising_model = Dict((1,) => 1.0, (1,2) => -1.0)

ρ = simulate_de(ising_model, 1.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 7.214235431263205e-9

ρ = simulate_de(ising_model, 10.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 8.508320654354584e-6

ρ = simulate_de(ising_model, 100.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.0003582731716551546

ρ = simulate_de(ising_model, 1000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.003550520566418469

ρ = simulate_de(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.03483614681737934

ρ = simulate_de(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 0.2990146971495551

ρ = simulate(ising_model, 10000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.3211653993039363e-14

ρ = simulate(ising_model, 100000.0, AS_LINEAR); pr_error = abs(1.0 - sum(z_measure_probabilities(ρ))) # 1.0436096431476471e-13

ρ_target = simulate(ising_model, 10000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 10000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.0348400079153087

ρ_target = simulate(ising_model, 100000.0, AS_LINEAR); ρ_de = simulate_de(ising_model, 100000.0, AS_LINEAR); error = sum(abs.(ρ_target .- ρ_de)) # 0.2990175297085176
``````

CC @zmorrell

I see from the source code that you are just calling `solve(prob)`, thus letting `DifferentialEquations` automagically choose the solver for you.
While that works out well for most simple cases, I think you might want to have more control on that.

I didn’t look into the details of the problem here but I believe, in general, a sensible first step when you have an unstable solution is to just try different solvers and timesteps (adaptive vs fixed), depending on your problem requirements. There’s a lot to choose from: ODE Solvers · DifferentialEquations.jl. It can really make a difference.

Also, a quick plot of the solution vs time rather than just the endpoint might give you some more insight (e.g., is it just straight blowing up, doing crazy oscillations, showing “non-physical” values…?)

I assume you’re tweaking tolerance?

Is your `f` function non-allocating and have you optimized solver choice? Is it stiff? Have you chosen linear solvers etc. effectively?

@ccoffrin I tried lowering the relative tolerance while specifying a non-stiff method and that seems to be giving adequate results for the two qubit test model shown here as well as a 3 qubit test case I ran as well. I suggest that we use Vern9() and fiddle with the default tolerances in the simulate_de function.

2 Likes

Thanks all for out input on this. We were able to get the `solve` method to produce correct in all of the cases I listed above by adjusting the `abstol` and `reltol`. This was a great help, thanks!