Solving differential equations with uncertain parameters

Specific Domains Modelling & Simulations
1

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, 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):

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:

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

Why are Intervals 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!

2

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

3

See DifferentialEquations.jl and Measurements.jl for a related question. (In that thread, they noted the the timespan needed to be modified too, because of adaptive timestepping.)

1 Like
4
  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 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 for another approach.

3 Likes
5

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.

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
image1794×1277 85.3 KB

1 Like
6

You may want to have a look at MonteCarloMeasurements.jl, and possibly also ReachabilityAnalysis.jl.

7

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

8

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 in there based on Picard iteration (details). This is what ReachabilityAnalysis.jl uses under the hood for the TMJet algorithms.

4 Likes
9

Here’s a worked out example I wrote a while ago > Introduction · ReachabilityAnalysis.jl

For a comparison of different tools you may check ARCH-COMP22 Category Report: Continuous and Hybrid Systems with Nonlinear Dynamics

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