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.
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 using the @ode_def
macro:
fitz = @ode_def_nohes FitzhughNagumo begin
dV = c*(V - V^3/3 + R)
dR = -(1/c)*(V - a - b*R)
end a=0.2 b=0.2 c=3.0
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
prob = ODEProblem(fitz,u0,tspan)
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 MonteCarloProblem for 100 trajectories:
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Euler(),num_monte=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)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Euler(),num_monte=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)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Euler(),num_monte=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())
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=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())
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=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:
g = @ode_def_bare LorenzExample begin
dx = σ*(y-x)
dy = x*(ρ-z) - y
dz = x*y - β*z
end σ=>10.0 ρ=>28.0 β=(8/3)
u0 = [1.0;0.0;0.0]
tspan = (0.0,30.0)
prob = ODEProblem(g,u0,tspan)
and then we build the ProbInts
type. Let's use the order 5 Tsit5
again.
cb = AdaptiveProbIntsUncertainty(5)
Then we solve the MonteCarloProblem
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Tsit5(),num_monte=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 junk. 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)
cb = AdaptiveProbIntsUncertainty(7)
monte_prob = MonteCarloProblem(prob)
sim = solve(monte_prob,Vern7(),num_monte=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 recieve junk.