Uncertainty Quantification

Uncertainty quantification allows a user to identify the uncertainty associated with the numerical approximation given by DifferentialEquations.jl. This page describes the different methods available for quantifying such uncertainties. Note that this requires one of the native Julia solvers like OrdinaryDiffEq.jl, StochasticDiffEq.jl, or DelayDiffEq.jl.

Installation

This functionality does not come standard with DifferentialEquations.jl. To use this functionality, you must install DiffEqUncertainty.jl:

]add DiffEqUncertainty
using DiffEqUncertainty

ProbInts

The ProbInts method for uncertainty quantification involves the transformation of an ODE into an associated SDE where the noise is related to the timesteps and the order of the algorithm. This is implemented into the DiffEq system via a callback function. The first form is:

ProbIntsUncertainty(σ,order,save=true)

σ is the noise scaling factor and order is the order of the algorithm. save is for choosing whether this callback should control the saving behavior. Generally this is true unless one is stacking callbacks in a CallbackSet. It is recommended that σ is representative of the size of the errors in a single step of the equation.

If you are using an adaptive algorithm, the callback

AdaptiveProbIntsUncertainty(order,save=true)

determines the noise scaling automatically using an internal error estimate.

Example 1: FitzHugh-Nagumo

In this example we will determine our uncertainty when solving the FitzHugh-Nagumo model with the Euler() method. We define the FitzHugh-Nagumo model:

function fitz(du,u,p,t)
  V,R = u
  a,b,c = p
  du[1] = c*(V - V^3/3 + R)
  du[2] = -(1/c)*(V -  a - b*R)
end
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
p = (0.2,0.2,3.0)
prob = ODEProblem(fitz,u0,tspan,p)

Now we define the ProbInts callback. In this case, our method is the Euler method and thus it is order 1. For the noise scaling, we will try a few different values and see how it changes. For σ=0.2, we define the callback as:

cb = ProbIntsUncertainty(0.2,1)

This is akin to having an error of approximately 0.2 at each step. We now build and solve a EnsembleProblem for 100 trajectories:

ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Euler(),trajectories=100,callback=cb,dt=1/10)

Now we can plot the resulting Monte Carlo solution:

using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

If we increase the amount of error, we see that some parts of the equation have less uncertainty than others. For example, at σ=0.5:

cb = ProbIntsUncertainty(0.5,1)
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Euler(),trajectories=100,callback=cb,dt=1/10)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

But at this amount of noise, we can see how we contract to the true solution by decreasing dt:

cb = ProbIntsUncertainty(0.5,1)
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Euler(),trajectories=100,callback=cb,dt=1/100)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

Example 2: Adaptive ProbInts on FitzHugh-Nagumo

While the first example is academic and shows how the ProbInts method scales, the fact that one should have some idea of the error in order to calibrate σ can lead to complications. Thus the more useful method in many cases is the AdaptiveProbIntsUncertainty version. In this version, no σ is required since this is calculated using an internal error estimate. Thus this gives an accurate representation of the possible error without user input.

Let's try this with the order 5 Tsit5() method on the same problem as before:

cb = AdaptiveProbIntsUncertainty(5)
sol = solve(prob,Tsit5())
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Tsit5(),trajectories=100,callback=cb)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

In this case, we see that the default tolerances give us a very good solution. However, if we increase the tolerance a lot:

cb = AdaptiveProbIntsUncertainty(5)
sol = solve(prob,Tsit5())
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Tsit5(),trajectories=100,callback=cb,abstol=1e-3,reltol=1e-1)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

we can see that the moments just after the rise can be uncertain.

Example 3: Adaptive ProbInts on the Lorenz Attractor

One very good use of uncertainty quantification is on chaotic models. Chaotic equations diverge from the true solution according to the error exponentially. This means that as time goes on, you get further and further from the solution. The ProbInts method can help diagnose how much of the timeseries is reliable.

As in the previous example, we first define the model:

function g(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,30.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(g,u0,tspan,p)

and then we build the ProbInts type. Let's use the order 5 Tsit5 again.

cb = AdaptiveProbIntsUncertainty(5)

Then we solve the MonteCarloProblem

ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Tsit5(),trajectories=100,callback=cb)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

Here we see that by t about 22 we start to receive strong deviations from the "true" solution. We can increase the amount of time before error explosion by using a higher order method with stricter tolerances:

tspan = (0.0,40.0)
prob = ODEProblem(g,u0,tspan,p)
cb = AdaptiveProbIntsUncertainty(7)
ensemble_prob = EnsembleProblem(prob)
sim = solve(ensemble_prob,Vern7(),trajectories=100,callback=cb,reltol=1e-6)
using Plots; plotly(); plot(sim,vars=(0,1),linealpha=0.4)

we see that we can extend the amount of time until we deviate strongly from the "true" solution. Of course, for a chaotic system like the Lorenz one presented here, it is impossible to follow the true solution for long times, due to the fact that the system is chaotic and unavoidable deviations due to the numerical precision of a computer get amplified exponentially.

However, not all hope is lost. The shadowing theorem is a strong statement for having confidence in numerical evolution of chaotic systems:

Although a numerically computed chaotic trajectory diverges exponentially from the true trajectory with the same initial coordinates, there exists an errorless trajectory with a slightly different initial condition that stays near ("shadows") the numerically computed one.

For more info on the shadowing theorem, please see the book Chaos in Dynamical Systems by E. Ott.