Download the notebook for this tutorial.
Following the book Validated Numerics (Princeton, 2011) by Warwick Tucker, we find rigorous (i.e., guaranteed, or validated) bounds on using interval arithmetic, via the
IntervalArithmetic.jl package. First, let's import the package
There are many ways to calculate . For illustrative purposes, we will use the following sum
It is known that the exact value is . Thus, if we can calculate rigorous bounds on , then we can find rigorous bounds on .
The idea is to split up into two parts, , with and .
By bounding using integrals from below and above, we can see that . Rigorous bounds on are found by calculating it using interval arithmetic.
may be calculated by summing either forwards or backwards. A naive (non-interval) version could be the following:
function forward_sum_naive(N) S_N = 0.0 for i in 1:N S_N += 1. / (i^2) end S_N end S = forward_sum_naive(10000) err = abs(S - pi^2/6.) # error S, err
(S, err) = (1.6448340718480652, 9.999500016122376e-5)
To find rigorous bounds for , we use interval arithmetic: each term is enclosed in an interval that is guaranteed to contain the true real value. A first idea is simply to wrap each term in the @interval macro, which converts its arguments into containing intervals:
function forward_S_N(N) S_N = (0.0) for i in 1:N S_N += ( 1. /(i^2) ) end S_N end N = 10^5 rigorous_approx_S_N = forward_S_N(N)
2.710170 seconds (16.48 M allocations: 705.208 MiB, 5.64% gc time) [1.64492, 1.64493]
We incorporate the respective bound on to obtain the bounds on , and hence on . We can also optimize the code by creating the interval
@interval(1) only once:
function forward_sum(N) S_N = (0.0) interval_one = (1.) for i in 1:N S_N += interval_one / (i^2) end T_N = interval_one / (N, N+1) S = S_N + T_N sqrt(6*S) end N = 10^6 S = forward_sum(N) S, diam(S)
0.657669 seconds (5.45 M allocations: 178.428 MiB, 7.03% gc time) (S, diam(S)) = ([3.14159, 3.1416], 2.1298829366855898e-10)
We can ask for the midpoint–radius representation, which shows that the calculated bounds are correct to around 10 decimal places:
We may check that the true value of is indeed contained in the interval:
big(π) ∈ S
We may repeat the calculation, but now summing in the opposite direction. Due to the way that floating-point arithmetic works, this gives a more accurate answer.
function reverse_sum(N) S_N = (0.0) interval_one = (1.) for i in N:-1:1 S_N += interval_one / (i^2) end T_N = interval_one / (N, N+1) S = S_N + T_N sqrt(6*S) end N = 10^6 S = reverse_sum(N) S, diam(S)
0.169420 seconds (3.03 M allocations: 77.610 MiB, 11.94% gc time) (S, diam(S)) = ([3.14159, 3.1416], 9.57456336436735e-13)
Note that the sqrt function is guaranteed (by the IEEE 754 standard) to give correctly-rounded results, so the resulting bounds on are guaranteed to be correct.
Note also that due to the way that the IntervalArithmetic package works, we can make the code simpler, at the expense of some performance. Only two explicit calls to the
@interval macro are now required:
function reverse_sum2(N) S_N = (0.0) for i in N:-1:1 S_N += 1 / (i^2) end T_N = 1 / (N, N+1) S = S_N + T_N sqrt(6*S) end N = 10^6 S = reverse_sum2(N) S, diam(S)
0.137809 seconds (2.06 M allocations: 48.472 MiB, 15.39% gc time) (S, diam(S)) = ([3.14159, 3.1416], 9.57456336436735e-13)