# Solving differential equations with uncertain parameters

**URL:** <https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [June 8, 2023, 3:57pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059 "2023-06-08T15:57:30Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![lmg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmg/32/32312_2.png) [@lmg](https://discourse.julialang.org/u/lmg)\
**Post date:** [June 8, 2023, 3:57pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/1 "2023-06-08T15:57:30Z")

</div>

Hello everybody,

I would like to solve IVPs with uncertain but bounded parameters. To achieve this, my plan was to use the `Interval` data type from the IntervalArithmetic.jl package. In [this post on custom number types](https://docs.sciml.ai/SciMLTutorialsOutput/html/type_handling/01-number_types.html), it is stated that DifferentialEquations.jl works with some Julia-defined types, such as ArbFloats.jl. However, this plan didn’t work out and I can’t find anything online on the compatibility of DifferentialEquations.jl and IntervalArithmetic.jl.

Here’s a minimial (not) working example, but it works with crisp parameters `p` (I’m using Julia 1.9.1 with DifferentialEquations.jl v7.8.0 and IntervalArithmetic v0.20.8):

```julia
using DifferentialEquations, IntervalArithmetic

function model(du, u, p, t)
   du[1] = -(p[3] + p[1]) * u[1] + (p[2] * u[2])
   du[2] = (p[1] * u[1]) - (p[2] * u[2])
end

tspan = (0.0, 10.0)
u0 = [1.0, 0.0]
p = [Interval(1.9, 2.0), Interval(0.2, 0.3), Interval(0.1, 0.2)]

prob = ODEProblem(model, u0, tspan, p)
sol = solve(prob, RK4())

```

This fails because of of the following error:

```julia
ERROR: MethodError: no method matching Float64(::Interval{Float64})

```

Why are `Interval`s converted to `Float64`? A conversion back to a floating point value would void the reason of using intervals in the first place, because I need to preserve the lower and upper bounds.

Are these two packages just not compatible, or am I doing something fundamentally wrong?

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 8, 2023, 4:00pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/2 "2023-06-08T16:00:24Z")

</div>

If I remember correctly, for this kind of stuff you may have to use uncertainty on `u0` and `tspan` as well

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 8, 2023, 4:01pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/3 "2023-06-08T16:01:23Z")

</div>

See [DifferentialEquations.jl and Measurements.jl](https://discourse.julialang.org/t/differentialequations-jl-and-measurements-jl/6350) for a related question. (In that thread, they noted the the timespan needed to be modified too, because of adaptive timestepping.)

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [June 8, 2023, 4:59pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/4 "2023-06-08T16:59:41Z")

</div>

1. The integrator assumes the type of `u` and `du` to be the same. However, multiplication by `p` will result in an `Interval`. Changing the initial conditions to also be intervals alleviates this problem.
2. The time steps also need to be intervals as mentioned.

In general, this would solve this type of issue, but `Intervals` will still not work. Although `Interval <: Number` it really represents a set. The solvers will require determining calling functions like `nextfloat(::Interval)` which doesn’t exist, and I don’t know if it makes sense anyway.

Additionally, if we could get past these types of issues, I suspect you would get a useless answer. Interval arithmetic is guaranteed to include the true propagated set, but suffers from a bloating of the interval that would most likely produces [-∞, ∞] for a state very quickly in an iterative algorithm like and ODE solve.

This stems from something called the [dependency problem](https://en.wikipedia.org/wiki/Interval_arithmetic#Dependency_problem) and the wrapping effect of interval arithmetic (essentially iteratively hitting the dependency problem).

Not all hope is lost though, checkout [GitHub - JuliaReach/ReachabilityAnalysis.jl: Methods to compute sets of states reachable by dynamical systems](https://github.com/JuliaReach/ReachabilityAnalysis.jl) for another approach.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [June 8, 2023, 5:07pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/5 "2023-06-08T17:07:11Z")

</div>

Here is a simple example that highlights the wrapping effect. This applies a pure rotation of -45 deg of line segment bounded by an interval box.

```julia
A = [0 1; -1 -1]
βint = -0.5 .. 0.5

x0 = [βint, -βint]

A*(A*(A*x0)) # [-1.5, 1.5] x [-2.5, 2.5]
A^3*x0 #[-0.5, 0.5] x[-0.5, 0.5]

```

The difference in the two approaches above is that the second delays the interval arithmetic and has a single interval vector multiplication vs. 3

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/2/f23407a08803c4b34395110a4cc647578cbec8b9.png)

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [June 8, 2023, 7:20pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/6 "2023-06-08T19:20:21Z")

</div>

You may want to have a look at [MonteCarloMeasurements.jl](https://baggepinnen.github.io/MonteCarloMeasurements.jl/latest/examples/#Differential-Equations-1), and possibly also [ReachabilityAnalysis.jl](https://juliareach.github.io/ReachabilityAnalysis.jl/dev/).

---

<div class="post-metadata">

**Author:** ![lmg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmg/32/32312_2.png) [@lmg](https://discourse.julialang.org/u/lmg)\
**Post date:** [June 9, 2023, 7:21am UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/7 "2023-06-09T07:21:54Z")

</div>

Thanks for the explanation and illustrative example!

I am aware of the downsides of using interval arithmetic. Nonetheless, I was curious if one could use intervals with DifferentialEquations.jl, because there are verified solvers like [DynIBEX](https://perso.ensta-paris.fr/~chapoutot/dynibex/) available for other languages and I wanted to know if DifferentialEquations.jl was able to provide verified solutions using intervals, too. I will have a look at the other suggestions from this thread.

ReachabilityAnalysis.jl seems like a suitable package for the original task, though.

---

<div class="post-metadata">

**Author:** ![DrPapa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drpapa/32/6835_2.png) [@DrPapa](https://discourse.julialang.org/u/DrPapa)\
**Post date:** [June 9, 2023, 1:10pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/8 "2023-06-09T13:10:41Z")

</div>

If you had something like a fixed time step RK method you could do verified integration if using intervals for time and state.

ReachabilityAnalysis.jl is probably what you want though for verified solutions to an ODE. Also take a look at TaylorModels.jl. It isn’t well documented but there is a [verified solver](https://github.com/JuliaIntervals/TaylorModels.jl/blob/b1697ac5d0c1d01d390fc54a7b3f8f719d1ae63f/src/validatedODEs.jl#LL1237C10-L1237C26) in there based on Picard iteration ([details](https://theses.hal.science/tel-00657843)). This is what ReachabilityAnalysis.jl uses under the hood for the `TMJet` algorithms.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [June 16, 2023, 11:31pm UTC](https://discourse.julialang.org/t/solving-differential-equations-with-uncertain-parameters/100059/9 "2023-06-16T23:31:33Z")

</div>

Here’s a worked out example I wrote a while ago \> [Introduction · ReachabilityAnalysis.jl](https://juliareach.github.io/ReachabilityAnalysis.jl/dev/tutorials/taylor_methods/introduction/)

For a comparison of different tools you may check [ARCH-COMP22 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics](https://easychair.org/publications/paper/JrQ4)

For nonlinear ODEs, besides the `TMJets` solver, thanks to agerlach it’s possible to use the Flow\* solver too.