# Seeking Tips for Solving an ODE with DifferentialEquations

**URL:** <https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, ode\
**Created:** [February 12, 2022, 6:32pm UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311 "2022-02-12T18:32:56Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [February 12, 2022, 6:32pm UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311/1 "2022-02-12T18:32:56Z")

</div>

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](https://github.com/lanl-ansi/QuantumAnnealing.jl/blob/main/src/simulate_de.jl).

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.

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

---

<div class="post-metadata">

**Author:** ![mastrof](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mastrof/32/17443_2.png) [@mastrof](https://discourse.julialang.org/u/mastrof)\
**Post date:** [February 14, 2022, 1:21am UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311/2 "2022-02-14T01:21:15Z")

</div>

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](https://diffeq.sciml.ai/stable/solvers/ode_solve/). 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…?)

---

<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:** [February 15, 2022, 6:35pm UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311/3 "2022-02-15T18:35:36Z")

</div>

> [@ccoffrin](#):
>
> - Are there changes we can make to improve the solution accuracy for long time horizons?

I assume you’re tweaking tolerance?

> [@ccoffrin](#):
>
> - Are there changes we can make to improve runtime performance?

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

---

<div class="post-metadata">

**Author:** ![zmorrell](https://avatars.discourse-cdn.com/v4/letter/z/a5b964/32.png) [@zmorrell](https://discourse.julialang.org/u/zmorrell)\
**Post date:** [February 16, 2022, 4:16pm UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311/4 "2022-02-16T16:16:09Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [February 19, 2022, 4:45pm UTC](https://discourse.julialang.org/t/seeking-tips-for-solving-an-ode-with-differentialequations/76311/5 "2022-02-19T16:45:33Z")

</div>

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!