# [show and tell] Blog post on interval arithmetic in Pluto.jl

I just wrote a blog post entirely in Julia and Pluto.jl! It’s the story of how I flew too close to the sun implementing my own interval arithmetic library, and lessons learned along the way. It’s got some math, some twists, love interests, action. A blog post for the whole family.

This is the first blog post I’ve ever published since my prehistoric wordpress days, and I have to say Pluto.jl was a dream to work with! It’s so good that it’s making me start to consider Julia for projects that I never would’ve otherwise done in Julia or in notebooks for that matter. It doesn’t have any of the same is-this-cell-stale anxiety that I had with Jupyter. Shout out to @fonsp and team!

27 Likes

Nice. Very clean. How did you host it? Is this more than just slapping the generated HTML onto a static site?

Thanks! I just exported to HTML and threw it onto a static site! I don’t know if this is the best solution, but it works for now…

5 Likes

It looks great. I guess I’m just amazed that it’s that easy to make that good looking of a site with Pluto. It makes it look like you did a ton of work, putting a working Binder integration button on the top and everything

3 Likes

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

7 Likes

Nice! What’s the trick here? What does @taylorize do?

Do you mean the math? You could check works by Berz and Makino, M. Joldes, X. Chen, F. Bünger to name just a few.

About the implementation, you could check a JuliaCon 2018 presentation by @lbenet and @dpsanders about Taylor models. The methods used above to propagate solutions of nonlinear ODEs is more recent though.

It is an optimization strategy used in TaylorIntegration.jl; I checked removing it in this example and there is no difference with and without it.

2 Likes

Great post!

For interactivity, you could use GitHub - JuliaPluto/PlutoSliderServer.jl: Web server to run just the @bind parts of a Pluto.jl notebook in addition. With it, users can manipulate PlutoUI elements on the web page itself (without running Pluto locally) and the shown notebook content updates in real-time.
The disadvantage is that it needs som (virtual) machine to run the server.

PlutoSliderServer is used e.g. here:

https://computationalthinking.mit.edu/Spring21/images/

Thanks! Yeah, I would’ve liked to add some interactivity but unfortunately the site is fully static, so I don’t have a VM to spare. It would really be nice if Julia could run in the browser via WebAssembly or something instead.

Wasm support for Julia would be great!
There are some efforts in this direction: GitHub - Keno/julia-wasm: Running julia on wasm, Web Platform Projects – Summer of Code

2 Likes