Nice post!
Let me share another example to further illustrate your statement that:
Practically speaking it means that although interval arithmetic generally works very well without supervision, you may want to write hand-crafted interval arithmetic functions for certain operations or consider different possible rewrites of your computation
This small example mixes intervals, the sine function… and differential equations, for the whole Julia 
familly. (Credits to @AnderGray for this nice example).
The function x(t) = x_0 \exp\left(\cos(t) - 1\right), x_0 \in \mathbb{R}, t \geq 0 solves the differential equation x'(t) = -x(t) ~ \sin(t).
Standard integration schemes fail to produce helpful solutions if the initial state is an interval, for reasons explained in your post, eg. dependency problems. This is easy to check:
using DifferentialEquations, IntervalArithmetic
# initial condition
x₀ = [-1 .. 1]
# define the problem
function f(dx, x, p, t)
dx[1] = -x[1] * sin(t)
end
# pass to solvers
prob = ODEProblem(f, x₀, (0.0, 2.0))
sol = solve(prob, Tsit5(), adaptive=false, dt=0.05, reltol=1e-6);
# there is no plot recipe readily available so we create a plot with the help of LazySets.jl
using LazySets, Plots
using LazySets: Interval
out = [Interval(sol.t[i]) × Interval(sol.u[i][1]) for i in 1:20]
plot(out, vars=(0, 1), xlab="t", ylab="x(t)", lw=3.0, alpha=1., c=:black, marker=:none, title="Standard integrator with interval initial condition")
See that the solution explodes already for 1 second. On the other hand, specialized algorithms can handle this case without noticeable wrapping effect, producing a sequence of sets whose union covers the true solution for all initial points.
using ReachabilityAnalysis, Plots
# define the model
@taylorize function f(dx, x, p, t)
dx[1] = -x[1] * sin(t)
end
# define the set of initial states; in this case, an interval
X0 = -1 .. 1
# define the initial-value problem
prob = @ivp(x' = f(x), x(0) ∈ X0, dim=1);
# solve it
sol = solve(prob, alg=TMJets21a(abstol=1e-10), T=15);
# it is illustrative to plot the computed flowpipe and the known analytic solution
# for a range of initial conditions that cover the range of `U0`.
analytic_sol(x0) = t -> x0 * exp(cos(t) - 1.0)
dt = range(0, 15, length=100)
x0vals = range(-1, 1, length=25)
fig = plot()
plot!(fig, sol, vars=(0, 1), lw=0.0, title="Specialized (Taylor-model based) integrator")
[plot!(fig, dt, analytic_sol(x0).(dt), lab="", c=:magenta, lw=2.0, xlab="t", ylab="x(t)") for x0 in x0vals]; #!jl
fig
We will discuss this and other examples at our JuliaCon workshop “It’s all set: a hands-on introduction to JuliaReach” next monday 
